# Atlas notebooks

***

> This notebook reproduces and extends parts of the figures and products of the AR6-WGI Atlas. It is part of a notebook collection available at https://github.com/IPCC-WG1/Atlas for reproducibility and reusability purposes. This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0).
>
> ![Creative Commons License >](https://i.creativecommons.org/l/by/4.0/88x31.png)

Using the 4th version of the IPCC WGI reference regions in Python#

Companion to the paper “An update of IPCC climate reference regions for subcontinental analysis of climate model data: Definition and aggregated datasets” (July 2020).

Mathias Hauser - Institute for Atmospheric and Climate Science, ETH Zurich, Zurich, Switzerland

A new set of reference regions is defined in Iturbide et al. (2020) for reporting sub-continental climate information. Officially named IPCC-WGI-reference-regions-v4 they will be referred to as AR6 regions here. These regions update the IPCC AR5 reference regions [AR5]. This tutorial shows how the python package regionmask [Create Masks of Geographical Regions — Regionmask 0.9.0 Documentation] can be used to plot these regions, and to create masks for arbitrary latitude and longitude grids (Version 0.9.0 was used [Hauser et al., 2022]). These masks can then be used to subset individual regions or to calculate regional averages. The package xarray [Xarray Documentation] is used to load netCDF files [Hoyer and Hamman, 2017].

This tutorial requires regionmask version 0.6.1 or later and xarray version 0.15.1 or later. Check the documentation of regionmask for details.

Imports#

Import regionmask and check the version:

import regionmask
regionmask.__version__
Hide code cell output 
'0.7.0'

Import xarray and check the version:

import xarray as xr
xr.set_options(display_style="text")
xr.__version__
Hide code cell output 
'0.19.0'

Import additional packages

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

np.set_printoptions(edgeitems=2)

Loading the reference regions#

The regions are available at regionmask.defined_regions.ar6. The whole set of 58 regions is available under (ar6.all). In addition the land (ar6.land) and ocean (ar6.ocean) regions are given separately. The numbering is kept consistent between the categories. Note that some regions are in the land and in the ocean categories (e.g. the Mediterranean).

regionmask.defined_regions.ar6

AR6 reference regions - Iturbide et al., 2020

Attributes
----------
all : Regions
    All regions (land + ocean).
land : Regions
    Land regions only
ocean : Regions
    Ocean regions only

We will illustrate the use of regionmask with regionmask.defined_regions.ar6.all

ar6_all = regionmask.defined_regions.ar6.all
ar6_all

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/regions.py:410: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for p in poly:

<regionmask.Regions>
Name:     AR6 reference regions
Source:   Iturbide et al., 2020 (Earth Syst. Sci. Data)

Regions:
 0 GIC      Greenland/Iceland
 1 NWN      N.W.North-America
 2 NEN      N.E.North-America
 3 WNA        W.North-America
 4 CNA        C.North-America
..  ..                    ...
53 ARS            Arabian-Sea
54 BOB          Bay-of-Bengal
55 EIO Equatorial.Indic-Ocean
56 SIO          S.Indic-Ocean
57 SOO         Southern-Ocean

[58 regions]

Plotting#

ar6_all.plot() creates an cartopy map plot including the outline of all regions. It returns an axes instance.

ax = ar6_all.plot()

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/plot.py:28: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.
  coords = _polygons_coords(polygons)
../_images/31338b9ecc442f2cba2a41512ddb3cf1f87de9ccb1e082b6b791f504bcf10a8d.png

The plot can also be customized:

f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.Robinson()))

text_kws = dict(color="#67000d", fontsize=8, bbox=dict(pad=0.2, color="w"))

ax = ar6_all.plot(
    ax=ax,
    add_ocean=False,
    line_kws=dict(linewidth=1),
    coastlines=False,
    text_kws=text_kws,
)

ax.coastlines(color="0.5", lw=0.5);

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
../_images/81f3f3aaddf541155f31bcc91c34bb8990603549894f9b59251dbbb96cc82236.png

