Spatial Regridding and Masking

Spatial Regridding and masking data cubes¶

A DeepESDL example notebook¶

This notebook demonstrates how to access two different data sets, transform one of them so that both share the same spatial grid, and mask one dataset by applying a condition from the other dataset.

Please, also refer to the DeepESDL documentation and visit the platform's website for further information!

Brockmann Consult, 2024

This notebook runs with the python environment deepesdl-xcube-1.7.0, please checkout the documentation for help on changing the environment.

In [1]:

Copied!





import datetime

import numpy as np
import pandas as pd
import shapely.geometry
import xarray as xr
from IPython.display import JSON
from xcube.core.gridmapping import GridMapping
from xcube.core.maskset import MaskSet
from xcube.core.resampling import resample_in_space
from xcube.core.store import new_data_store
import datetime

import numpy as np
import pandas as pd
import shapely.geometry
import xarray as xr
from IPython.display import JSON
from xcube.core.gridmapping import GridMapping
from xcube.core.maskset import MaskSet
from xcube.core.resampling import resample_in_space
from xcube.core.store import new_data_store

Open CCI store

In [2]:

Copied!

cci_store = new_data_store("ccizarr")
cci_store = new_data_store("ccizarr")

Open s3 store and list all datasets

In [3]:

Copied!

root = "deep-esdl-public"
root = "deep-esdl-public"

In [4]:

Copied!

s3_store = new_data_store("s3", root=root)
s3_store = new_data_store("s3", root=root)

In [5]:

Copied!

list(s3_store.get_data_ids())
list(s3_store.get_data_ids())

Out[5]:

['LC-1x2160x2160-1.0.0.levels',
 'SMOS-L2C-OS-20230101-20231231-1W-res0-1x1000x1000.levels',
 'SMOS-L2C-OS-20230101-20231231-1W-res0-53x120x120.zarr',
 'SMOS-L2C-OS-20230101-20231231-1W-res0.zarr',
 'SMOS-L2C-SM-20230101-20231231-1W-res0-1x1000x1000.levels',
 'SMOS-L2C-SM-20230101-20231231-1W-res0-53x120x120.zarr',
 'SMOS-L2C-SM-20230101-20231231-1W-res0.zarr',
 'SMOS-freezethaw-1x720x720-1.0.1.zarr',
 'SMOS-freezethaw-4267x10x10-1.0.1.zarr',
 'SeasFireCube_v3.zarr',
 'black-sea-1x1024x1024.levels',
 'black-sea-256x128x128.zarr',
 'esa-cci-permafrost-1x1151x1641-0.0.2.levels',
 'esdc-8d-0.25deg-1x720x1440-3.0.1.zarr',
 'esdc-8d-0.25deg-256x128x128-3.0.1.zarr',
 'hydrology-1D-0.009deg-100x60x60-3.0.2.zarr',
 'hydrology-1D-0.009deg-1418x70x76-2.0.0.zarr',
 'hydrology-1D-0.009deg-1x1102x2415-2.0.0.levels',
 'hydrology-1D-0.009deg-1x1102x966-3.0.2.levels',
 'ocean-1M-9km-1x1080x1080-1.4.0.levels',
 'ocean-1M-9km-64x256x256-1.4.0.zarr',
 'polar-100m-1x2048x2048-1.0.1.zarr']

Open Land Cover (LC) dataset from s3 bucket, which is saved as a multilevel dataset:

In [6]:

Copied!

ml_LC = s3_store.open_data("LC-1x2160x2160-1.0.0.levels", decode_cf=True)
ml_LC = s3_store.open_data("LC-1x2160x2160-1.0.0.levels", decode_cf=True)

Lets open level 0, which is the base level and therefore has the highest resolution:

In [7]:

Copied!

LC = ml_LC.get_dataset(0)
LC = ml_LC.get_dataset(0)

Now let's search for soil datasets provided via the xcube cci store:

