Skip to content

Instantly share code, notes, and snippets.

@mdsumner

mdsumner/erddap_vs_rangedownload.md

Last active July 3, 2025 20:31
Show Gist options
  • Save mdsumner/b157c366ec73376b8b6d7e915b8501bd to your computer and use it in GitHub Desktop.
Save mdsumner/b157c366ec73376b8b6d7e915b8501bd to your computer and use it in GitHub Desktop.
Download ZIP
Raw
erddap_vs_rangedownload.md

ERDDAP doesn't support range downloading

export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR

gdalinfo /vsicurl/https://coastwatch.pfeg.noaa.gov/erddap/files/jplMURSST41/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc
ERROR 1: 416: Range downloading not supported by this server: Error {
    code=416;
    message="Requested Range Not Satisfiable: REQUESTED_RANGE_NOT_SATISFIABLE: Don't try to connect to .nc or .hdf files on ERDDAP's /files/ system as if they were local files. It is horribly inefficient and often causes other problems. Instead: a) Use (OPeN)DAP client software to connect to ERDDAP's DAP services for this dataset (which have /griddap/ or /tabledap/ in the URL). That's what DAP is for. b) Or, use the dataset's Data Access Form to request a subset of data. c) Or, if you need the entire file or repeated access over a long period of time, use curl, wget, or your browser to download the entire file, then access the data from your local copy of the file.";
}

The files are actually sourced from here:

https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc

And for this we need EARTHDATA creds, I have my 'Authorization: Bearer ' token in ~/earthdata

export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR  ## this prevents GDAL scanning for sidecar files which we know we don't need
export GDAL_HTTP_HEADER_FILE=~/earthdata
gdalinfo vrt:///vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?sd_name=analysed_sst

I use vrt:// syntax to pick just one variable, that gives us this output:

Driver: VRT/Virtual Raster
Files: NETCDF:"/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc":analysed_sst
Size is 36000, 17999
Origin = (-179.995000000000005,89.994999999999990)
Pixel Size = (0.010000000000000,-0.010000000000000)
Metadata:
  analysed_sst#add_offset=298.15
  analysed_sst#comment=Interim near-real-time (nrt) version using Multi-Resolution Variational Analysis (MRVA) method for interpolation; to be replaced by Final version
  analysed_sst#coordinates=lon lat
  analysed_sst#long_name=analysed sea surface temperature
  analysed_sst#scale_factor=0.001
  analysed_sst#source=MODIS_T-JPL, MODIS_A-JPL, AMSR2-REMSS, AVHRRMTB_G-NAVO, iQUAM-NOAA/NESDIS, Ice_Conc-OSISAF
  analysed_sst#standard_name=sea_surface_foundation_temperature
  analysed_sst#units=kelvin
  analysed_sst#valid_max=32767
  analysed_sst#valid_min=-32767
  analysed_sst#_FillValue=-32768
  lat#axis=Y
  lat#comment=geolocations inherited from the input data without correction
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lat#valid_max=90
  lat#valid_min=-90
  lon#axis=X
  lon#comment=geolocations inherited from the input data without correction
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  lon#valid_max=180
  lon#valid_min=-180
  NC_GLOBAL#acknowledgment=Please acknowledge the use of these data with the following statement:  These data were provided by JPL under support by NASA MEaSUREs program.
  NC_GLOBAL#cdm_data_type=grid
  NC_GLOBAL#comment=Interim-MUR(nrt) will be replaced by MUR-Final in about 3 days; MUR = "Multi-scale Ultra-high Resolution"
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#[email protected]
  NC_GLOBAL#creator_name=JPL MUR SST project
  NC_GLOBAL#creator_url=http://mur.jpl.nasa.gov
  NC_GLOBAL#date_created=20250702T084701Z
  NC_GLOBAL#easternmost_longitude=180
  NC_GLOBAL#file_quality_level=3
  NC_GLOBAL#gds_version_id=2.0
  NC_GLOBAL#geospatial_lat_resolution=0.0099999998
  NC_GLOBAL#geospatial_lat_units=degrees north
  NC_GLOBAL#geospatial_lon_resolution=0.0099999998
  NC_GLOBAL#geospatial_lon_units=degrees east
  NC_GLOBAL#history=near real time (nrt) version created at nominal 1-day latency.
  NC_GLOBAL#id=MUR-JPL-L4-GLOB-v04.1
  NC_GLOBAL#institution=Jet Propulsion Laboratory
  NC_GLOBAL#keywords=Oceans > Ocean Temperature > Sea Surface Temperature
  NC_GLOBAL#keywords_vocabulary=NASA Global Change Master Directory (GCMD) Science Keywords
  NC_GLOBAL#license=These data are available free of charge under data policy of JPL PO.DAAC.
  NC_GLOBAL#Metadata_Conventions=Unidata Observation Dataset v1.0
  NC_GLOBAL#metadata_link=http://podaac.jpl.nasa.gov/ws/metadata/dataset/?format=iso&shortName=MUR-JPL-L4-GLOB-v04.1
  NC_GLOBAL#naming_authority=org.ghrsst
  NC_GLOBAL#netcdf_version_id=4.1
  NC_GLOBAL#northernmost_latitude=90
  NC_GLOBAL#platform=Terra, Aqua, GCOM-W, MetOp-B, Buoys/Ships
  NC_GLOBAL#processing_level=L4
  NC_GLOBAL#product_version=04.1nrt
  NC_GLOBAL#project=NASA Making Earth Science Data Records for Use in Research Environments (MEaSUREs) Program
  NC_GLOBAL#[email protected]
  NC_GLOBAL#publisher_name=GHRSST Project Office
  NC_GLOBAL#publisher_url=http://www.ghrsst.org
  NC_GLOBAL#references=http://podaac.jpl.nasa.gov/Multi-scale_Ultra-high_Resolution_MUR-SST
  NC_GLOBAL#sensor=MODIS, AMSR2, AVHRR, in-situ
  NC_GLOBAL#source=MODIS_T-JPL, MODIS_A-JPL, AMSR2-REMSS, AVHRRMTB_G-NAVO, iQUAM-NOAA/NESDIS, Ice_Conc-OSISAF
  NC_GLOBAL#southernmost_latitude=-90
  NC_GLOBAL#spatial_resolution=0.01 degrees
  NC_GLOBAL#standard_name_vocabulary=NetCDF Climate and Forecast (CF) Metadata Convention
  NC_GLOBAL#start_time=20250701T090000Z
  NC_GLOBAL#stop_time=20250701T090000Z
  NC_GLOBAL#summary=A merged, multi-sensor L4 Foundation SST analysis product from JPL.
  NC_GLOBAL#time_coverage_end=20250701T210000Z
  NC_GLOBAL#time_coverage_start=20250630T210000Z
  NC_GLOBAL#title=Daily MUR SST, Interim near-real-time (nrt) product
  NC_GLOBAL#uuid=27665bc0-d5fc-11e1-9b23-0800200c9a66
  NC_GLOBAL#westernmost_longitude=-180
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={1,4}
  NETCDF_DIM_time_VALUES=1404205200
  time#axis=T
  time#comment=Nominal time of analyzed fields
  time#long_name=reference time of sst field
  time#standard_name=time
  time#units=seconds since 1981-01-01 00:00:00 UTC