Selecting Regions#

Selecting regions is done by indexing (note the double brackets):

new_zealand = ar6_all[[43]]
new_zealand

<regionmask.Regions>
Name:     AR6 reference regions
Source:   Iturbide et al., 2020 (Earth Syst. Sci. Data)

Regions:
43 NZ New-Zealand

[1 regions]

We plot the selected region and color the ocean blue with add_ocean=True:

projection = ccrs.PlateCarree(central_longitude=180)

ax = new_zealand.plot(proj=projection, label="abbrev", add_ocean=True)

ax.set_extent([120, 185, -20, -60], ccrs.PlateCarree())

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/plot.py:28: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.
  coords = _polygons_coords(polygons)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
../_images/109d5556f6c58ffa7f28faa2c90859f5901e033f4bef88e6a4d24ecc4bee6dd3.png

Regions can also be selected by their abbreviation or name, let’s select several regions at once:

australasia = ar6_all[["NZ", "SEA", "NAU", "C.Australia", "SAU"]]

ax = australasia.plot(proj=projection, label="abbrev", add_ocean=True)
../_images/3921e21c331427a5d1244e4919bc827898f5a052c2eaf7cfe208c012bb7e8979.png

Example Dataset#

In this example we will consider a NetCDF file with historical temperature simulations and show several analysis and visualizations at a global scale or filtering by regions. Therefore, we load an example dataset from the CMIP6 archive using xarray. The dataset corresponds to the following ESGF request:

  • MIP Era: CMIP6

  • Source ID: CESM2

  • Experiment ID: historical

  • Variant Label: r1i1p1f1

  • Table ID: Amon

  • Frequency: mon

  • Variable: tas

Note: if you do not have access to a tas dataset you can replace the code cell below with the following:

ds = xr.tutorial.load_dataset("air_temperature")
tas = ds.air
# calculate monthly mean
tas = tas.resample(time="m").mean()

fN = ("auxiliary-material/CMIP6Amon_tas_CanESM5_r1i1p1f1_historical_gn_185001-201412.nc") 

ds = xr.open_dataset(fN)

tas = ds.tas

# calculate annual mean
tas = tas.groupby("time.year").mean("time")
tas = tas.rename(year="time")

# convert to celsius
tas = tas - 273.15

tas
Hide code cell output 
<xarray.DataArray 'tas' (time: 165, lat: 64, lon: 128)>
array([[[-45.430008, -45.761246, ..., -44.75174 , -45.091705],
        [-45.896744, -46.590225, ..., -44.218903, -45.101715],
        ...,
        [-17.36702 , -17.171188, ..., -17.895187, -17.624313],
        [-18.64714 , -18.557785, ..., -18.7753  , -18.694626]],

       [[-44.4169  , -44.72325 , ..., -43.792145, -44.099274],
        [-44.805893, -45.43962 , ..., -43.2538  , -44.080124],
        ...,
        [-16.245056, -16.032959, ..., -16.703278, -16.45105 ],
        [-17.691727, -17.641739, ..., -17.772202, -17.741959]],

       ...,

       [[-42.882904, -43.175674, ..., -42.251236, -42.573135],
        [-43.605988, -44.256134, ..., -41.999786, -42.84581 ],
        ...,
        [-14.617554, -14.334442, ..., -15.218384, -14.910461],
        [-15.229614, -15.190063, ..., -15.348999, -15.286438]],

       [[-42.611115, -42.900314, ..., -42.01175 , -42.32065 ],
        [-42.784317, -43.42357 , ..., -41.23117 , -42.06967 ],
        ...,
        [-13.792969, -13.55722 , ..., -14.30069 , -14.045441],
        [-14.829285, -14.778809, ..., -14.893127, -14.84668 ]]],
      dtype=float32)
Coordinates:
  * lat      (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
  * lon      (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
    height   float64 ...
  * time     (time) int64 1850 1851 1852 1853 1854 ... 2010 2011 2012 2013 2014

The example data is a temperature field. Let’s plot the first time step:

proj = ccrs.Robinson()

f, ax = plt.subplots(subplot_kw=dict(projection=proj))

h = tas.isel(time=0).plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), robust=True, center=0
)

