With rasterio geometry_mask

Rasterix includes a dask-aware wrapper of rasterio.features.geometry_mask and rioxarray’s geometry_clip

Read data

Read in some raster data

import xarray as xr
import xproj  # noqa

ds = xr.tutorial.open_dataset("eraint_uvz")[["u"]]
ds = ds.proj.assign_crs(spatial_ref="epsg:4326")
ds

/home/docs/checkouts/readthedocs.org/user_builds/rasterix/envs/latest/lib/python3.12/site-packages/xarray/conventions.py:205: SerializationWarning: variable 'z' has non-conforming '_FillValue' np.float64(nan) defined, dropping '_FillValue' entirely.
  var = coder.decode(var, name=name)
/home/docs/checkouts/readthedocs.org/user_builds/rasterix/envs/latest/lib/python3.12/site-packages/xarray/conventions.py:205: SerializationWarning: variable 'u' has non-conforming '_FillValue' np.float64(nan) defined, dropping '_FillValue' entirely.
  var = coder.decode(var, name=name)
/home/docs/checkouts/readthedocs.org/user_builds/rasterix/envs/latest/lib/python3.12/site-packages/xarray/conventions.py:205: SerializationWarning: variable 'v' has non-conforming '_FillValue' np.float64(nan) defined, dropping '_FillValue' entirely.
  var = coder.decode(var, name=name)

<xarray.Dataset> Size: 6MB
Dimensions:      (month: 2, level: 3, latitude: 241, longitude: 480)
Coordinates:
  * month        (month) int32 8B 1 7
  * level        (level) int32 12B 200 500 850
  * latitude     (latitude) float32 964B 90.0 89.25 88.5 ... -88.5 -89.25 -90.0
  * longitude    (longitude) float32 2kB -180.0 -179.2 -178.5 ... 178.5 179.2
  * spatial_ref  int64 8B 0
Data variables:
    u            (month, level, latitude, longitude) float64 6MB ...
Indexes:
    spatial_ref  CRSIndex (crs=EPSG:4326)
Attributes:
    Conventions:  CF-1.0
    Info:         Monthly ERA-Interim data. Downloaded and edited by fabien.m...

xarray.Dataset

• Dimensions:
  • month: 2
  • level: 3
  • latitude: 241
  • longitude: 480
• Coordinates: (5)
  • month
    (month)
    int32
    1 7

    array([1, 7], dtype=int32)

  • level
    (level)
    int32
    200 500 850

    units :

    millibars

    long_name :

    pressure_level

    array([200, 500, 850], dtype=int32)

  • latitude
    (latitude)
    float32
    90.0 89.25 88.5 ... -89.25 -90.0

    units :

    degrees_north

    long_name :

    latitude

    array([ 90.  ,  89.25,  88.5 , ..., -88.5 , -89.25, -90.  ],
          shape=(241,), dtype=float32)

  • longitude
    (longitude)
    float32
    -180.0 -179.2 ... 178.5 179.2

    units :

    degrees_east

    long_name :

    longitude

    array([-180.  , -179.25, -178.5 , ...,  177.75,  178.5 ,  179.25],
          shape=(480,), dtype=float32)

  • spatial_ref
    ()
    int64
    0

    array(0)

• Data variables: (1)
  • u
    (month, level, latitude, longitude)
    float64
    ...

    number_of_significant_digits :

    2

    units :

    m s**-1

    long_name :

    U component of wind

    standard_name :

    eastward_wind

    [694080 values with dtype=float64]

• Indexes: (1)
  • spatial_ref
    CRSIndex (crs=EPSG:4326)

    CRSIndex
    <Geographic 2D CRS: EPSG:4326>
    Name: WGS 84
    Axis Info [ellipsoidal]:
    - Lat[north]: Geodetic latitude (degree)
    - Lon[east]: Geodetic longitude (degree)
    Area of Use:
    - name: World.
    - bounds: (-180.0, -90.0, 180.0, 90.0)
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

• Attributes: (2)

  Conventions :

  CF-1.0

  Info :

  Monthly ERA-Interim data. Downloaded and edited by fabien.maussion@uibk.ac.at

Read in example geometries

import geodatasets
import geopandas as gpd

