

In [ ]:

    
%matplotlib inline
%config InlineBackend.figure_format = "retina"

from matplotlib import rcParams

rcParams["savefig.dpi"] = 200
rcParams["font.size"] = 8

import warnings

warnings.filterwarnings("ignore")

Create 2D integer masks

In this tutorial we will show how to create 2D integer mask for arbitrary latitude and longitude grids.

.. note:: 2D masks are good for plotting. However, to calculate weighted regional averages 3D boolean masks are more convenient. See the :doc:`tutorial on 3D masks`.

Import regionmask and check the version:



In [ ]:

    
import regionmask

regionmask.__version__

Load xarray and the tutorial data:



In [ ]:

    
import xarray as xr
import numpy as np

Creating a mask

Define a lon/ lat grid with a 1° grid spacing, where the points define the center of the grid.



In [ ]:

    
lon = np.arange(-179.5, 180)
lat = np.arange(-89.5, 90)

We will create a mask with the SREX regions (Seneviratne et al., 2012).



In [ ]:

    
regionmask.defined_regions.srex

The function mask determines which gripoints lie within the polygon making up the each region:



In [ ]:

    
mask = regionmask.defined_regions.srex.mask(lon, lat)
mask

mask is now a xarray.Dataset with shape lat x lon (if you need a numpy array use mask.values). Gridpoints that do not fall in a region are NaN, the gridpoints that fall in a region are encoded with the number of the region (here 1 to 26).

We can now plot the mask:



In [ ]:

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

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

mask.plot(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False);

Working with a mask

masks can be used to select data in a certain region and to calculate regional averages - let's illustrate this with a 'real' dataset:



In [ ]:

    
airtemps = xr.tutorial.load_dataset("air_temperature")

The example data is a temperature field over North America. Let's plot the first time step:



In [ ]:

    
# choose a good projection for regional maps
proj = ccrs.LambertConformal(central_longitude=-100)

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

airtemps.isel(time=1).air.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())

ax.coastlines();

Conviniently we can directly pass an xarray object to the mask function. It gets the longitude and latitude from the DataArray/ Dataset and creates the mask. If the longitude and latitude in the xarray object are not called "lon" and "lat", respectively; you can pass their name via the lon_name and lat_name keyword.



In [ ]:

    
mask = regionmask.defined_regions.srex.mask(airtemps)

.. note:: From version 0.5 ``regionmask`` automatically detects wether the longitude needs to be wrapped around, i.e. if the regions extend from -180° E to 180° W, while the grid goes from 0° to 360° W as in our example:



In [ ]:

    
lon = airtemps.lon.values
print("Grid extent:    {:3.0f}°E to {:3.0f}°E".format(lon.min(), lon.max()))

bounds = regionmask.defined_regions.srex.bounds_global
print("Region extent: {:3.0f}°E to {:3.0f}°E".format(bounds[0], bounds[2]))

Let's plot the mask of the regions:



In [ ]:

    
proj = ccrs.LambertConformal(central_longitude=-100)
ax = plt.subplot(111, projection=proj)

low = mask.min()
high = mask.max()

levels = np.arange(low - 0.5, high + 1)

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

# for colorbar: find abbreviations of all regions that were selected
reg = np.unique(mask.values)
reg = reg[~np.isnan(reg)]
abbrevs = regionmask.defined_regions.srex[reg].abbrevs

cbar = plt.colorbar(h, orientation="horizontal", fraction=0.075, pad=0.05)

cbar.set_ticks(reg)
cbar.set_ticklabels(abbrevs)
cbar.set_label("Region")

ax.coastlines()

# fine tune the extent
ax.set_extent([200, 330, 10, 75], crs=ccrs.PlateCarree())

We want to select the region 'Central North America'. Thus we first need to find out which number this is:



In [ ]:

    
CNA_index = regionmask.defined_regions.srex.map_keys("C. North America")
CNA_index

Mask out a region

xarray provides the handy where function:



In [ ]:

    
airtemps_CNA = airtemps.where(mask == CNA_index)

Check everything went well by repeating the first plot with the selected region:



In [ ]:

    
# choose a good projection for regional maps
proj = ccrs.LambertConformal(central_longitude=-100)

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

regionmask.defined_regions.srex[["CNA"]].plot(ax=ax, add_label=False)

airtemps_CNA.isel(time=1).air.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())

ax.coastlines();

Looks good - with this we can calculate the region average.

Calculate weighted regional average

From version 0.15.1 xarray includes a function to calculate the weighted mean - we use cos(lat) as proxy of the grid cell area

.. note:: It is better to use a model's original grid cell area (e.g. areacella). ``cos(lat)`` works reasonably well for regular lat/ lon grids. For irregular grids (regional models, ocean models, ...) it is not appropriate.



In [ ]:

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

ts_airtemps_CNA = airtemps_CNA.weighted(weights).mean(dim=("lat", "lon")) - 273.15

We plot the resulting time series:



In [ ]:

    
f, ax = plt.subplots()
ts_airtemps_CNA.air.plot.line(ax=ax, label="Central North America")

ax.axhline(0, color="0.1", lw=0.5)

plt.legend();

To get the regional average for each region you would need to loop over them. However, it's easier to use a 3D mask.

Calculate regional averages using `groupby`

.. warning:: Using ``groupby`` offers some convenience and is faster than using ``where`` and a loop. However, it cannot currently be combinded with ``weighted`` (xarray `GH3937 `_). Therefore, I recommend working with a :doc:`3D mask`.



In [ ]:

    
# you can group over all integer values of the mask
airtemps_all = airtemps.groupby(mask).mean()
airtemps_all

Multidimensional coordinates

Regionmask can also handle mutltidimensional longitude/ latitude grids (e.g. from a regional climate model). As xarray provides such an example dataset, we will use it to illustrate it. See also in the xarray documentation.

Load the tutorial data:



In [ ]:

    
rasm = xr.tutorial.load_dataset("rasm")

The example data is a temperature field over the Northern Hemisphere. Let's plot the first time step:



In [ ]:

    
# choose a projection
proj = ccrs.NorthPolarStereo()

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

rasm.isel(time=1).Tair.plot.pcolormesh(
    ax=ax, x="xc", y="yc", transform=ccrs.PlateCarree()
)

# add the abbreviation of the regions
regionmask.defined_regions.srex.plot(
    ax=ax, regions=[1, 2, 11, 12, 18], coastlines=False, label="abbrev"
)

ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree())

ax.coastlines();

Again we pass the xarray object to regionmask. We have to specify "xc" and "yc" as the longitude and latitude coordinates of the array:



In [ ]:

    
mask = regionmask.defined_regions.srex.mask(rasm, lon_name="xc", lat_name="yc")
mask

We want to select the region 'NAS' (Northern Asia).

Select using `where`

We have to select by index (the number of the region), we thus map from the abbreviation to the index.



In [ ]:

    
rasm_NAS = rasm.where(mask == regionmask.defined_regions.srex.map_keys("NAS"))