ax.coastlines();

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
../_images/6c4cd9348a45ba6fb65ee01a34c4fed1529d840d32676cbafd30eba3094d3bb9.png

Creating a 2D integer mask#

We illustrate the creation of masks using the land regions:

ar6_land = regionmask.defined_regions.ar6.land

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/regions.py:410: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for p in poly:

Using ar6_land.mask(lon, lat) we can create an array where each grid cell is numbered according to the region it belongs to. Grid cells that are not part of any region are NaN. We can directly pass an xarray object containing lon and lat coordinates to the mask function.

mask_2D = ar6_land.mask(tas)

The resulting mask indicates which region each grid point belongs to and looks as follows:

proj = ccrs.Robinson()
f, ax = plt.subplots(subplot_kw=dict(projection=proj))

h = mask_2D.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)

ax.coastlines()

ar6_land.plot_regions(line_kws=dict(lw=0.5), add_label=False);

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/plot.py:28: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.
  coords = _polygons_coords(polygons)
../_images/39ed583ce6c5c0334eb30ef6b8bfa893773f0529e7e457a80a7e8b735821e0b9.png

Creating a 3D boolean mask#

2D masks are good for plotting. However, to calculate weighted regional averages 3D boolean masks are more convenient. 3D masks are created as follows:

mask_3D = ar6_land.mask_3D(tas)
mask_3D
Hide code cell output 
<xarray.DataArray 'region' (region: 46, lat: 64, lon: 128)>
array([[[False, False, ..., False, False],
        [False, False, ..., False, False],
        ...,
        [False, False, ..., False, False],
        [False, False, ..., False, False]],

       [[False, False, ..., False, False],
        [False, False, ..., False, False],
        ...,
        [False, False, ..., False, False],
        [False, False, ..., False, False]],

       ...,

       [[ True,  True, ...,  True,  True],
        [ True,  True, ...,  True,  True],
        ...,
        [False, False, ..., False, False],
        [False, False, ..., False, False]],

       [[False, False, ..., False, False],
        [False, False, ..., False, False],
        ...,
        [False, False, ..., False, False],
        [False, False, ..., False, False]]])
