Lazy resampling

Example on how to lazily resample a xarray.Dataset or xarray.DataArray with dask using odc.geo.xr.xr_reproject

import sys
import numpy as np
import xarray as xr
from rasterio.enums import Resampling
import odc.geo
from odc.geo.xr import xr_reproject
from odc.geo.geobox import GeoBox, BoundingBox

print(sys.version)
print(xr.__version__)
print(odc.geo.__version__)

3.13.1 | packaged by conda-forge | (main, Dec  5 2024, 21:23:54) [GCC 13.3.0]
2025.3.1
0.4.10

The following function is only to simluate that the crs, while initially present on the opened dataset, can get lost easily along the way in workflows, leading to the reprojection to fail. This is not a shortcoming of xarray or xr_reproject, but I want to document it, as this happened to me several times before.

def replace_2_by_0(da: xr.DataArray) -> xr.DataArray:
    """Replace 2s by 0"""
    return da.where(da != 2, 0)

ds_webcam = xr.open_zarr("./data/webcam_snow_cover.zarr")

ds_webcam["nodata"] = np.nan

Check that important info on crs is present

assert ds_webcam.odc.geobox

assert ds_webcam.odc.crs

Construct geobox needed for xr_reproject

xmin, xmax = 656320.0, 661020.0
ymax, ymin = 5213480.0, 5209240.0

lower_res = 20 # target resolution in meters

epsg_code = 25832

bbox = BoundingBox(left=xmin, top=ymax, bottom=ymin, right=xmax)
geobox = GeoBox.from_bbox(bbox=bbox, crs=epsg_code, resolution=lower_res)

Perform computation on dataset

I sometimes use xr.DataArray objects in computations. It is easy to oversee that the data array (here ds_webcam.snow_cover) does not inherit or carry over all the important information (.odc.crs and .odc.geobox) present on the parent dataset object (ds_webcam) when used in computations.

# perform the dummy operation
snow_cover = replace_2_by_0(ds_webcam.snow_cover)

# geobox is present!
assert snow_cover.odc.geobox

/home/<user>/miniconda3/envs/satpy/lib/python3.13/site-packages/odc/geo/_xr_interop.py:503: UserWarning: grid_mapping=spatial_ref is not pointing to valid coordinate
  warnings.warn(

# but the crs is not (as the above warning indicates)
snow_cover.odc.crs

To me this was initially counter intuitive, as the parent dataset has the complete information:

ds_webcam.odc.crs

CRS('PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]]')

Let’s fix that:

snow_cover = snow_cover.odc.assign_crs("epsg:25832")

Check source no data value and set correctly

# show that nodata attribute is missing
assert getattr(snow_cover, "nodata", None) is None

# add nodata attribute
snow_cover["nodata"] = np.nan

Perform lazy resampling

snow_cover_downsampled = xr_reproject(snow_cover, geobox, resampling=Resampling.mode)
snow_cover_downsampled.name = "snow_cover"

Which returns a dask-backed dataarray:

snow_cover_downsampled

<xarray.DataArray 'snow_cover' (time: 1, y: 212, x: 235)> Size: 399kB
dask.array<reproject, shape=(1, 212, 235), dtype=float64, chunksize=(1, 212, 235), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 8B 2021-07-15
    nodata       float64 8B nan
  * y            (y) float64 2kB 5.213e+06 5.213e+06 ... 5.209e+06 5.209e+06
  * x            (x) float64 2kB 6.563e+05 6.564e+05 ... 6.61e+05 6.61e+05
    spatial_ref  int32 4B 25832
Attributes:
    AREA_OR_POINT:  Area

xarray.DataArray

'snow_cover'

• time: 1
• y: 212
• x: 235

• dask.array<chunksize=(1, 212, 235), meta=np.ndarray>

  Array Chunk Bytes 389.22 kiB 389.22 kiB Shape (1, 212, 235) (1, 212, 235) Dask graph 1 chunks in 5 graph layers Data type float64 numpy.ndarray

• Coordinates: (5)
  • time
    (time)
    datetime64[ns]
    2021-07-15

    array(['2021-07-15T00:00:00.000000000'], dtype='datetime64[ns]')

  • nodata
    ()
    float64
    nan

    array(nan)

  • y
    (y)
    float64
    5.213e+06 5.213e+06 ... 5.209e+06

    units :

    metre

    resolution :

    -20.0

    crs :

    EPSG:25832

    array([5213470., 5213450., 5213430., ..., 5209290., 5209270., 5209250.],
          shape=(212,))

  • x
    (x)
    float64
    6.563e+05 6.564e+05 ... 6.61e+05

    units :

    metre

    resolution :

    20.0

    crs :

    EPSG:25832

    array([656330., 656350., 656370., ..., 660970., 660990., 661010.], shape=(235,))

  • spatial_ref
    ()
    int32
    25832

    spatial_ref :

    PROJCRS["ETRS89 / UTM zone 32N",BASEGEOGCRS["ETRS89",DATUM["European Terrestrial Reference System 1989",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4258]],CONVERSION["UTM zone 32N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Europe between 6°E and 12°E: Austria; Denmark - onshore and offshore; Germany - onshore and offshore; Italy - onshore and offshore; Norway including Svalbard - onshore and offshore; Spain - offshore."],BBOX[36.53,6,84.01,12.01]],USAGE[SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],AREA["Europe between 6°E and 12°E and approximately 36°30'N to 84°N."],BBOX[36.53,6,84.01,12.01]],ID["EPSG",25832]]

    crs_wkt :

    PROJCRS["ETRS89 / UTM zone 32N",BASEGEOGCRS["ETRS89",DATUM["European Terrestrial Reference System 1989",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4258]],CONVERSION["UTM zone 32N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Europe between 6°E and 12°E: Austria; Denmark - onshore and offshore; Germany - onshore and offshore; Italy - onshore and offshore; Norway including Svalbard - onshore and offshore; Spain - offshore."],BBOX[36.53,6,84.01,12.01]],USAGE[SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],AREA["Europe between 6°E and 12°E and approximately 36°30'N to 84°N."],BBOX[36.53,6,84.01,12.01]],ID["EPSG",25832]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    ETRS89

    horizontal_datum_name :

    European Terrestrial Reference System 1989

    projected_crs_name :

    ETRS89 / UTM zone 32N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    9.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    GeoTransform :

    656320 20 0 5213480 0 -20

    array(25832, dtype=int32)

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2021-07-15'], dtype='datetime64[ns]', name='time', freq=None))

  • y
    PandasIndex

    PandasIndex(Index([5213470.0, 5213450.0, 5213430.0, 5213410.0, 5213390.0, 5213370.0,
           5213350.0, 5213330.0, 5213310.0, 5213290.0,
           ...
           5209430.0, 5209410.0, 5209390.0, 5209370.0, 5209350.0, 5209330.0,
           5209310.0, 5209290.0, 5209270.0, 5209250.0],
          dtype='float64', name='y', length=212))

  • x
    PandasIndex

    PandasIndex(Index([656330.0, 656350.0, 656370.0, 656390.0, 656410.0, 656430.0, 656450.0,
           656470.0, 656490.0, 656510.0,
           ...
           660830.0, 660850.0, 660870.0, 660890.0, 660910.0, 660930.0, 660950.0,
           660970.0, 660990.0, 661010.0],
          dtype='float64', name='x', length=235))

• Attributes: (1)

  AREA_OR_POINT :

  Area

snow_cover_downsampled.plot()