Geolocation:
  GEOREFERENCING_CONVENTION=PIXEL_CENTER
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=NETCDF:"/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc":lon
  Y_BAND=1
  Y_DATASET=NETCDF:"/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc":lat
Corner Coordinates:
Upper Left  (-179.9950000,  89.9950000) 
Lower Left  (-179.9950000, -89.9950000) 
Upper Right ( 180.0050000,  89.9950000) 
Lower Right ( 180.0050000, -89.9950000) 
Center      (   0.0050000,  -0.0000000) 
Band 1 Block=2047x1023 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Unit Type: kelvin
  Offset: 298.15,   Scale:0.001
  Metadata:
    add_offset=298.15
    comment=Interim near-real-time (nrt) version using Multi-Resolution Variational Analysis (MRVA) method for interpolation; to be replaced by Final version
    coordinates=lon lat
    long_name=analysed sea surface temperature
    NETCDF_DIM_time=1404205200
    NETCDF_VARNAME=analysed_sst
    scale_factor=0.001
    source=MODIS_T-JPL, MODIS_A-JPL, AMSR2-REMSS, AVHRRMTB_G-NAVO, iQUAM-NOAA/NESDIS, Ice_Conc-OSISAF
    standard_name=sea_surface_foundation_temperature
    units=kelvin
    valid_max=32767
    valid_min=-32767
    _FillValue=-32768

