Show/Hide data repr

.. note:: This tutorial was generated from an IPython notebook that can be downloaded `here <../../notebooks/mask_2D.ipynb>`_. .. _mask_2D: Create 2D region masks ====================== In this tutorial we will show how to create 2D mask for arbitrary latitude and longitude grids. Import regionmask and check the version: .. code:: python import regionmask regionmask.__version__ .. parsed-literal:: '0.5.0+dev' Load xarray and the tutorial data: .. code:: python import xarray as xr import numpy as np Creating a mask --------------- We define a lon/ lat grid with a 1° grid spacing, where the points define the center of the grid. .. code:: python lon = np.arange(-179.5, 180) lat = np.arange(-89.5, 90) We create a mask with the SREX regions (Seneviratne et al., 2012). .. code:: python regionmask.defined_regions.srex .. parsed-literal:: Name: SREX Source: Seneviratne et al., 2012 (https://www.ipcc.ch/site/assets/uploads/2... 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] The function ``mask`` determines which gripoints lie within the polygon making up the each region: .. code:: python mask = regionmask.defined_regions.srex.mask(lon, lat) mask .. raw:: html

xarray.DataArray

'region'

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

Attributes: (0)

``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``: .. code:: python 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); .. image:: mask_2D_files/mask_2D_13_0.png 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: .. code:: python airtemps = xr.tutorial.load_dataset("air_temperature") The example data is a temperature field over North America. Let’s plot the first time step: .. code:: python # load plotting libraries import matplotlib.pyplot as plt import cartopy.crs as ccrs # 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(); .. image:: mask_2D_files/mask_2D_17_0.png 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. .. code:: python 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: .. code:: python lon = airtemps.lon print( "Grid extent: {:3.0f}°E to {:3.0f}°E".format(lon.values.min(), lon.values.max()) ) bounds = regionmask.defined_regions.srex.bounds_global print("Region extent: {:3.0f}°E to {:3.0f}°E".format(bounds[0], bounds[2])) .. parsed-literal:: Grid extent: 200°E to 330°E Region extent: -168°E to 180°E Let’s plot the mask of the regions: .. code:: python 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()) .. image:: mask_2D_files/mask_2D_23_0.png We want to select the region ‘Central North America’. Thus we first need to find out which number this is: .. code:: python CNA_index = regionmask.defined_regions.srex.map_keys("C. North America") CNA_index .. parsed-literal:: 4 Select using ``where`` ~~~~~~~~~~~~~~~~~~~~~~ ``xarray`` provides the handy ``where`` function: .. code:: python airtemps_CNA = airtemps.where(mask == CNA_index) Check everything went well by repeating the first plot with the selected region: .. code:: python # 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_ocean=False, add_label=False) airtemps_CNA.isel(time=1).air.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree()) ax.coastlines(); .. image:: mask_2D_files/mask_2D_30_0.png Looks good - let’s take the area average and plot the time series. 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 .. code:: python weights = np.cos(np.deg2rad(airtemps.lat)) ts_airtemps_CNA = airtemps_CNA.weighted(weights).mean(dim=("lat", "lon")) - 273.15 ts_airtemps = airtemps.weighted(weights).mean(dim=("lat", "lon")) - 273.15 # and the line plot f, ax = plt.subplots() ts_airtemps_CNA.air.plot.line(ax=ax, label="Central North America") ts_airtemps.air.plot(ax=ax, label="Entire Domain") ax.axhline(0, color="0.1", lw=0.5) plt.legend(ncol=2); .. image:: mask_2D_files/mask_2D_32_0.png To get the regional average over each region you need to loop over them. Select using ``groupby`` ~~~~~~~~~~~~~~~~~~~~~~~~ .. warning:: Using ``groupby`` offers some convinience and is faster than using ``where`` and a loop. However, it cannot currently be combinded with ``weighted`` (xarray `GH3937 `_). Therefore I recommend using ``where``. .. code:: python # you can group over all integer values of the mask airtemps_all = airtemps.groupby(mask).mean() we can add the abbreviations and names of the regions to the DataArray .. code:: python # extract the abbreviations and the names of the regions from regionmask abbrevs = regionmask.defined_regions.srex[airtemps_all.region.values].abbrevs names = regionmask.defined_regions.srex[airtemps_all.region.values].names airtemps_all.coords["abbrevs"] = ("region", abbrevs) airtemps_all.coords["names"] = ("region", names) now we can select the regions in many ways .. code:: python # as before, by the index of the region r1 = airtemps_all.sel(region=CNA_index).air # with the abbreviation r2 = airtemps_all.isel(region=(airtemps_all.abbrevs == "WNA")).air.squeeze() # with the long name r3 = airtemps_all.isel(region=(airtemps_all.names == "E. North America")).air.squeeze() regs = [r1, r2, r3] Now, let’s plot the three selected regions: .. code:: python f, axes = plt.subplots(3, 1, sharex=True) for i, reg in enumerate(regs): ax = axes[i] reg.plot(ax=ax) ax.set_title(reg.names.values) plt.setp(axes, xlabel="") f.tight_layout() .. image:: mask_2D_files/mask_2D_41_0.png 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: .. code:: python rasm = xr.tutorial.load_dataset("rasm") The example data is a temperature field over the Northern Hemisphere. Let’s plot the first time step: .. code:: python # load plotting libraries import matplotlib.pyplot as plt import cartopy.crs as ccrs # 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], add_ocean=False, coastlines=False, label="abbrev" ) ax.set_extent([-180, 180, 43, 90], ccrs.PlateCarree()) ax.coastlines(); .. image:: mask_2D_files/mask_2D_45_0.png Again we pass the xarray object to regionmask. We have to specify ``"xc"`` and ``"yc"`` as the longitude and latitude coordinates of the array: .. code:: python mask = regionmask.defined_regions.srex.mask(rasm, lon_name="xc", lat_name="yc") mask .. raw:: html

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. .. code:: python rasm_NAS = rasm.where(mask == regionmask.defined_regions.srex.map_keys("NAS")) Check everything went well by repeating the first plot with the selected region: .. code:: python # choose a projection proj = ccrs.NorthPolarStereo() ax = plt.subplot(111, projection=proj) ax.set_global() rasm_NAS.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=["NAS"], add_ocean=False, coastlines=False, label="abbrev" ) ax.set_extent([-180, 180, 45, 90], ccrs.PlateCarree()) ax.coastlines(); .. image:: mask_2D_files/mask_2D_51_0.png References ---------- - Special Report on Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation (SREX, Seneviratne et al., 2012: https://www.ipcc.ch/site/assets/uploads/2018/03/SREX-Ch3-Supplement_FINAL-1.pdf)