With rasterio rasterize

Rasterix includes a dask-aware wrapper for rasterio.features.rasterize.

Read data

Read in some raster data

import xarray as xr
import xproj  # noqa

ds = xr.tutorial.open_dataset("eraint_uvz")
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: 17MB
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:
    z            (month, level, latitude, longitude) float64 6MB ...
    u            (month, level, latitude, longitude) float64 6MB ...
    v            (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: (3)
  • z
    (month, level, latitude, longitude)
    float64
    ...

    number_of_significant_digits :

    5

    units :

    m**2 s**-2

    long_name :

    Geopotential

    standard_name :

    geopotential

    [694080 values with dtype=float64]

  • 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]

  • v
    (month, level, latitude, longitude)
    float64
    ...

    number_of_significant_digits :

    2

    units :

    m s**-1

    long_name :

    V component of wind

    standard_name :

    northward_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"))
world

	featurecla	scalerank	min_zoom	geometry
0	Land	1	1.0	POLYGON ((-59.57209 -80.04018, -59.86585 -80.5...
1	Land	1	1.0	POLYGON ((-159.20818 -79.49706, -161.1276 -79....
2	Land	1	0.0	POLYGON ((-45.15476 -78.04707, -43.92083 -78.4...
3	Land	1	1.0	POLYGON ((-121.21151 -73.50099, -119.91885 -73...
4	Land	1	1.0	POLYGON ((-125.55957 -73.48135, -124.03188 -73...
...	...	...	...	...
122	Land	1	1.0	POLYGON ((51.13619 80.54728, 49.79368 80.41543...
123	Land	0	0.0	POLYGON ((99.93976 78.88094, 97.75794 78.7562,...
124	Land	0	0.0	POLYGON ((-87.02 79.66, -85.81435 79.3369, -87...
125	Land	0	0.0	POLYGON ((-68.5 83.10632, -65.82735 83.02801, ...
126	Land	0	0.0	POLYGON ((-27.10046 83.51966, -20.84539 82.726...

127 rows × 4 columns

from rasterix.rasterize import rasterize

n = rasterize(ds, world[["geometry"]], xdim="longitude", ydim="latitude", engine="rasterio")
n.plot()

<matplotlib.collections.QuadMesh at 0x784c44b47620>

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 = rasterize(chunked, world[["geometry"]], xdim="longitude", ydim="latitude", engine="rasterio")
d

<xarray.DataArray 'rasterized' (latitude: 241, longitude: 480)> Size: 116kB
dask.array<replace_values, shape=(241, 480), dtype=uint8, chunksize=(241, 120), 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

'rasterized'

• latitude: 241
• longitude: 480

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

  Array Chunk Bytes 112.97 kiB 28.24 kiB Shape (241, 480) (241, 120) Dask graph 4 chunks in 16 graph layers Data type uint8 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)

xr.testing.assert_identical(n, d)

dask.array + dask-geopandas

import dask_geopandas as dgpd

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

<xarray.DataArray 'rasterized' (latitude: 241, longitude: 480)> Size: 116kB
dask.array<replace_values, shape=(241, 480), dtype=uint8, 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

'rasterized'

• 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 17 graph layers Data type uint8 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)

xr.testing.assert_identical(dd, d)