Regionmask¶

By Mazvita Chikomo 11/20/2023¶

Creates Masks for geographic regions¶

When working with gridded data such as climate model output or reanalysis data it is often important to create regional averages, e.g., over countries, continents or regions defined in the literature. To do so we need to know for each grid point to which region it belongs to. RegionMask allows us to do just that!

Regionmask is a package that:¶

Can be used sed to create masks of geographic regions for arbitrary longitude and latitude grids.
These masks indicate which region a gridpoint belongs to.
They come in two varian-2D integer Masks-3D Boolean Masks.
takes great care to consistently treat gridpoints and overlapping regions,
contains a number of defined regions, including:
countries, landmasks, regions used in the scientific literature, can plot figures of these regions
supports using arbitrary existing or user-defined region definitions:regions defined as shapefiles can be accessed via geopandas, user-defined regions can be created via numpy or shapely

This tutorial will highlight some of the utilities of using the regionmask package with additional information at the end of the tutorial

References: https://regionmask.readthedocs.io/en/stable/index.html#s://github.com/hydpy-dev/hydpyub.com/hydpy-dev/hydpy

Where to start:¶

1. Installation Instructions¶

In this tutorial there are two ways to install regionmask:

i. One can use conda, which at time can be the easiest way to install, however in some cases you may run into this issue:

"Solving environment: failed with initial frozen solve. Retrying with flexible solve. Solving environment: / failed with repodata from current_repodata.json, will retry with next repodata source."

In [ ]:

Copied!

# Code to install using conda:
conda install -c conda-forge regionmask cartopy
# Code to install using conda:
conda install -c conda-forge regionmask cartopy

ii. Alternatively one can use pip to install the regionmask package, as well as the other packages used in this tutorial: cartopy.crs, xarray, numpy and matplotlib

In [ ]:

Copied!

#Code to install all required dependencies:
pip install regionmask
#Code to install all required dependencies:
pip install regionmask

2. Importing various packages required for this tutorial:¶

Useful Packages:

regionmask: Creates masks for different regions.

cartopy: Cartopy is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses.

xarray: An open source project and Python package that introduces labels in the form of dimensions, coordinates, and attributes on top of raw NumPy-like arrays, which allows for more intuitive, more concise, and less error-prone user experience.

numpy: A Python library that provides a multidimensional array object, various derived objects (such as masked arrays and matrices), and an assortment of routines for fast operations on arrays, including mathematical, logical, shape manipulation, sorting, selecting, I/O, discrete Fourier transforms, basic linear algebra, basic statistical operations, random simulation and much more.

matplotlib: a comprehensive library for creating plots and statistical visualizations in Python

In [8]:

Copied!





# Importing packages that will be necessary for this tutorial
import regionmask
import cartopy.crs as ccrs 
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
# Importing packages that will be necessary for this tutorial
import regionmask
import cartopy.crs as ccrs 
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe

Example 1: Creating a mask using predefined regions using defined_regions function¶

Using the function regionmask.defined_regions can be made from the outline of the countries obtained from Natural Earth. They are automatically downloaded, cached1 and opened with geopandas. The following countries and regions are defined in regionmask.

Countries 1:110m
Countries 1:50
Countries 1:10m
US States 1:50m
US States 1:10m

For this example SREX regions will be used as the predefined regions(Seneviratne et al., 2012)

In [55]:

Copied!

#Create an object with various predefined regions
srex = regionmask.defined_regions.srex

#Showing what srex looks like
print(srex)
#Create an object with various predefined regions
srex = regionmask.defined_regions.srex

#Showing what srex looks like
print(srex)