world = gpd.read_file(geodatasets.get_path("naturalearth land"))
asia = world.iloc[slice(112, 113)]
asia

	featurecla	scalerank	min_zoom	geometry
112	Land	0	0.0	POLYGON ((106.97028 76.9743, 107.24011 76.4801...

from rasterix.rasterize import geometry_mask

geometry_mask(ds, asia[["geometry"]], xdim="longitude", ydim="latitude").plot()

<matplotlib.collections.QuadMesh at 0x7dc2af8a7a40>

By default, geometry_clip clips to the total_bounds of the provided geometries

from rasterix.rasterize import geometry_clip

n = geometry_clip(ds, asia[["geometry"]], xdim="longitude", ydim="latitude")
n.u.isel(month=1, level=0).plot()

<matplotlib.collections.QuadMesh at 0x7dc2af671f10>

Out-of-core support

All combinations of chunked and in-memory arrays and geometries are supported.

dask.array + geopandas

chunked = ds.chunk({"latitude": -1, "longitude": 120})
d = geometry_clip(chunked, asia[["geometry"]], xdim="longitude", ydim="latitude")
d

<xarray.Dataset> Size: 2MB
Dimensions:      (month: 2, level: 3, latitude: 150, longitude: 263)
Coordinates:
  * month        (month) int32 8B 1 7
  * level        (level) int32 12B 200 500 850
  * latitude     (latitude) float32 600B 77.25 76.5 75.75 ... -33.0 -33.75 -34.5
  * longitude    (longitude) float32 1kB -17.25 -16.5 -15.75 ... 178.5 179.2
  * spatial_ref  int64 8B 0
Data variables:
    u            (month, level, latitude, longitude) float64 2MB dask.array<chunksize=(2, 3, 150, 23), meta=np.ndarray>
Indexes:
    spatial_ref  CRSIndex (crs=EPSG:4326)
Attributes:
    Conventions:  CF-1.0
    Info:         Monthly ERA-Interim data. Downloaded and edited by fabien.m...

xarray.Dataset

• Dimensions:
  • month: 2
  • level: 3
  • latitude: 150
  • longitude: 263
• Coordinates: (5)
  • month
    (month)
    int32
    1 7

    array([1, 7], dtype=int32)

  • level
    (level)
    int32
    200 500 850

    units :

    millibars

    long_name :

    pressure_level

    array([200, 500, 850], dtype=int32)

  • latitude
    (latitude)
    float32
    77.25 76.5 75.75 ... -33.75 -34.5

    units :

    degrees_north

    long_name :

    latitude

    array([ 77.25,  76.5 ,  75.75,  75.  ,  74.25,  73.5 ,  72.75,  72.  ,  71.25,
            70.5 ,  69.75,  69.  ,  68.25,  67.5 ,  66.75,  66.  ,  65.25,  64.5 ,
            63.75,  63.  ,  62.25,  61.5 ,  60.75,  60.  ,  59.25,  58.5 ,  57.75,
            57.  ,  56.25,  55.5 ,  54.75,  54.  ,  53.25,  52.5 ,  51.75,  51.  ,
            50.25,  49.5 ,  48.75,  48.  ,  47.25,  46.5 ,  45.75,  45.  ,  44.25,
            43.5 ,  42.75,  42.  ,  41.25,  40.5 ,  39.75,  39.  ,  38.25,  37.5 ,
            36.75,  36.  ,  35.25,  34.5 ,  33.75,  33.  ,  32.25,  31.5 ,  30.75,
            30.  ,  29.25,  28.5 ,  27.75,  27.  ,  26.25,  25.5 ,  24.75,  24.  ,
            23.25,  22.5 ,  21.75,  21.  ,  20.25,  19.5 ,  18.75,  18.  ,  17.25,
            16.5 ,  15.75,  15.  ,  14.25,  13.5 ,  12.75,  12.  ,  11.25,  10.5 ,
             9.75,   9.  ,   8.25,   7.5 ,   6.75,   6.  ,   5.25,   4.5 ,   3.75,
             3.  ,   2.25,   1.5 ,   0.75,   0.  ,  -0.75,  -1.5 ,  -2.25,  -3.  ,
            -3.75,  -4.5 ,  -5.25,  -6.  ,  -6.75,  -7.5 ,  -8.25,  -9.  ,  -9.75,
           -10.5 , -11.25, -12.  , -12.75, -13.5 , -14.25, -15.  , -15.75, -16.5 ,
           -17.25, -18.  , -18.75, -19.5 , -20.25, -21.  , -21.75, -22.5 , -23.25,
           -24.  , -24.75, -25.5 , -26.25, -27.  , -27.75, -28.5 , -29.25, -30.  ,
           -30.75, -31.5 , -32.25, -33.  , -33.75, -34.5 ], dtype=float32)

  • longitude
    (longitude)
    float32
    -17.25 -16.5 -15.75 ... 178.5 179.2

    units :

    degrees_east

    long_name :

    longitude

    array([-17.25, -16.5 , -15.75, ..., 177.75, 178.5 , 179.25],
          shape=(263,), dtype=float32)

  • spatial_ref
    ()
    int64
    0

    array(0)

• Data variables: (1)
  • u
    (month, level, latitude, longitude)
    float64
    dask.array<chunksize=(2, 3, 150, 23), meta=np.ndarray>

    number_of_significant_digits :

    2

    units :

    m s**-1

    long_name :

    U component of wind

    standard_name :

    eastward_wind

    Array Chunk Bytes 1.81 MiB 843.75 kiB Shape (2, 3, 150, 263) (2, 3, 150, 120) Dask graph 3 chunks in 18 graph layers Data type float64 numpy.ndarray

• Indexes: (1)
  • spatial_ref
    CRSIndex (crs=EPSG:4326)

    CRSIndex
    <Geographic 2D CRS: EPSG:4326>
    Name: WGS 84
    Axis Info [ellipsoidal]:
    - Lat[north]: Geodetic latitude (degree)
    - Lon[east]: Geodetic longitude (degree)
    Area of Use:
    - name: World.
    - bounds: (-180.0, -90.0, 180.0, 90.0)
    Datum: World Geodetic System 1984 ensemble
    - Ellipsoid: WGS 84
    - Prime Meridian: Greenwich
    

• Attributes: (2)

  Conventions :

  CF-1.0

  Info :

  Monthly ERA-Interim data. Downloaded and edited by fabien.maussion@uibk.ac.at

xr.testing.assert_identical(n, d)

dask.array + dask-geopandas

import dask_geopandas as dgpd

dd = geometry_mask(
    ds.chunk({"latitude": -1, "longitude": 240}),
    dgpd.from_geopandas(world[["geometry"]], npartitions=3),
    xdim="longitude",
    ydim="latitude",
)
dd

<xarray.DataArray 'mask' (latitude: 241, longitude: 480)> Size: 116kB
dask.array<all-aggregate, shape=(241, 480), dtype=bool, chunksize=(241, 240), chunktype=numpy.ndarray>
Coordinates:
  * latitude     (latitude) float32 964B 90.0 89.25 88.5 ... -88.5 -89.25 -90.0
  * longitude    (longitude) float32 2kB -180.0 -179.2 -178.5 ... 178.5 179.2
    spatial_ref  int64 8B 0

xarray.DataArray

'mask'

• latitude: 241
• longitude: 480

• dask.array<chunksize=(241, 240), meta=np.ndarray>

  Array Chunk Bytes 112.97 kiB 56.48 kiB Shape (241, 480) (241, 240) Dask graph 2 chunks in 14 graph layers Data type bool numpy.ndarray

• Coordinates: (3)
  • latitude
    (latitude)
    float32
    90.0 89.25 88.5 ... -89.25 -90.0

    units :

    degrees_north

    long_name :

    latitude

    array([ 90.  ,  89.25,  88.5 , ..., -88.5 , -89.25, -90.  ],
          shape=(241,), dtype=float32)

  • longitude
    (longitude)
    float32
    -180.0 -179.2 ... 178.5 179.2

    units :

    degrees_east

    long_name :

    longitude

    array([-180.  , -179.25, -178.5 , ...,  177.75,  178.5 ,  179.25],
          shape=(480,), dtype=float32)

  • spatial_ref
    ()
    int64
    0

    array(0)

dd.plot()

<matplotlib.collections.QuadMesh at 0x7dc2a65933e0>