VRT warp for CERSAT (augment NetCDFs that are in a projection)

[mdsumner]

What follows is a note-to-self, about how we could streamline the fudge-fixes in raadtools as formal parameter sets, so we can still have our handy raad functions but also analogous virtual raster definitions that we can stream through the raster readers (and GDALwarp, QGIS, Manifold, all the same).

This VRT text, created for a specific subdataset for a daily file:

SDS_datasource <- 'NETCDF:"/rdsi/PUBLIC/raad/data/ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-concentration/data/antarctic/daily/netcdf/2021/20210317.nc":concentration'

sprintf('
<VRTDataset rasterXSize="632" rasterYSize="664" subClass="VRTWarpedDataset">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / Antarctic Polar Stereographic",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-71],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",NORTH],AXIS["Northing",NORTH],AUTHORITY["EPSG","3031"]]</SRS>
  <GeoTransform> -3.9500000000000000e+06,  1.2500000000000002e+04,  0.0000000000000000e+00,  4.3500000000000000e+06,  0.0000000000000000e+00, -1.2500000000000002e+04</GeoTransform>
  <Metadata>
    <MDI key="concentration#add_offset">0</MDI>
    <MDI key="concentration#long_name">sea-ice concentration</MDI>
    <MDI key="concentration#missing_value">-128</MDI>
    <MDI key="concentration#scale_factor">1</MDI>
    <MDI key="concentration#units">percent</MDI>
    <MDI key="concentration#_FillValue">-128</MDI>
    <MDI key="NC_GLOBAL#CONVENTIONS">COARDS</MDI>
    <MDI key="NC_GLOBAL#creation_time">2007-016T10:51:07.000</MDI>
    <MDI key="NC_GLOBAL#grid">NSIDC</MDI>
    <MDI key="NC_GLOBAL#instrument">SSM/I</MDI>
    <MDI key="NC_GLOBAL#long_name">Sea-ice concentration as observed by SSM/I</MDI>
    <MDI key="NC_GLOBAL#netcdf_version_id">3.4</MDI>
    <MDI key="NC_GLOBAL#platform_id">F13</MDI>
    <MDI key="NC_GLOBAL#pole">south</MDI>
    <MDI key="NC_GLOBAL#producer_agency">IFREMER</MDI>
    <MDI key="NC_GLOBAL#producer_institution">CERSAT</MDI>
    <MDI key="NC_GLOBAL#product_version">2.0</MDI>
    <MDI key="NC_GLOBAL#short_name">PSI-F13-Concentration</MDI>
    <MDI key="NC_GLOBAL#spatial_resolution">12.5 km</MDI>
    <MDI key="NC_GLOBAL#time_resolution">daily</MDI>
    <MDI key="NETCDF_DIM_EXTRA">{time}</MDI>
    <MDI key="NETCDF_DIM_time_DEF">{1,4}</MDI>
    <MDI key="NETCDF_DIM_time_VALUES">805764</MDI>
    <MDI key="time#long_name">time</MDI>
    <MDI key="time#units">hours since 1900-1-1 0:0:0</MDI>
  </Metadata>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
    <Metadata>
      <MDI key="add_offset">0</MDI>
      <MDI key="long_name">sea-ice concentration</MDI>
      <MDI key="missing_value">-128</MDI>
      <MDI key="NETCDF_DIM_time">805764</MDI>
      <MDI key="NETCDF_VARNAME">concentration</MDI>
      <MDI key="scale_factor">1</MDI>
      <MDI key="units">percent</MDI>
      <MDI key="_FillValue">-128</MDI>
    </Metadata>
    <NoDataValue>-128</NoDataValue>
    <UnitType>percent</UnitType>
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>Bilinear</ResampleAlg>
    <WorkingDataType>Float32</WorkingDataType>
    <Option name="INIT_DEST">NO_DATA</Option>
    <Option name="ERROR_OUT_IF_EMPTY_SOURCE_WINDOW">FALSE</Option>
    <SourceDataset relativeToVRT="0">%s</SourceDataset>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoTransform>-3950000,12500,0,-3950000,0,12500</SrcGeoTransform>
            <SrcInvGeoTransform>316,8.00000000000000065e-05,0,316,0,8.00000000000000065e-05</SrcInvGeoTransform>
            <DstGeoTransform>-3950000,12500.0000000000018,0,4350000,0,-12500.0000000000018</DstGeoTransform>
            <DstInvGeoTransform>315.999999999999943,7.9999999999999993e-05,0,347.999999999999943,0,-7.9999999999999993e-05</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1">
        <SrcNoDataReal>-128</SrcNoDataReal>
        <SrcNoDataImag>0</SrcNoDataImag>
        <DstNoDataReal>-128</DstNoDataReal>
        <DstNoDataImag>0</DstNoDataImag>
      </BandMapping>
    </BandList>
  </GDALWarpOptions>
</VRTDataset>
', SDS_datasource)

resolves to a fully-fledged (lazy) raster object:

class      : RasterLayer 
dimensions : 664, 632, 419648  (nrow, ncol, ncell)
resolution : 12500, 12500  (x, y)
extent     : -3950000, 3950000, -3950000, 4350000  (xmin, xmax, ymin, ymax)
crs        : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
source     : .../geo.vrt 
names      : geo 
values     : 0, 255  (min, max)

created by

• generating VRT from the source NetCDF with gdal_translate
• modifying the GeoTransform: negating the y-scale and setting ymax to the ymin (because the data is flipped in y-positive
  convention)
• gdalwarp that VRT to a new vrt (and now we can forget the first one), replace the 'SourceDataset' tag with the original subdataset string (not the temp .vrt)

The flip is encoded in the Src and Dst GeoTransforms. So we can

• generate this warping VRT on the fly (with our known geotransform-flip values and CRS)
• stream that through raster/terra/stars or GDALwarp directly for any given target grid.

Notes

• don't think there's any file specific details in the VRT beyond SourceDataset (no date etc), so we could add that for the date possibly
• not using 'glue' because there's other curly braces and I don't know about escaping those

[mdsumner]

couple of VRTs and a rough nsidc example using it

vrt.zip

library(raster)

prj <- "+proj=laea +lon_0=175 +lat_0=-72"
r0 <- projectExtent(raster::raster(raster::extent(100, 250, -90, -50), crs = "+proj=longlat", res = 0.2), prj)
#res(r0) <- res(r0)*4

#e <- new("Extent", xmin = -729574.370995987, xmax = 1006027.80116313, 
#         ymin = -812611.063431228, ymax = 765209.093077062)
#r0 <- raster(e, res = 30000, crs = prj)

files <- raadtools::icefiles()


gt <- affinity::extent_dim_to_gt(raster::extent(r0), dim(r0)[2:1])
wkt <- sf::st_crs(projection(r0))$wkt


tfile <- tempfile(fileext = ".vrt")
file <- sample(files$fullname, 1)
writeLines(sprintf(readLines("~/vrt0/nsidc_south.vrt"), file), tfile)

r <- raster::setValues(r0, vapour:::vapour_warp_raster(tfile, geotransform = gt, dimension = dim(r0)[2:1], wkt = wkt)[[1]]/2.5)
nobig <- function(x) x[x < 101]
brks <- sort(unique(quantile(nobig(values(r)), seq(0, 1, length.out = 24))))
plot(r, col = grey.colors(length(brks) - 1), breaks = brks)