vrt_warp.md

This url has a netcdf:

"https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc"

it has correct georeferencing (bbox 0, -90, 360, 90) but doesn't have the crs set ("EPSG:4326" is appropriate)

it also has multiple variables, which traditionally requires this verbose "subdataset syntax" for GDAL:

"NETCDF:\"/vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc\":sst"

enter vrt:// protocol

we can augment this with "dynamic" and very lightweight GDAL vrt, to ensure the CRS is set and we only are interested in "sst":

url = "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc"
dsn = f"vrt:///vsicurl/{url}?sd_name=sst&a_srs=EPSG:4326"

dsn
'vrt:///vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc?sd_name=sst&a_srs=EPSG:4326'

This works in many contexts.

dsn = 'vrt:///vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc?sd_name=sst&a_srs=EPSG:4326'


from osgeo import gdal
gdal.UseExceptions()
ds = gdal.Open(dsn)
ds.GetSpatialRef().ExportToWkt()
#'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"]]'

We can use the warp api to resolve to -180,180 and crop in the same step, with a temporary file we don't have to manage.

dsn = 'vrt:///vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc?sd_name=sst&a_srs=EPSG:4326'


from osgeo import gdal

gdal.UseExceptions()
ds = gdal.Open(dsn)
lon180ds = gdal.Warp("/vsimem/lon180.tif", ds, outputBounds = [-20, -30, 40, 25])
lon180ds.RasterXSize
#240
lon180ds.RasterYSize
#220
lon180ds.GetGeoTransform()
#(-20.0, 0.25, 0.0, 25.0, 0.0, -0.25)

This works the same in other languages (though, vsicurl+netcdf cannot work in windows or macos sadly).

dsn <- 'vrt:///vsicurl/https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/198110/oisst-avhrr-v02r01.19811002.nc?sd_name=sst&a_srs=EPSG:4326'

library(gdalraster)
warp(dsn, "/vsimem/lon180.tif", "", cl_arg = c("-te", -20, -30, 40, 25))
ds <- new(GDALRaster, "/vsimem/lon180.tif")
plot_raster(ds)
maps::map(add = TRUE)