gdalxarray-warp.md

vrt = gdalxarray.warp("/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif", 
     cutlineDSName =  "/vsicurl/https://github.com/mdsumner/geoboundaries/releases/download/latest/geoBoundariesCGAZ_ADM0.parquet",      
     cutlineSQL = "SELECT shapeGroup FROM geoBoundariesCGAZ_ADM0 WHERE shapeGroup IN ('AUS')",  
     cropToCutline = True, resolution = 0.1)
ds = xarray.open_dataset(vrt, engine = "gdalxarray", multidim = False)
ds
<xarray.Dataset> Size: 858kB
Dimensions:    (band: 1, y: 447, x: 947)
Coordinates:
  * band       (band) int64 8B 1
  * y          (y) float64 4kB -10.11 -10.21 -10.31 ... -54.51 -54.61 -54.71
  * x          (x) float64 8kB 73.29 73.39 73.49 73.59 ... 167.7 167.8 167.9
  * crs        int64 8B 0
Data variables:
    band_data  (band, y, x) int16 847kB ...
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
    crs      CRSIndex (crs=EPSG:4326)

composable from the above, we add two new arguments - crs and resolution - the output is a fully lazy xarray bound to GEBZCO 2024 via a vector crop contract that specifies a country, but the input, cutline, and output are entirely independent of each other in terms of crs, shape, bbox, resolution. (1653 m is what survives with nearest neighbour sampling,

 vrt = gdalxarray.warp("/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif", cutlineDSName =  "/vsicurl/https://github.com/mdsumner/geoboundaries/releases/download/latest/geoBoundariesCGAZ_ADM0.parquet", cutlineSQL = "SELECT shapeGroup FROM geoBoundariesCGAZ_ADM0 WHERE shapeGroup IN ('AUS')",  cropToCutline = True, resolution = 25000, crs = "+proj=laea")
ds = xarray.open_dataset(vrt, engine = "gdalxarray", multidim = False)
ds
<xarray.Dataset> Size: 275kB
Dimensions:    (band: 1, y: 383, x: 351)
Coordinates:
  * band       (band) int64 8B 1
  * y          (y) float64 3kB -1.903e+06 -1.928e+06 ... -1.143e+07 -1.145e+07
  * x          (x) float64 3kB 2.786e+06 2.811e+06 ... 1.151e+07 1.154e+07
  * crs        int64 8B 0
Data variables:
    band_data  (band, y, x) int16 269kB ...
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y
    crs      CRSIndex (crs=PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84 ...)
ds.band_data.values.max()
np.int16(1653)