Coordinates:
  * lat      (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86
  * lon      (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
  * region   (region) int64 0 1 2 3 4 5 6 7 8 9 ... 37 38 39 40 41 42 43 44 45
    abbrevs  (region) <U4 'GIC' 'NWN' 'NEN' 'WNA' ... 'SAU' 'NZ' 'EAN' 'WAN'
    names    (region) <U22 'Greenland/Iceland' ... 'W.Antarctica'

ar6_land.mask_3D uses the same algorithm to determine if a gridpoint is in a region as for the 2D mask. However, it returns a xarray.Dataset with shape region x lat x lon. Gridpoints that belong to a region are True, all others are False.

The first two regions look as follows:

from matplotlib import colors as mplc

cmap1 = mplc.ListedColormap(["none", "#9ecae1"])

fg = mask_3D.isel(region=slice(2)).plot(
    subplot_kws=dict(projection=ccrs.PlateCarree()),
    col="region",
    col_wrap=2,
    transform=ccrs.PlateCarree(),
    add_colorbar=False,
    aspect=1.5,
    cmap=cmap1,
)

for ax in fg.axes.flatten():
    ax.coastlines(lw=0.5, color="0.5")

fg.fig.subplots_adjust(hspace=0, wspace=0.1);

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
../_images/fa4c37084c2615e6e6d5e7d286c9c655babc535bdce4098c7d9d5c6fd14d9e6e.png

Select a region of a 3D mask#

As mask_3D contains region, abbrevs, and names as (non-dimension) coordinates we can use each of those to select an individual region:

# 1) by the index of the region:
r1 = mask_3D.sel(region=2)

# 2) with the abbreviation
r2 = mask_3D.isel(region=(mask_3D.abbrevs == "CNA"))

# 3) with the long name:
r3 = mask_3D.isel(region=(mask_3D.names == "C.North-America"))

Mask out a region#

Using tas.where a specific region can be ‘masked out’ (i.e. all data points outside of the region become NaN):

tas_CNA = tas.where(r1)

proj = ccrs.Robinson()

ax = plt.subplot(111, projection=proj)

tas_CNA.isel(time=0).plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())

ax.set_title("Central North America")

ax.coastlines();

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
../_images/516ff201bda6e298931c6c47d3fdcc34a9974019aa7551f8ea22eaf9e66d359a.png

Note this also works with a 2D mask, using tas.where(mask_2D == 2).

Calculate weighted regional averages#

Using the 3-dimensional mask it is possible to calculate weighted averages of all regions in one go, using the weighted method (requires xarray 0.15.1 or later). Here, we use cos(lat) as proxy of the grid area. This works well for the rectangular grid in our example. In general, it is better to the model’s own grid area.

weights = np.cos(np.deg2rad(tas.lat))

tas_regional = tas.weighted(mask_3D * weights).mean(dim=("lat", "lon"))

Let’s break down what happens here. By multiplying mask_3D * weights we get a DataArray where gridpoints within a region get a weight proportional to the gridcell area and all others get a weight of 0. airtemps.weighted(mask_3D * weights) creates an xarray object which can be used for weighted operations. From this we calculate the weighted mean over the lat and lon dimensions. The resulting dataarray has the dimensions region x time:

tas_regional
Hide code cell output 
<xarray.DataArray 'tas' (time: 165, region: 46)>
array([[-14.34435989,  -4.57573997, ..., -32.31694793, -20.89343577],
       [-14.05052335,  -5.25600597, ..., -31.61877546, -21.7574261 ],
       ...,
       [-12.404766  ,  -2.08977824, ..., -29.92128283, -18.47642827],
       [-11.14616891,  -2.25170586, ..., -29.97843364, -19.29482384]])
Coordinates:
    height   float64 2.0
  * time     (time) int64 1850 1851 1852 1853 1854 ... 2010 2011 2012 2013 2014
  * region   (region) int64 0 1 2 3 4 5 6 7 8 9 ... 37 38 39 40 41 42 43 44 45
    abbrevs  (region) <U4 'GIC' 'NWN' 'NEN' 'WNA' ... 'SAU' 'NZ' 'EAN' 'WAN'
    names    (region) <U22 'Greenland/Iceland' ... 'W.Antarctica'

The regionally-averaged time series can be plotted:

tas_regional.isel(region=slice(6)).plot(col="region", col_wrap=3, sharey=False);
../_images/720a2cc116944acd2b7b72cb631bf7c78835f40783f6ffc2a7b87ce0e0e0396a.png

Restrict the mask to land points#

Combining the mask of the regions with a land-sea mask we can create a mask containing only land points. For simplicity we use the natural_earth.land_110 regions to create a land mask but generally it is recommended to use the model’s original land/ sea mask.

With this caveat in mind we can create the land-sea mask and plot it:

land_110 = regionmask.defined_regions.natural_earth.land_110

land_mask = land_110.mask_3D(tas)

# add a plot
ax = plt.axes(projection=ccrs.Robinson())
land_mask.plot(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False);

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/regions.py:410: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for p in poly:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
../_images/a07d4ee9be44a951aaf6a5027c74b34afe419391b875970603ad91125157caeb.png

To create the combined mask we multiply the two:

mask_lsm = mask_3D * land_mask.squeeze(drop=True)

Note the .squeeze(drop=True). This is required to remove the region dimension from land_mask.

The resulting mask_lsm is also a 3D mask which cannot be directly plotted. For this regionmask.plot_3D_mask is given as convenience function to flatten and plot 3D masks:

ar6_land.plot(add_label=False, line_kws=dict(lw=1), proj=ccrs.Robinson())

regionmask.plot_3D_mask(mask_lsm, transform=ccrs.PlateCarree(), add_colorbar=False);

/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/utils.py:15: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  polys += list(p)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/regionmask/core/plot.py:28: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.
  coords = _polygons_coords(polygons)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:228: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(multi_line_string) > 1:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:239: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings = list(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:280: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  for line in multi_line_string:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:347: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  if len(p_mline) > 0:
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
/home/phanaur/mambaforge/envs/tfg/lib/python3.7/site-packages/cartopy/crs.py:385: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
  line_strings.extend(multi_line_string)
../_images/140a3b8ea7433f5b37c083f7ebe9e80d5e5dc4fed3c2eeeeaf8635fd547f6dd7.png

We can now again calculate the regional mean temperature for all land points:

tas_regional_land = tas.weighted(mask_lsm * weights).mean(dim=("lat", "lon"))

and compare it to the temperature of the land and ocean:

f, axes = plt.subplots(2, 3, sharex=True, figsize=(9, 6))
axes = axes.flatten()

for i, ax in enumerate(axes):

    ds = tas_regional.isel(region=i)
    ds.plot(ax=ax, label="Land + Ocean")

    ds = tas_regional_land.isel(region=i)
    ds.plot(ax=ax, label="Land only")

    ax.set_title(ds.names.values)
    if i < 3:
        ax.set_xlabel("")

axes[2].legend()
plt.tight_layout()
../_images/339217aaa64421aa585cf90bf0ec7b82782bd4122c9931fcb7e812ac2eb6a4c7.png

Calculating weighed averages in xarray prior to version 0.15.1#

xarray only offers weighted averages since version 0.15.1. For users of older versions the a weighted_mean function is provided below. This function can also compute the weighed mean over all regions (weighted_mean(tas, mask_3D * weights, ("lon", "lat"))), however, it is less memory efficient and much slower than the one implemented in xarray.

def weighted_mean(da, weights, dim):
    """Reduce da by a weighted mean along some dimension(s).

    Parameters
    ----------
    da : DataArray
        Object over which the weighted reduction operation is applied.    
    weights : DataArray
        An array of weights associated with the values in this Dataset.
    dim : str or sequence of str, optional
        Dimension(s) over which to apply the weighted `mean`.
        
    Returns
    -------
    weighted_mean : DataArray
        New DataArray with weighted mean applied to its data and
        the indicated dimension(s) removed.
    """

    weighted_sum = (da * weights).sum(dim=dim, skipna=True)
    # need to mask weights where data is not valid
    masked_weights = weights.where(da.notnull())
    sum_of_weights = masked_weights.sum(dim=dim, skipna=True)
    valid_weights = sum_of_weights != 0
    sum_of_weights = sum_of_weights.where(valid_weights)

    return weighted_sum / sum_of_weights

Requirements#

Running this notebook requires the regionmask, netCDF4 and jupyter packages (all other packages are dependencies). regionmask itself is a pure Python package, but its dependencies are not. The easiest way to get them installed is to use conda. The package is avilable on the conda-forge channel.

conda install -c conda-forge regionmask xarray jupyter netCDF4

References#

reg

Create masks of geographical regions — regionmask 0.9.0 documentation. https://regionmask.readthedocs.io/en/stable/.

Xar

Xarray documentation. https://docs.xarray.dev/en/latest/index.html.

AR5

ipcc AR5. AR5 Regions. https://www.ipcc-data.org/guidelines/pages/ar5_regions.html.

HSB+22

Mathias Hauser, Aaron Spring, Julius Busecke, Martin van Driel, Ruth Lorenz, and readthedocs-assistant. Regionmask/regionmask: version 0.9.0. Zenodo, March 2022. doi:10.5281/zenodo.6324240.

HH17

Stephan Hoyer and Joe Hamman. Xarray: N-D labeled Arrays and Datasets in Python. Journal of Open Research Software, 5(1):10, April 2017. doi:10.5334/jors.148.