<regionmask.Regions 'SREX'>
Source:   Seneviratne et al., 2012 (https://www.ipcc.ch/site/assets/uploads/2...
overlap:  False

Regions:
 1 ALA       Alaska/N.W. Canada
 2 CGI     Canada/Greenl./Icel.
 3 WNA         W. North America
 4 CNA         C. North America
 5 ENA         E. North America
..  ..                      ...
22 EAS                  E. Asia
23 SAS                  S. Asia
24 SEA                S.E. Asia
25 NAU             N. Australia
26 SAU S. Australia/New Zealand

[26 regions]

Every region has two plotting functions that draws the outlines of all regions:¶

plot(): draws the region polygons on a cartopy GeoAxes (map)
plot_regions(): draws the the region polygons onlynes:

Visualizing the different regions using plot_regions()¶

Calling .plot_regions() on any region without any arguments draws the region polygons only

In [12]:

Copied!

srex.plot_regions()
srex.plot_regions()

Out[12]:

<AxesSubplot:>

No description has been provided for this image

Adding a map to the polygons using plot_regions and matplotlib.patheffects:¶

To create a map one can use the subplots and plot method from Matplotlib, which has a large number of arguments to adjust the layout of the axes. For example, you can pass a custom projection, the labels can display the abbreviation insead of the region number, one can add coastlines etc.

This example also shows how to use matplotlib.patheffects to ensure the labels are easily readable without covering too much of the map (compare to the map above):

import matplotlib.patheffect as pe (if you haven't already)

We will also use Cartopy to produce maps and other geospatial data analyses.

import cartopy.crs as ccrs (if you haven't already)

In [27]:

Copied!





#Import cartopy if it was not imported earlier
import cartopy.crs as ccrs

#Creating a figure (map canvas) with projection Robisnon
figure, ax = plt.subplots(subplot_kw=dict(projection=ccrs.Robinson()))

#Creating text labels for each region, but since path effects is being used the labels will not be in big gray boxes covering important information
text_kws = dict(
    bbox=dict(color="none"),
    path_effects=[pe.withStroke(linewidth=2, foreground="w")],
    color="#67000d",
    fontsize=8,)

# Plotting the predefined srex regions onto the map canvas above and adding the text labels
srex.plot_regions(ax=ax, line_kws=dict(lw=1), text_kws=text_kws)

#Adding coastlines to the map
ax.coastlines()
# Creating a global map
ax.set_global()
#Import cartopy if it was not imported earlier
import cartopy.crs as ccrs

#Creating a figure (map canvas) with projection Robisnon
figure, ax = plt.subplots(subplot_kw=dict(projection=ccrs.Robinson()))

#Creating text labels for each region, but since path effects is being used the labels will not be in big gray boxes covering important information
text_kws = dict(
    bbox=dict(color="none"),
    path_effects=[pe.withStroke(linewidth=2, foreground="w")],
    color="#67000d",
    fontsize=8,)

# Plotting the predefined srex regions onto the map canvas above and adding the text labels
srex.plot_regions(ax=ax, line_kws=dict(lw=1), text_kws=text_kws)

#Adding coastlines to the map
ax.coastlines()
# Creating a global map
ax.set_global()

Visualizing the different regions using plot()¶

Calling .plot() on any region without any arguments draws the default map with a PlateCarree() projection and includes the coastlines:

In [15]:

Copied!

srex.plot();
srex.plot();

Creating a more detailed global map using .plot() and matplotlib.patheffects¶

In [26]:

Copied!





#import matplotlib.patheffect as pe if not done already
import matplotlib.patheffects as pe

#The same text labels (text_kws) will be used as before 
ax = srex.plot(
    projection=ccrs.Robinson(), label="abbrev", add_ocean=True, text_kws=text_kws)
#Creating a global map
ax.set_global()
#import matplotlib.patheffect as pe if not done already
import matplotlib.patheffects as pe

#The same text labels (text_kws) will be used as before 
ax = srex.plot(
    projection=ccrs.Robinson(), label="abbrev", add_ocean=True, text_kws=text_kws)
#Creating a global map
ax.set_global()

Plotting only a subset of the regions:¶

To plot a selection of regions subset them using indexing:

In [29]:

Copied!





# regions can be selected by number, abbreviation or long name
regions = [11, "CEU", "S. Europe/Mediterranean"]

# choose a good projection for regional maps
projection = ccrs.LambertConformal(central_longitude=15)

#Creating the plot
ax = srex[regions].plot(
    add_ocean=True,
    resolution="50m",
    projection=projection,
    label="abbrev",
    text_kws=text_kws,)

# fine tune the extent
ax.set_extent([-15, 45, 28, 76], crs=ccrs.PlateCarree())
# regions can be selected by number, abbreviation or long name
regions = [11, "CEU", "S. Europe/Mediterranean"]

# choose a good projection for regional maps
projection = ccrs.LambertConformal(central_longitude=15)

#Creating the plot
ax = srex[regions].plot(
    add_ocean=True,
    resolution="50m",
    projection=projection,
    label="abbrev",
    text_kws=text_kws,)

# fine tune the extent
ax.set_extent([-15, 45, 28, 76], crs=ccrs.PlateCarree())

Example 2: Creating your own regions using the Regions Function¶

Assume you have two custom regions in the US, you can easily use these to create Regions:

In [53]:

Copied!





# Create two numpy arrays, one for each region
US1 = np.array([[-100.0, 30], [-100, 40], [-120, 35]])
US2 = np.array([[-100.0, 30], [-80, 30], [-80, 40], [-100, 40]])

#Use the Regions function to convert the new variables into regions 
regionmask.Regions([US1, US2])
# Create two numpy arrays, one for each region
US1 = np.array([[-100.0, 30], [-100, 40], [-120, 35]])
US2 = np.array([[-100.0, 30], [-80, 30], [-80, 40], [-100, 40]])

#Use the Regions function to convert the new variables into regions 
regionmask.Regions([US1, US2])

Out[53]:

<regionmask.Regions 'unnamed'>
overlap:  None

Regions:
0 r0 Region0
1 r1 Region1

[2 regions]

This creates two unnamed regions, but giving the regions additional information will make them more useful. To set names and abbrevs:

In [50]:

Copied!





# Creating a names variable for each region
names = ["US_west", "US_east"]
# Creating an abbreviation variable for each region
abbrevs = ["USw", "USe"]


# Setting the name and abbreviation for each region
USregions = regionmask.Regions([US1, US2], names=names, abbrevs=abbrevs, name="US")

#Calling the new variable USregions to see if the names and abbreviations have been set
USregions
# Creating a names variable for each region
names = ["US_west", "US_east"]
# Creating an abbreviation variable for each region
abbrevs = ["USw", "USe"]


# Setting the name and abbreviation for each region
USregions = regionmask.Regions([US1, US2], names=names, abbrevs=abbrevs, name="US")

#Calling the new variable USregions to see if the names and abbreviations have been set
USregions

Out[50]:

<regionmask.Regions 'US'>
overlap:  None

Regions:
0 USw US_west
1 USe US_east

[2 regions]

Plotting the assigned regions onto a map:¶

In [52]:

Copied!





#Plotting the regions variable USregions using the plot() method from matplotlib and assigning abbreviations as the label
ax = USregions.plot(label="abbrev")

#Zooming into the two regions
# fine tune the extent
ax.set_extent([225, 300, 25, 45], crs=ccrs.PlateCarree())
#Plotting the regions variable USregions using the plot() method from matplotlib and assigning abbreviations as the label
ax = USregions.plot(label="abbrev")

#Zooming into the two regions
# fine tune the extent
ax.set_extent([225, 300, 25, 45], crs=ccrs.PlateCarree())

Example 3. Using the Mask Function to create a mask to define points to each regions¶

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

In [56]:

Copied!





# Start by defining a lon/lat grid with 1 degree grid spacing, where the points define the center of the grid. 
lon = np.arange(-179.5, 180)
lat = np.arange(-89.5, 90)

#Creating the mask that uses the srex predefined regions with 1 degree grid spacing
mask = regionmask.defined_regions.srex.mask(lon, lat)

mask
# Start by defining a lon/lat grid with 1 degree grid spacing, where the points define the center of the grid. 
lon = np.arange(-179.5, 180)
lat = np.arange(-89.5, 90)

#Creating the mask that uses the srex predefined regions with 1 degree grid spacing
mask = regionmask.defined_regions.srex.mask(lon, lat)

mask

Out[56]:

<xarray.DataArray 'mask' (lat: 180, lon: 360)>
array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])
Coordinates:
  * lat      (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
  * lon      (lon) float64 -179.5 -178.5 -177.5 -176.5 ... 177.5 178.5 179.5
Attributes:
    standard_name:  region
    flag_values:    [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19...
    flag_meanings:  ALA CGI WNA CNA ENA CAM AMZ NEB WSA SSA NEU CEU MED SAH W...

xarray.DataArray

'mask'

lat: 180
lon: 360

nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

array([[nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       ...,
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan],
       [nan, nan, nan, ..., nan, nan, nan]])

Coordinates: (2)

lat

(lat)

float64

-89.5 -88.5 -87.5 ... 88.5 89.5

array([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
       -79.5, -78.5, -77.5, -76.5, -75.5, -74.5, -73.5, -72.5, -71.5, -70.5,
       -69.5, -68.5, -67.5, -66.5, -65.5, -64.5, -63.5, -62.5, -61.5, -60.5,
       -59.5, -58.5, -57.5, -56.5, -55.5, -54.5, -53.5, -52.5, -51.5, -50.5,
       -49.5, -48.5, -47.5, -46.5, -45.5, -44.5, -43.5, -42.5, -41.5, -40.5,
       -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5, -30.5,
       -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5,
       -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5,
        -9.5,  -8.5,  -7.5,  -6.5,  -5.5,  -4.5,  -3.5,  -2.5,  -1.5,  -0.5,
         0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
        10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,
        20.5,  21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,
        30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,
        40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5,
        50.5,  51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,  58.5,  59.5,
        60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,  67.5,  68.5,  69.5,
        70.5,  71.5,  72.5,  73.5,  74.5,  75.5,  76.5,  77.5,  78.5,  79.5,
        80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5])

lon

(lon)

float64

-179.5 -178.5 ... 178.5 179.5

array([-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5])

Indexes: (2)

lat

PandasIndex

PandasIndex(Index([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
       ...
        80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5],
      dtype='float64', name='lat', length=180))

lon

PandasIndex

PandasIndex(Index([-179.5, -178.5, -177.5, -176.5, -175.5, -174.5, -173.5, -172.5, -171.5,
       -170.5,
       ...
        170.5,  171.5,  172.5,  173.5,  174.5,  175.5,  176.5,  177.5,  178.5,
        179.5],
      dtype='float64', name='lon', length=360))

Attributes: (3)
standard_name :
region
flag_values :
[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26]
flag_meanings :
ALA CGI WNA CNA ENA CAM AMZ NEB WSA SSA NEU CEU MED SAH WAF EAF SAF NAS WAS CAS TIB EAS SAS SEA NAU SAU

- Mask is now an xarray.Dataset with shape (lat x lon).¶

- If a point lies outside of a region a value of NAN is given¶

In [38]:

Copied!





#Verifying the shape of mask
print('Lat shape is:', lat.shape)
print('Lon shape is:', lon.shape)
print('Mask shape is:', mask.shape)
#Verifying the shape of mask
print('Lat shape is:', lat.shape)
print('Lon shape is:', lon.shape)
print('Mask shape is:', mask.shape)

Lat shape is: (180,)
Lon shape is: (360,)
Mask shape is: (180, 360)

Visualizing the masks created using cartopy and matplotlib¶

In [54]:

Copied!





# Import if you have not done so before
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

#Creating the mapping canvas and setting the projection to PlateCarree
f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
#Adding Coastlines
ax.coastlines()

#Adding outlines for each srex region to see if the mask lines us
regionmask.defined_regions.srex.plot(ax=ax, add_label=False, line_kws=dict(lw=0.5, color="0.2"))

#Plotting the new masks
mask.plot(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)
# Import if you have not done so before
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

#Creating the mapping canvas and setting the projection to PlateCarree
f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
#Adding Coastlines
ax.coastlines()

#Adding outlines for each srex region to see if the mask lines us
regionmask.defined_regions.srex.plot(ax=ax, add_label=False, line_kws=dict(lw=0.5, color="0.2"))

#Plotting the new masks
mask.plot(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)

Out[54]:

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7f5c29e809d0>

Part 4: Applying regionmask to real data¶

For this tutorial we will use the NCEP reanalysis air temperature subsets which creates an air temperature field over North America. This data is freely available as an open dataset as a part of the Xarray tutorial library.

https://docs.xarray.dev/en/stable/generated/xarray.tutorial.open_dataset.html

In [80]:

Copied!

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

In [81]:

Copied!

#Displaying the data
airtemps
#Displaying the data
airtemps

Out[81]:

<xarray.Dataset>
Dimensions:  (lat: 25, time: 2920, lon: 53)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air      (time, lat, lon) float32 241.2 242.5 243.5 ... 296.5 296.2 295.7
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...

xarray.Dataset

Dimensions:
- lat: 25
- time: 2920
- lon: 53

Coordinates: (3)

lat

(lat)

float32

75.0 72.5 70.0 ... 20.0 17.5 15.0

standard_name :: latitude
long_name :: Latitude
units :: degrees_north
axis :: Y

array([75. , 72.5, 70. , 67.5, 65. , 62.5, 60. , 57.5, 55. , 52.5, 50. , 47.5,
       45. , 42.5, 40. , 37.5, 35. , 32.5, 30. , 27.5, 25. , 22.5, 20. , 17.5,
       15. ], dtype=float32)

lon

(lon)

float32

200.0 202.5 205.0 ... 327.5 330.0

standard_name :: longitude
long_name :: Longitude
units :: degrees_east
axis :: X

array([200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
       225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
       250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
       275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
       300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
       325. , 327.5, 330. ], dtype=float32)

time

(time)

datetime64[ns]

2013-01-01 ... 2014-12-31T18:00:00

standard_name :: time
long_name :: Time

array(['2013-01-01T00:00:00.000000000', '2013-01-01T06:00:00.000000000',
       '2013-01-01T12:00:00.000000000', ..., '2014-12-31T06:00:00.000000000',
       '2014-12-31T12:00:00.000000000', '2014-12-31T18:00:00.000000000'],
      dtype='datetime64[ns]')

Data variables: (1)

air

(time, lat, lon)

float32

241.2 242.5 243.5 ... 296.2 295.7

long_name :: 4xDaily Air temperature at sigma level 995
units :: degK
precision :: 2
GRIB_id :: 11
GRIB_name :: TMP
var_desc :: Air temperature
dataset :: NMC Reanalysis
level_desc :: Surface
statistic :: Individual Obs
parent_stat :: Other
actual_range :: [185.16 322.1 ]

array([[[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    ,
         238.59999],
        [243.79999, 244.5    , 244.7    , ..., 232.79999, 235.29999,
         239.29999],
        [250.     , 249.79999, 248.89   , ..., 233.2    , 236.39   ,
         241.7    ],
        ...,
        [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    ,
         294.69998],
        [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    ,
         295.19998],
        [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   ,
         296.6    ]],

       [[242.09999, 242.7    , 243.09999, ..., 232.     , 233.59999,
         235.79999],
        [243.59999, 244.09999, 244.2    , ..., 231.     , 232.5    ,
         235.7    ],
        [253.2    , 252.89   , 252.09999, ..., 230.79999, 233.39   ,
         238.5    ],
...
        [293.69   , 293.88998, 295.38998, ..., 295.09   , 294.69   ,
         294.29   ],
        [296.29   , 297.19   , 297.59   , ..., 295.29   , 295.09   ,
         294.38998],
        [297.79   , 298.38998, 298.49   , ..., 295.69   , 295.49   ,
         295.19   ]],

       [[245.09   , 244.29   , 243.29   , ..., 241.68999, 241.48999,
         241.79   ],
        [249.89   , 249.29   , 248.39   , ..., 239.59   , 240.29   ,
         241.68999],
        [262.99   , 262.19   , 261.38998, ..., 239.89   , 242.59   ,
         246.29   ],
        ...,
        [293.79   , 293.69   , 295.09   , ..., 295.29   , 295.09   ,
         294.69   ],
        [296.09   , 296.88998, 297.19   , ..., 295.69   , 295.69   ,
         295.19   ],
        [297.69   , 298.09   , 298.09   , ..., 296.49   , 296.19   ,
         295.69   ]]], dtype=float32)

Indexes: (3)

lat

PandasIndex

PandasIndex(Index([75.0, 72.5, 70.0, 67.5, 65.0, 62.5, 60.0, 57.5, 55.0, 52.5, 50.0, 47.5,
       45.0, 42.5, 40.0, 37.5, 35.0, 32.5, 30.0, 27.5, 25.0, 22.5, 20.0, 17.5,
       15.0],
      dtype='float32', name='lat'))

lon

PandasIndex

PandasIndex(Index([200.0, 202.5, 205.0, 207.5, 210.0, 212.5, 215.0, 217.5, 220.0, 222.5,
       225.0, 227.5, 230.0, 232.5, 235.0, 237.5, 240.0, 242.5, 245.0, 247.5,
       250.0, 252.5, 255.0, 257.5, 260.0, 262.5, 265.0, 267.5, 270.0, 272.5,
       275.0, 277.5, 280.0, 282.5, 285.0, 287.5, 290.0, 292.5, 295.0, 297.5,
       300.0, 302.5, 305.0, 307.5, 310.0, 312.5, 315.0, 317.5, 320.0, 322.5,
       325.0, 327.5, 330.0],
      dtype='float32', name='lon'))

time

PandasIndex

PandasIndex(DatetimeIndex(['2013-01-01 00:00:00', '2013-01-01 06:00:00',
               '2013-01-01 12:00:00', '2013-01-01 18:00:00',
               '2013-01-02 00:00:00', '2013-01-02 06:00:00',
               '2013-01-02 12:00:00', '2013-01-02 18:00:00',
               '2013-01-03 00:00:00', '2013-01-03 06:00:00',
               ...
               '2014-12-29 12:00:00', '2014-12-29 18:00:00',
               '2014-12-30 00:00:00', '2014-12-30 06:00:00',
               '2014-12-30 12:00:00', '2014-12-30 18:00:00',
               '2014-12-31 00:00:00', '2014-12-31 06:00:00',
               '2014-12-31 12:00:00', '2014-12-31 18:00:00'],
              dtype='datetime64[ns]', name='time', length=2920, freq=None))

Attributes: (5)
Conventions :
COARDS
title :
4x daily NMC reanalysis (1948)
description :
Data is from NMC initialized reanalysis (4x/day). These are the 0.9950 sigma level values.
platform :
Model
references :
http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html

Visualizing what the first time step in the dataset looks like using cartopy and indexing¶

In [86]:

Copied!





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

#Creating the map canvas
ax = plt.subplot(111, projection=proj)

#Plotting the first time step an dtransforming the data projection from PlateCarree to the map canvas projection
airtemps.isel(time=1).air.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())
#Adding coastline
ax.coastlines();
# choose a good projection for regional maps
proj = ccrs.LambertConformal(central_longitude=-100)

#Creating the map canvas
ax = plt.subplot(111, projection=proj)

#Plotting the first time step an dtransforming the data projection from PlateCarree to the map canvas projection
airtemps.isel(time=1).air.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree())
#Adding coastline
ax.coastlines();

Creating a mask that matches the srex predefined regions to the air temperature data¶

In [83]:

Copied!

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

Region mask is able to automatically detect whether 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 [87]:

Copied!

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

Grid extent:    200°E to 330°E
Region extent: -168°E to 180°E

Now plotting the mask over the data¶

In [88]:

Copied!





# Creating the map canvas
ax = plt.subplot(111, projection = ccrs.AlbersEqualArea(central_lon, central_lat))

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

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

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


# for colorbar: find abbreviations of all regions that were selected

#1. Assigning a unique value for each region
reg = np.unique(mask.values)
#2. Ensuring to only collect values for areas within a mask (isnan Tests element-wise for NaN and return result as a boolean array)
reg = reg[~np.isnan(reg)]
#3. Assigning the abbreviation of each srex region to the corresponding mask regions 
abbrevs = regionmask.defined_regions.srex[reg].abbrevs

#Assigning the colorbar to the figure
cbar = plt.colorbar(h, orientation="horizontal", fraction=0.075, pad=0.05)

#Editing the colorbar and adding ticks to go with each region
cbar.set_ticks(reg)
cbar.set_ticklabels(abbrevs)
cbar.set_label("Region")

#Adding coastlines to the map canvas
ax.coastlines()

#Zooming into the area and fine tune the extent
ax.set_extent([200, 330, 10, 75], crs=ccrs.PlateCarree())
# Creating the map canvas
ax = plt.subplot(111, projection = ccrs.AlbersEqualArea(central_lon, central_lat))

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

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

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


# for colorbar: find abbreviations of all regions that were selected

#1. Assigning a unique value for each region
reg = np.unique(mask.values)
#2. Ensuring to only collect values for areas within a mask (isnan Tests element-wise for NaN and return result as a boolean array)
reg = reg[~np.isnan(reg)]
#3. Assigning the abbreviation of each srex region to the corresponding mask regions 
abbrevs = regionmask.defined_regions.srex[reg].abbrevs

#Assigning the colorbar to the figure
cbar = plt.colorbar(h, orientation="horizontal", fraction=0.075, pad=0.05)

#Editing the colorbar and adding ticks to go with each region
cbar.set_ticks(reg)
cbar.set_ticklabels(abbrevs)
cbar.set_label("Region")

#Adding coastlines to the map canvas
ax.coastlines()

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

Additional Information¶

Useful API information from the regionmask package:¶

Top Level Functions¶

mask_geopandas(geodataframe, lon_or_obj[, ...]): Creates a 2D float mask of a set of regions for the given lat/ lon grid
(geodataframe, lon_or_obj): Creates a 3D boolean mask of a set of regions for the given lat/ lon grid
from_geopandas(geodataframe, *[, numbers, ...]): Creates regionmask.Regions from a geopandas.GeoDataFrame.
flatten_3D_mask(mask_3D):Flattens 3D masks
plot_3D_mask(mask_3D, **kwargs): Flattens and plot 3D masks
set_options(**kwargs):Sets options for regionmask in a controlled context.
get_options():Gets options for regionmask.et options for regionmask.

Source: https://regionmask.readthedocs.io/en/stable/api.html

References:¶

https://regionmask.readthedocs.io/en/stable/notebooks/plotting.html

https://matplotlib.org/

https://docs.xarray.dev/en/stable/

https://numpy.org/doc/stable/

https://www.ibm.com/docs/en/watson-studio-local/1.2.3?topic=notebooks-markdown-jupyter-cheatsheet

https://docs.xarray.dev/en/stable/generated/xarray.tutorial.open_dataset.html https://numpy.org/doc/stable/reference/generated/numpy.isnan.html

In [ ]:

Copied!