But, note how we have geolocation arrays (that means GDAL is referencing the longitude and latitude 1D arrays to determine the extent of the grid, and these are problematic 1d arrays, AND the girid has no CRS (the CRS that is there refers to the geolocation arrays).

There's no "geotransform", i.e. no "Origin" and "Pixel Size"

We want gdal to ignore the longitude,latitude arrays and use the correct extent instead so we add some more vrt:// args. We also add offset/scale so that the values come out scaled to Celsius rather than Kelvin.

Now we get geotransform and crs.

gdalinfo "vrt:///vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?sd_name=analysed_sst&a_ullr=-180,89.995,180,-89.995&a_offset=25&a_scale=0.001&a_srs=EPSG:4326"

Driver: VRT/Virtual Raster
Files: NETCDF:"/vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc":analysed_sst
Size is 36000, 17999
Coordinate System is:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,89.995000000000005)
Pixel Size = (0.010000000000000,-0.010000000000000)
Corner Coordinates:
Upper Left  (-180.0000000,  89.9950000) (180d 0' 0.00"W, 89d59'42.00"N)
Lower Left  (-180.0000000, -89.9950000) (180d 0' 0.00"W, 89d59'42.00"S)
Upper Right ( 180.0000000,  89.9950000) (180d 0' 0.00"E, 89d59'42.00"N)
Lower Right ( 180.0000000, -89.9950000) (180d 0' 0.00"E, 89d59'42.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=2047x1023 Type=Int16, ColorInterp=Undefined
  NoData Value=-32768
  Offset: 25,   Scale:0.001

Now that we have a fully sanitized NetCDF layer we can remap it easily, or read any subset we like.

## these could be set with bash+export or with terra::set_gdal_config
Sys.setenv("GDAL_DISABLE_READDIR_ON_OPEN" = "EMPTY_DIR")
Sys.setenv("GDAL_HTTP_HEADER_FILE" = "~/earthdata")

library(terra)
dsn <- "vrt:///vsicurl/https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?sd_name=analysed_sst&
a_ullr=-180,89.995,180,-89.995&a_offset=25&a_scale=0.001&a_srs=EPSG:4326"

(r <- rast(dsn))

crop(r, ext(140, 150, -50, -40))
class       : SpatRaster 
size        : 1000, 1000, 1  (nrow, ncol, nlyr)
resolution  : 0.01, 0.01  (x, y)
extent      : 140, 150, -50.005, -40.005  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
varname     : 20250701090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?sd_name=analysed_sst&a_ullr=-180,89.995,180,-89.995&a_offset=25&a_scale=0 
name        : 20250701090000-JPL-L4_GHRSST-S~,-89.995&a_offset=25&a_scale=0 
min value   :                                                         7.102 
max value   :                                                        16.515

image

Or, a different way to achieve similar is use a local projection.

rr <- project(r, rast(ext(c(-1, 1, -1, 1) * 1e6), res = 1000, crs = "+proj=laea +lon_0=145 +lat_0=-45"), by_util = TRUE)

## cos(lat) to get 10*1degree at this latitude (approx)
rr <- project(r, rast(ext(c(-1, 1, -1, 1) * 1e6 * cos(45 * pi/180)), res = 1000, crs = "+proj=laea +lon_0=145 +lat_0=-45"), by_util = TRUE)

image

@mdsumner
Copy link
Author

mdsumner commented Jul 3, 2025

Here's a bit on using stars for some of this example: https://nmfs-opensci.github.io/NMFSHackDays-2025/topics-2025/2025-02-21-earthdata/4-data-cubes.html

I don't understand how to use aggregate()

library(stars)
vrtarg <- "?sd_name=analysed_sst&a_ullr=-180,89.995,180,-89.995&a_offset=25&a_scale=0.001&a_srs=EPSG:4326"

results <- c("https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200116090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc", 
             "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200117090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc",
             "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200118090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc")

dsn <- sprintf("vrt:///vsicurl/%s%s", results, vrtarg)
time <- as.POSIXct(stringr::str_extract(basename(results), "^[0-9]{8}"), "%Y%m%d", tz = "UTC")
st <- stars::read_stars(dsn, along = "time")
attr(st, "dimensions")$time$values <- time
attr(st, "dimensions")$time$offset <- time[1]
attr(st, "dimensions")$time$delta <- diff(time[1:2])


stcrop <- sf::st_crop(st, sf::st_bbox(c(xmin=-75.5, xmax=-73.5,  ymin=33.5, ymax=35.5 ), crs = sf::st_crs(st)))
plot(stcrop)

image

@mdsumner
Copy link
Author

mdsumner commented Jul 3, 2025 •
edited
Loading

Use 'apply' rather than aggregate

st_as_stars(st_apply(stcrop, 1:2, mean))
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
stars object with 2 dimensions and 1 attribute
attribute(s):
          Min.  1st Qu.   Median     Mean 3rd Qu.     Max.
mean  16.99267 23.54833 24.01417 23.97449   24.84 25.61933
dimension(s):
   from    to offset delta refsys x/y
x 10451 10650   -180  0.01 WGS 84 [x]
y  5450  5650     90 -0.01 WGS 84 [y]

or crop and take temporal mean with terra

library(terra)
vrtarg <- "?sd_name=analysed_sst&a_ullr=-180,89.995,180,-89.995&a_offset=25&a_scale=0.001&a_srs=EPSG:4326"

results <- c("https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200116090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc",
             "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200117090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc",
             "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20200118090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc")

dsn <- sprintf("vrt:///vsicurl/%s%s", results, vrtarg)
mean(crop(rast(dsn), ext(c(xmin=-75.5, xmax=-73.5,  ymin=33.5, ymax=35.5 ))))


class       : SpatRaster
size        : 200, 200, 1  (nrow, ncol, nlyr)
resolution  : 0.01, 0.01  (x, y)
extent      : -75.5, -73.5, 33.495, 35.495  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s)   : memory
name        :     mean
min value   : 17.10200
max value   : 25.61933

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment