Demo on the “most common” regridder

For regridding categorical data to a coarser resolution, normal regridding methods don’t work. This can be data like land cover or vegetation type.

For this reason we implemented a “most common” regridder, which takes the most common value of the source data cells that overlap with the target grid cells.

If you want to regrid data to a finer resolution, the “nearest-neighbor” regridder can be used.

We start with importing the required packages:

import xarray as xr
import xarray_regrid

import dask.distributed

client = dask.distributed.Client()

Next we need a high resolution dataset to regrid. We used the LCCS land cover data which is available from the Climate Data Store.

Note the data is loaded in as a dask arrays, allowing for lazy computation.

ds = xr.open_dataset(
    "/data/C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc",
    chunks="auto",
)

da = ds["lccs_class"]  # Only take the class variable.
da = da.sortby(["lat", "lon"])
da = da.rename({"lat": "latitude", "lon": "longitude"})
da

<xarray.DataArray 'lccs_class' (time: 1, latitude: 64800, longitude: 129600)> Size: 8GB
dask.array<getitem, shape=(1, 64800, 129600), dtype=uint8, chunksize=(1, 9257, 10125), chunktype=numpy.ndarray>
Coordinates:
  * latitude   (latitude) float64 518kB -90.0 -90.0 -89.99 ... 89.99 90.0 90.0
  * longitude  (longitude) float64 1MB -180.0 -180.0 -180.0 ... 180.0 180.0
  * time       (time) datetime64[ns] 8B 2020-01-01
Attributes:
    standard_name:        land_cover_lccs
    flag_colors:          #ffff64 #ffff64 #ffff00 #aaf0f0 #dcf064 #c8c864 #00...
    long_name:            Land cover class defined in LCCS
    valid_min:            1
    valid_max:            220
    ancillary_variables:  processed_flag current_pixel_state observation_coun...
    flag_meanings:        no_data cropland_rainfed cropland_rainfed_herbaceou...
    flag_values:          [  0  10  11  12  20  30  40  50  60  61  62  70  7...

xarray.DataArray

'lccs_class'

• time: 1
• latitude: 64800
• longitude: 129600

• dask.array<chunksize=(1, 9257, 10125), meta=np.ndarray>

  Array Chunk Bytes 7.82 GiB 89.39 MiB Shape (1, 64800, 129600) (1, 9257, 10125) Dask graph 104 chunks in 4 graph layers Data type uint8 numpy.ndarray

• Coordinates: (3)
  • latitude
    (latitude)
    float64
    -90.0 -90.0 -89.99 ... 90.0 90.0

    units :

    degrees_north

    long_name :

    latitude

    standard_name :

    latitude

    valid_min :

    -90.0

    valid_max :

    90.0

    bounds :

    lat_bounds

    axis :

    Y

    array([-89.998611, -89.995833, -89.993056, ...,  89.993056,  89.995833,
            89.998611])

  • longitude
    (longitude)
    float64
    -180.0 -180.0 ... 180.0 180.0

    units :

    degrees_east

    long_name :

    longitude

    standard_name :

    longitude

    valid_min :

    -180.0

    valid_max :

    180.0

    bounds :

    lon_bounds

    axis :

    X

    array([-179.998611, -179.995833, -179.993056, ...,  179.993056,  179.995833,
            179.998611])

  • time
    (time)
    datetime64[ns]
    2020-01-01

    standard_name :

    time

    long_name :

    time

    axis :

    T

    bounds :

    time_bounds

    array(['2020-01-01T00:00:00.000000000'], dtype='datetime64[ns]')

• Indexes: (3)
  • latitude
    PandasIndex

    PandasIndex(Index([-89.99861111111112, -89.99583333333334, -89.99305555555556,
           -89.99027777777778,           -89.9875, -89.98472222222222,
           -89.98194444444445, -89.97916666666667, -89.97638888888889,
           -89.97361111111111,
           ...
            89.97361111111113,  89.97638888888889,  89.97916666666669,
            89.98194444444445,  89.98472222222222,  89.98750000000001,
            89.99027777777778,  89.99305555555557,  89.99583333333334,
            89.99861111111113],
          dtype='float64', name='latitude', length=64800))

  • longitude
    PandasIndex

    PandasIndex(Index([ -179.9986111111111, -179.99583333333334, -179.99305555555554,
           -179.99027777777778,           -179.9875, -179.98472222222222,
           -179.98194444444445, -179.97916666666666,  -179.9763888888889,
            -179.9736111111111,
           ...
             179.9736111111111,   179.9763888888889,  179.97916666666669,
            179.98194444444448,  179.98472222222222,            179.9875,
             179.9902777777778,  179.99305555555554,  179.99583333333334,
            179.99861111111113],
          dtype='float64', name='longitude', length=129600))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-01'], dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (8)

  standard_name :

  land_cover_lccs

  flag_colors :

  #ffff64 #ffff64 #ffff00 #aaf0f0 #dcf064 #c8c864 #006400 #00a000 #00a000 #aac800 #003c00 #003c00 #005000 #285000 #285000 #286400 #788200 #8ca000 #be9600 #966400 #966400 #966400 #ffb432 #ffdcd2 #ffebaf #ffc864 #ffd278 #ffebaf #00785a #009678 #00dc82 #c31400 #fff5d7 #dcdcdc #fff5d7 #0046c8 #ffffff

  long_name :

  Land cover class defined in LCCS

  valid_min :

  1

  valid_max :

  220

  ancillary_variables :

  processed_flag current_pixel_state observation_count change_count

  flag_meanings :

  no_data cropland_rainfed cropland_rainfed_herbaceous_cover cropland_rainfed_tree_or_shrub_cover cropland_irrigated mosaic_cropland mosaic_natural_vegetation tree_broadleaved_evergreen_closed_to_open tree_broadleaved_deciduous_closed_to_open tree_broadleaved_deciduous_closed tree_broadleaved_deciduous_open tree_needleleaved_evergreen_closed_to_open tree_needleleaved_evergreen_closed tree_needleleaved_evergreen_open tree_needleleaved_deciduous_closed_to_open tree_needleleaved_deciduous_closed tree_needleleaved_deciduous_open tree_mixed mosaic_tree_and_shrub mosaic_herbaceous shrubland shrubland_evergreen shrubland_deciduous grassland lichens_and_mosses sparse_vegetation sparse_tree sparse_shrub sparse_herbaceous tree_cover_flooded_fresh_or_brakish_water tree_cover_flooded_saline_water shrub_or_herbaceous_cover_flooded urban bare_areas bare_areas_consolidated bare_areas_unconsolidated water snow_and_ice

  flag_values :

  [ 0 10 11 12 20 30 40 50 60 61 62 70 71 72 80 81 82 90 100 110 120 121 122 130 140 150 151 152 153 160 170 180 190 200 201 202 210 220]

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

colors = da.attrs["flag_colors"].split(" ")
cmap = ListedColormap(colors)

ax = da.sel(latitude=slice(51, 54), longitude=slice(3.4, 6.4)).plot(cmap=cmap, vmin=10, vmax=220)
ax = plt.gca()
ax.set_aspect('equal')

We will also define our target grid:

from xarray_regrid import Grid

target_dataset = Grid(
    north=90,
    east=90,
    south=0,
    west=0,
    resolution_lat=1,
    resolution_lon=1,
).create_regridding_dataset()

The default chunks are a bit large for this regridding operation, so we need to rechunk before continuing to avoid memory issues:

da = da.chunk({"time": -1, "latitude": 4050, "longitude": 4050})
da.data

Array Chunk Bytes 7.82 GiB 15.64 MiB Shape (1, 64800, 129600) (1, 4050, 4050) Dask graph 512 chunks in 5 graph layers Data type uint8 numpy.ndarray

Using regrid.most_common you can now regrid the data. This is currently only implemented for DataArrays, not xr.Dataset.

Note that we have to provide the values of the expected labels in the data. This dataset already conventiently stores these in the attributes.

da_regrid = da.regrid.most_common(
    target_dataset, values=da.attrs["flag_values"], time_dim="time"
)

When we call .plot on the DataArray, computation will begin.

da_regrid.plot(x="longitude", cmap=cmap, vmin=10, vmax=220)

<matplotlib.collections.QuadMesh at 0x7fdec58fdaf0>