In [8]:

Copied!

cci_store.list_data_ids()
cci_store.list_data_ids()

Out[8]:

['ESACCI-BIOMASS-L4-AGB-MERGED-100m-2010-2018-fv2.0.zarr',
 'ESACCI-BIOMASS-L4-AGB-MERGED-100m-2010-2020-fv4.0.zarr',
 'ESACCI-GHG-L2-CH4-SCIAMACHY-WFMD-2002-2011-fv1.zarr',
 'ESACCI-GHG-L2-CO2-OCO-2-FOCAL-2014-2021-v10.zarr',
 'ESACCI-GHG-L2-CO2-SCIAMACHY-WFMD-2002-2012-fv1.zarr',
 'ESACCI-ICESHEETS_Antarctica_GMB-2002-2016-v1.1.zarr',
 'ESACCI-ICESHEETS_Greenland_GMB-2003-2016-v1.1.zarr',
 'ESACCI-L3C_CLOUD-CLD_PRODUCTS-AVHRR_NOAA-1982-2016-fv3.0.zarr',
 'ESACCI-L3C_SNOW-SWE-1979-2018-fv1.0.zarr',
 'ESACCI-L3C_SNOW-SWE-1979-2020-fv2.0.zarr',
 'ESACCI-L4_GHRSST-SST-GMPE-GLOB_CDR2.0-1981-2016-v02.0-fv01.0.zarr',
 'ESACCI-LC-L4-LCCS-Map-300m-P1Y-1992-2015-v2.0.7b.zarr',
 'ESACCI-LST-L3C-LST-MODISA-0.01deg_1DAILY_DAY-2002-2018-fv3.00.zarr',
 'ESACCI-LST-L3C-LST-MODISA-0.01deg_1DAILY_NIGHT-2002-2018-fv3.00.zarr',
 'ESACCI-LST-L3S-LST-IRCDR_-0.01deg_1DAILY_DAY-1995-2020-fv3.00.zarr',
 'ESACCI-LST-L3S-LST-IRCDR_-0.01deg_1DAILY_NIGHT-1995-2020-fv3.00.zarr',
 'ESACCI-LST-L3S-LST-IRCDR_-0.01deg_1MONTHLY_DAY-1995-2020-fv3.00.zarr',
 'ESACCI-LST-L3S-LST-IRCDR_-0.01deg_1MONTHLY_NIGHT-1995-2020-fv3.00.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-1M_MONTHLY_4km_GEO_PML_OCx_QAA-1997-2020-fv5.0.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-1M_MONTHLY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-1D_DAILY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-1Y_YEARLY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-5D_DAILY_4km_GEO_PML_OCx_QAA-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-IOP-MERGED-8D_DAILY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-OC_PRODUCTS-MERGED-1D_DAILY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-OC_PRODUCTS-MERGED-1M_MONTHLY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-OC_PRODUCTS-MERGED-1Y_YEARLY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-OC_PRODUCTS-MERGED-5D_DAILY_4km_GEO_PML_OCx_QAA-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-OC_PRODUCTS-MERGED-8D_DAILY_4km_GEO_PML_OCx_QAA-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-RRS-MERGED-1M_MONTHLY_4km_GEO_PML_RRS-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-RRS-MERGED-1Y_YEARLY_4km_GEO_PML_RRS-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-RRS-MERGED-1D_DAILY_4km_GEO_PML_RRS-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-RRS-MERGED-5D_DAILY_4km_GEO_PML_RRS-1997-2022-fv6.0.zarr',
 'ESACCI-OC-L3S-RRS-MERGED-8D_DAILY_4km_GEO_PML_RRS-1997-2022-fv6.0.zarr',
 'ESACCI-PERMAFROST-L4-ALT-MODISLST-AREA4_PP-1997-2018-fv02.0.zarr',
 'ESACCI-SEAICE-L3C-SITHICK-RA2_ENVISAT-NH25KMEASE2-2002-2012-fv2.0.zarr',
 'ESACCI-SEAICE-L3C-SITHICK-SIRAL_CRYOSAT2-NH25KMEASE2-2010-2017-fv2.0.zarr',
 'ESACCI-SEAICE-L4-SICONC-AMSR_50.0kmEASE2-NH-2002-2017-fv2.1.zarr',
 'ESACCI-SEALEVEL-IND-MSLTR-MERGED-1993-2016-fv02.zarr',
 'ESACCI-SEALEVEL-L4-MSLA-MERGED-1993-2015-fv02.zarr',
 'ESACCI-SOILMOISTURE-L3S-SSMV-COMBINED-1978-2020-fv05.3.zarr',
 'ESACCI-SOILMOISTURE-L3S-SSMV-COMBINED-1978-2021-fv07.1.zarr']

Open Soil Moisture (SM) dataset from xcube cci store:

In [9]:

Copied!

SM = cci_store.open_data("ESACCI-SOILMOISTURE-L3S-SSMV-COMBINED-1978-2021-fv07.1.zarr")
SM = cci_store.open_data("ESACCI-SOILMOISTURE-L3S-SSMV-COMBINED-1978-2021-fv07.1.zarr")

In [10]:

Copied!

SM
SM

Out[10]:

<xarray.Dataset> Size: 654GB
Dimensions:         (time: 15767, lat: 720, lon: 1440)
Coordinates:
  * lat             (lat) float64 6kB 89.88 89.62 89.38 ... -89.38 -89.62 -89.88
  * lon             (lon) float64 12kB -179.9 -179.6 -179.4 ... 179.6 179.9
  * time            (time) datetime64[ns] 126kB 1978-11-01 ... 2021-12-31
Data variables:
    dnflag          (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    flag            (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    freqbandID      (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    mode            (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    sensor          (time, lat, lon) float64 131GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    sm              (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    sm_uncertainty  (time, lat, lon) float32 65GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
    t0              (time, lat, lon) float64 131GB dask.array<chunksize=(16, 720, 720), meta=np.ndarray>
Attributes: (12/44)
    Conventions:                  CF-1.9
    cdm_data_type:                Grid
    comment:                      This dataset was produced with funding of t...
    contact:                      cci_sm_contact@eodc.eu
    creator_email:                cci_sm_developer@eodc.eu
    creator_name:                 Department of Geodesy and Geoinformation, V...
    ...                           ...
    time_coverage_end_product:    20211231T235959Z
    time_coverage_resolution:     P1D
    time_coverage_start:          1978-11-01 00:00:00
    time_coverage_start_product:  19781101T000000Z
    title:                        ESA CCI Surface Soil Moisture COMBINED acti...
    tracking_id:                  ad35798e-58e0-488f-b5b9-593874a47700

Subsetting time and space for the sake of an efficient example

Temporal subset, entire 2020:

In [11]:

Copied!

start_date = datetime.datetime(2020, 1, 1)
stop_date = datetime.datetime(2020, 12, 31)
start_date = datetime.datetime(2020, 1, 1)
stop_date = datetime.datetime(2020, 12, 31)

Spatial subset:

In [12]:

Copied!





min_lat = 70.0
max_lat = 60.0
min_lon = 15
max_lon = 25.0
min_lat = 70.0
max_lat = 60.0
min_lon = 15
max_lon = 25.0

Subsetting the Land Cover dataset, results in only one time slice per year for LC:

In [13]:

Copied!





LC = LC.sel(
    lat=slice(min_lat, max_lat),
    lon=slice(min_lon, max_lon),
    time=slice(start_date, stop_date),
)
LC = LC.sel(
    lat=slice(min_lat, max_lat),
    lon=slice(min_lon, max_lon),
    time=slice(start_date, stop_date),
)

In [14]:

Copied!

LC
LC

Out[14]:

<xarray.Dataset> Size: 156MB
Dimensions:              (time: 1, lat: 3600, lon: 3600, bounds: 2)
Coordinates:
  * lat                  (lat) float64 29kB 70.0 70.0 69.99 ... 60.01 60.0 60.0
  * lon                  (lon) float64 29kB 15.0 15.0 15.01 ... 24.99 25.0 25.0
  * time                 (time) datetime64[ns] 8B 2020-01-01
Dimensions without coordinates: bounds
Data variables:
    change_count         (time, lat, lon) uint8 13MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>
    crs                  (time) int32 4B dask.array<chunksize=(1,), meta=np.ndarray>
    current_pixel_state  (time, lat, lon) float32 52MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>
    lat_bounds           (time, lat, bounds) float64 58kB dask.array<chunksize=(1, 1440, 2), meta=np.ndarray>
    lccs_class           (time, lat, lon) uint8 13MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>
    lon_bounds           (time, lon, bounds) float64 58kB dask.array<chunksize=(1, 1080, 2), meta=np.ndarray>
    observation_count    (time, lat, lon) uint16 26MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>
    processed_flag       (time, lat, lon) float32 52MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>
    time_bounds          (time, bounds) datetime64[ns] 16B dask.array<chunksize=(1, 2), meta=np.ndarray>
Attributes: (12/38)
    Conventions:                CF-1.6
    TileSize:                   2025:2025
    cdm_data_type:              grid
    comment:                    
    contact:                    https://www.ecmwf.int/en/about/contact-us/get...
    creation_date:              20181130T095431Z
    ...                         ...
    time_coverage_end:          20101231
    time_coverage_resolution:   P1Y
    time_coverage_start:        20100101
    title:                      Land Cover Map of ESA CCI brokered by CDS
    tracking_id:                96ac9aca-1ca7-45c6-b4a5-ab448c692646
    type:                       ESACCI-LC-L4-LCCS-Map-300m-P1Y

Subsetting the Soil Moisture dataset, choosing only one month here:

In [15]:

Copied!

start_date = datetime.datetime(2020, 8, 1)
stop_date = datetime.datetime(2020, 8, 31)
start_date = datetime.datetime(2020, 8, 1)
stop_date = datetime.datetime(2020, 8, 31)

In [16]:

Copied!





SM = SM.sel(
    lat=slice(min_lat, max_lat),
    lon=slice(min_lon, max_lon),
    time=slice(start_date, stop_date),
)
SM = SM.sel(
    lat=slice(min_lat, max_lat),
    lon=slice(min_lon, max_lon),
    time=slice(start_date, stop_date),
)

In [17]:

Copied!

SM
SM

Out[17]:

<xarray.Dataset> Size: 2MB
Dimensions:         (time: 31, lat: 40, lon: 40)
Coordinates:
  * lat             (lat) float64 320B 69.88 69.62 69.38 ... 60.62 60.38 60.12
  * lon             (lon) float64 320B 15.12 15.38 15.62 ... 24.38 24.62 24.88
  * time            (time) datetime64[ns] 248B 2020-08-01 ... 2020-08-31
Data variables:
    dnflag          (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    flag            (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    freqbandID      (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    mode            (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    sensor          (time, lat, lon) float64 397kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    sm              (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    sm_uncertainty  (time, lat, lon) float32 198kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
    t0              (time, lat, lon) float64 397kB dask.array<chunksize=(15, 40, 40), meta=np.ndarray>
Attributes: (12/44)
    Conventions:                  CF-1.9
    cdm_data_type:                Grid
    comment:                      This dataset was produced with funding of t...
    contact:                      cci_sm_contact@eodc.eu
    creator_email:                cci_sm_developer@eodc.eu
    creator_name:                 Department of Geodesy and Geoinformation, V...
    ...                           ...
    time_coverage_end_product:    20211231T235959Z
    time_coverage_resolution:     P1D
    time_coverage_start:          1978-11-01 00:00:00
    time_coverage_start_product:  19781101T000000Z
    title:                        ESA CCI Surface Soil Moisture COMBINED acti...
    tracking_id:                  ad35798e-58e0-488f-b5b9-593874a47700

Resample to the same Grid

Here we use xcube's GridMapping method to extract the specification of both grids. Soil Moisture is the grid to be transformed.

In [18]:

Copied!

source_gm = GridMapping.from_dataset(SM)
source_gm = GridMapping.from_dataset(SM)

In [19]:

Copied!

source_gm
source_gm

Out[19]:

is_regular: True
is_j_axis_up: False
is_lon_360: False
crs: EPSG:4326
xy_res: (0.25, 0.25)
xy_bbox: (15, 60, 25, 70)
ij_bbox: (0, 0, 40, 40)
xy_dim_names: ('lon', 'lat')
xy_var_names: ('lon', 'lat')
size: (40, 40)
tile_size: (40, 40)

The target grid mapping is the one from the Land Cover dataset:

In [20]:

Copied!

target_gm = GridMapping.from_dataset(LC)
target_gm = GridMapping.from_dataset(LC)

In [21]:

Copied!

target_gm
target_gm

Out[21]:

is_regular: True
is_j_axis_up: False
is_lon_360: False
crs: EPSG:4326
xy_res: (0.002777775, 0.002777775)
xy_bbox: (15, 60, 25, 70)
ij_bbox: (0, 0, 3600, 3600)
xy_dim_names: ('lon', 'lat')
xy_var_names: ('lon', 'lat')
size: (3600, 3600)
tile_size: (2160, 2160)

Now we resample Soil Moisture to the grid provided by Land Cover:

In [22]:

Copied!

resampled_SM = resample_in_space(SM, source_gm=source_gm, target_gm=target_gm)
resampled_SM = resample_in_space(SM, source_gm=source_gm, target_gm=target_gm)

Lets compare the different grid mappings, to see whether the resampled Soil Moisture datasets has the desired grid mapping now:

The grid mapping of Land Cover which was the target grid mapping:

In [23]:

Copied!

target_gm
target_gm

Out[23]:

is_regular: True
is_j_axis_up: False
is_lon_360: False
crs: EPSG:4326
xy_res: (0.002777775, 0.002777775)
xy_bbox: (15, 60, 25, 70)
ij_bbox: (0, 0, 3600, 3600)
xy_dim_names: ('lon', 'lat')
xy_var_names: ('lon', 'lat')
size: (3600, 3600)
tile_size: (2160, 2160)

The grid mapping of our resampled Soil Moisture, which is as expeced the same as the target grid mapping:

In [24]:

Copied!

GridMapping.from_dataset(resampled_SM)
GridMapping.from_dataset(resampled_SM)

Out[24]:

is_regular: True
is_j_axis_up: False
is_lon_360: False
crs: EPSG:4326
xy_res: (0.002777775, 0.002777775)
xy_bbox: (15, 60, 25, 70)
ij_bbox: (0, 0, 3600, 3600)
xy_dim_names: ('lon', 'lat')
xy_var_names: ('lon', 'lat')
size: (3600, 3600)
tile_size: (2160, 2160)

Create a mask from LC classes

Now that we have our two datasets with the same grid mapping we can use the land cover classes to mask the Soil Moisture datasets.

Converting as they need to be an numpy array to be used with MaskSet:

In [25]:

Copied!

LC.lccs_class.attrs["flag_values"] = np.array(LC.lccs_class.attrs["flag_values"])
LC.lccs_class.attrs["flag_values"] = np.array(LC.lccs_class.attrs["flag_values"])

Creating a Land Cover mask from the classes using xcube's MaskSet:

In [26]:

Copied!

LC_mask = MaskSet(LC.lccs_class)
LC_mask = MaskSet(LC.lccs_class)

In [27]:

Copied!

LC_mask
LC_mask

Out[27]:

Flag name	Mask	Value
no_data	None	0
cropland_rainfed	None	10
cropland_rainfed_herbaceous_cover	None	11
cropland_rainfed_tree_or_shrub_cover	None	12
cropland_irrigated	None	20
mosaic_cropland	None	30
mosaic_natural_vegetation	None	40
tree_broadleaved_evergreen_closed_to_open	None	50
tree_broadleaved_deciduous_closed_to_open	None	60
tree_broadleaved_deciduous_closed	None	61
tree_broadleaved_deciduous_open	None	62
tree_needleleaved_evergreen_closed_to_open	None	70
tree_needleleaved_evergreen_closed	None	71
tree_needleleaved_evergreen_open	None	72
tree_needleleaved_deciduous_closed_to_open	None	80
tree_needleleaved_deciduous_closed	None	81
tree_needleleaved_deciduous_open	None	82
tree_mixed	None	90
mosaic_tree_and_shrub	None	100
mosaic_herbaceous	None	110
shrubland	None	120
shrubland_evergreen	None	121
shrubland_deciduous	None	122
grassland	None	130
lichens_and_mosses	None	140
sparse_vegetation	None	150
sparse_tree	None	151
sparse_shrub	None	152
sparse_herbaceous	None	153
tree_cover_flooded_fresh_or_brakish_water	None	160
tree_cover_flooded_saline_water	None	170
shrub_or_herbaceous_cover_flooded	None	180
urban	None	190
bare_areas	None	200
bare_areas_consolidated	None	201
bare_areas_unconsolidated	None	202
water	None	210
snow_and_ice	None	220

For the masking to work, both data arrays must have identical coordinates, we thus enforce this condition by assigning coordinates and merging into one dataset

For the masking example, one specific class, rainfed cropland, is choosen here:

In [28]:

Copied!

LC_cropmask = LC_mask.cropland_rainfed.to_dataset()
LC_cropmask = LC_mask.cropland_rainfed.to_dataset()

Although the coordinates of the LC dataset and the resampled_SM dataset are almost identical, tiny numerical differences would prevent the masking from working. Thus, the coordinates from the resampled_SM dataset are assigned as the coordinates to the LC_cropmask

In [29]:

Copied!

LC_cropmask = LC_cropmask.assign_coords(lat=resampled_SM.lat, lon=resampled_SM.lon)
LC_cropmask = LC_cropmask.assign_coords(lat=resampled_SM.lat, lon=resampled_SM.lon)

The resulting dataset is again a data cube with only one time slices

In [30]:

Copied!

LC_cropmask
LC_cropmask

Out[30]:

<xarray.Dataset> Size: 13MB
Dimensions:           (time: 1, lat: 3600, lon: 3600)
Coordinates:
  * time              (time) datetime64[ns] 8B 2020-01-01
  * lat               (lat) float64 29kB 70.0 70.0 69.99 ... 60.01 60.0 60.0
  * lon               (lon) float64 29kB 15.0 15.0 15.01 ... 24.99 25.0 25.0
Data variables:
    cropland_rainfed  (time, lat, lon) uint8 13MB dask.array<chunksize=(1, 1440, 1080), meta=np.ndarray>

Because both datasets share the same spatial grid, we can now insert the cropland_rainfed mask as an additional data variable into the resampled_SM datacube. By selecting the only time slice, we remove the time information here. The new data variable crop_mask in the resulting dataset thus only depending on lat and lon.

In [31]:

Copied!

resampled_SM["crop_mask"] = LC_cropmask.cropland_rainfed.isel(time=0)
resampled_SM["crop_mask"] = LC_cropmask.cropland_rainfed.isel(time=0)

In [32]:

Copied!

resampled_SM
resampled_SM

Out[32]:

<xarray.Dataset> Size: 16GB
Dimensions:         (time: 31, lat: 3600, lon: 3600)
Coordinates:
  * time            (time) datetime64[ns] 248B 2020-08-01 ... 2020-08-31
  * lon             (lon) float64 29kB 15.0 15.0 15.01 15.01 ... 24.99 25.0 25.0
  * lat             (lat) float64 29kB 70.0 70.0 69.99 69.99 ... 60.01 60.0 60.0
Data variables:
    dnflag          (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    flag            (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    freqbandID      (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    mode            (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    sensor          (time, lat, lon) float64 3GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    sm              (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    sm_uncertainty  (time, lat, lon) float32 2GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    t0              (time, lat, lon) float64 3GB dask.array<chunksize=(1, 2160, 2160), meta=np.ndarray>
    crop_mask       (lat, lon) uint8 13MB dask.array<chunksize=(1440, 1080), meta=np.ndarray>
Attributes: (12/44)
    Conventions:                  CF-1.9
    cdm_data_type:                Grid
    comment:                      This dataset was produced with funding of t...
    contact:                      cci_sm_contact@eodc.eu
    creator_email:                cci_sm_developer@eodc.eu
    creator_name:                 Department of Geodesy and Geoinformation, V...
    ...                           ...
    time_coverage_end_product:    20211231T235959Z
    time_coverage_resolution:     P1D
    time_coverage_start:          1978-11-01 00:00:00
    time_coverage_start_product:  19781101T000000Z
    title:                        ESA CCI Surface Soil Moisture COMBINED acti...
    tracking_id:                  ad35798e-58e0-488f-b5b9-593874a47700

Plot the mask of rainfed croplands (1 True, 0 False).

In [33]:

Copied!

resampled_SM.crop_mask.plot()
resampled_SM.crop_mask.plot()

Out[33]:

<matplotlib.collections.QuadMesh at 0x7fbca8296850>

No description has been provided for this image

Mask one data array with a condition from another, inserted into one dataset

Select Soil Moisture (sm) only for rainfed croplands. The result is again a datacube of the same dimensions as the inputs, because we removed the time information from the cropmask before.

In [34]:

Copied!

SM_crop = resampled_SM.sm.where(
    resampled_SM.crop_mask == 1
)

SM_crop
SM_crop = resampled_SM.sm.where(
    resampled_SM.crop_mask == 1
)

SM_crop

Out[34]:

<xarray.DataArray 'sm' (time: 31, lat: 3600, lon: 3600)> Size: 2GB
dask.array<where, shape=(31, 3600, 3600), dtype=float32, chunksize=(1, 1440, 1080), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 248B 2020-08-01 2020-08-02 ... 2020-08-31
  * lon      (lon) float64 29kB 15.0 15.0 15.01 15.01 ... 24.99 24.99 25.0 25.0
  * lat      (lat) float64 29kB 70.0 70.0 69.99 69.99 ... 60.01 60.01 60.0 60.0
Attributes:
    _CoordinateAxes:      time lat lon
    ancillary_variables:  sm_uncertainty flag t0
    dtype:                float32
    long_name:            Volumetric Soil Moisture
    standard_name:        soil_moisture_content
    units:                m3 m-3
    valid_range:          [0.0, 1.0]

Visualise result: sm only on rainfed cropland

In [35]:

Copied!

SM_crop.isel(time=0).plot()
SM_crop.isel(time=0).plot()

Out[35]:

<matplotlib.collections.QuadMesh at 0x7fbca830dc90>

The histogram of the sm variable on rainfed croplands:

In [36]:

Copied!

SM_crop.isel(time=0).plot.hist(bins=100)
SM_crop.isel(time=0).plot.hist(bins=100)

Out[36]:

(array([1.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 0.0000e+00, 2.0000e+00, 0.0000e+00,
        1.0000e+00, 1.0000e+00, 0.0000e+00, 0.0000e+00, 6.0000e+00,
        4.0000e+00, 0.0000e+00, 1.0000e+01, 1.7000e+01, 2.5000e+01,
        3.0000e+01, 2.1000e+01, 8.2000e+01, 2.0000e+02, 2.7700e+02,
        4.1500e+02, 5.9700e+02, 5.1100e+02, 9.4800e+02, 1.2230e+03,
        1.2290e+03, 1.3000e+03, 1.2380e+03, 1.9540e+03, 4.3730e+03,
        3.7250e+03, 5.9280e+03, 6.6450e+03, 7.2750e+03, 8.0110e+03,
        8.6910e+03, 9.9490e+03, 1.0385e+04, 1.0353e+04, 1.0969e+04,
        1.0769e+04, 8.9880e+03, 8.6950e+03, 9.2770e+03, 8.6980e+03,
        7.1390e+03, 4.7680e+03, 2.8230e+03, 1.0370e+03, 9.5100e+02,
        9.8600e+02, 9.1500e+02, 7.7800e+02, 6.8000e+02, 7.4100e+02,
        7.4200e+02, 7.6600e+02, 7.8800e+02, 8.9100e+02, 9.5500e+02,
        9.9200e+02, 1.1020e+03, 1.2780e+03, 1.1710e+03, 2.1290e+03,
        3.6200e+03, 2.6400e+03, 3.0710e+03, 3.4170e+03, 8.3040e+03,
        2.4350e+03, 2.2000e+03, 2.1750e+03, 2.2330e+03, 2.1700e+03,
        2.0010e+03, 1.7240e+03, 1.2250e+03, 1.3030e+03, 1.2080e+03,
        1.1760e+03, 7.2000e+02, 5.9800e+02, 1.9600e+02, 2.9000e+01]),
 array([0.17464125, 0.17603008, 0.1774189 , 0.17880774, 0.18019657,
        0.18158539, 0.18297422, 0.18436304, 0.18575187, 0.1871407 ,
        0.18852952, 0.18991835, 0.19130719, 0.19269601, 0.19408484,
        0.19547367, 0.19686249, 0.19825132, 0.19964015, 0.20102897,
        0.20241781, 0.20380664, 0.20519546, 0.20658429, 0.20797312,
        0.20936194, 0.21075077, 0.21213959, 0.21352842, 0.21491724,
        0.21630608, 0.21769491, 0.21908373, 0.22047256, 0.22186139,
        0.22325021, 0.22463904, 0.22602788, 0.22741669, 0.22880553,
        0.23019436, 0.23158318, 0.23297201, 0.23436084, 0.23574966,
        0.23713849, 0.23852733, 0.23991615, 0.24130498, 0.24269381,
        0.24408263, 0.24547145, 0.2468603 , 0.24824911, 0.24963793,
        0.25102675, 0.2524156 , 0.25380442, 0.25519323, 0.25658208,
        0.2579709 , 0.25935972, 0.26074857, 0.26213738, 0.2635262 ,
        0.26491505, 0.26630387, 0.26769269, 0.26908153, 0.27047035,
        0.27185917, 0.27324802, 0.27463683, 0.27602565, 0.2774145 ,
        0.27880332, 0.28019214, 0.28158098, 0.2829698 , 0.28435862,
        0.28574747, 0.28713629, 0.2885251 , 0.28991395, 0.29130277,
        0.29269159, 0.29408044, 0.29546925, 0.29685807, 0.29824692,
        0.29963574, 0.30102456, 0.3024134 , 0.30380222, 0.30519104,
        0.30657989, 0.30796871, 0.30935752, 0.31074637, 0.31213516,
        0.31352401]),
 <BarContainer object of 100 artists>)

For comparison, we also visualise the full soil moisture product, without masking condition.

In [37]:

Copied!

SM.sm.isel(time=0).plot()
SM.sm.isel(time=0).plot()

Out[37]:

<matplotlib.collections.QuadMesh at 0x7fbca810dd50>

In [ ]: