<jupyter_start><jupyter_text>Create xarray region maskIn this tutorial we will show how to create a mask for arbitrary latitude and longitude grids using xarray. It is very similar to the tutorial Create Mask (numpy). Import regionmask and check the version:<jupyter_code>import regionmask
regionmask.__version__<jupyter_output><empty_output><jupyter_text>Load xarray and the tutorial data:<jupyter_code>import xarray as xr
import numpy as np
airtemps = xr.tutorial.load_dataset('air_temperature')<jupyter_output><empty_output><jupyter_text>The example data is a temperature field over North America. Let's plot the first time step:<jupyter_code># 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();<jupyter_output><empty_output><jupyter_text>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 longituda 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. Here we use the Giorgi regions.The extent of the grid and the regions is not the same:<jupyter_code>lon = airtemps.lon
print("Grid extent:    ", lon.values.min(), lon.values.max())

bounds = regionmask.defined_regions.giorgi.bounds_global
print("Region extent: ", bounds[0], bounds[2])<jupyter_output><empty_output><jupyter_text>.. 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.<jupyter_code>mask = regionmask.defined_regions.giorgi.mask(airtemps)<jupyter_output><empty_output><jupyter_text>This did not work in earlier versions. We had to set `wrap_lon=True`. Let's plot the mask of the regions:<jupyter_code>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.giorgi[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());<jupyter_output><empty_output><jupyter_text>We want to select the region 'Central North America'. Thus we first need to find out which number this is:<jupyter_code>regionmask.defined_regions.giorgi.map_keys('Central North America')<jupyter_output><empty_output><jupyter_text>Select using `where` `xarray` provides the handy `where` function:<jupyter_code>airtemps_CNA = airtemps.where(mask == 6)<jupyter_output><empty_output><jupyter_text>Check everything went well by repeating the first plot with the selected region:<jupyter_code># choose a good projection for regional maps
proj=ccrs.LambertConformal(central_longitude=-100)

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

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

ax.coastlines();<jupyter_output><empty_output><jupyter_text>Looks good - let's take the area average and plot the time series. (Note: you should use `cos(lat)` weights to correctly calculate an area average. Unfortunately this is not yet (as of version 0.13) implemented in xarray.)<jupyter_code>ts_airtemps_CNA = airtemps_CNA.mean(dim=('lat', 'lon')) - 273.15
ts_airtemps = airtemps.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')

plt.legend(ncol=2);<jupyter_output><empty_output><jupyter_text>Select using `groupby`<jupyter_code># you can group over all integer values of the mask
# you have to take the mean over `stacked_lat_lon`
airtemps_all = airtemps.groupby(mask).mean('stacked_lat_lon')<jupyter_output><empty_output><jupyter_text>we can add the abbreviations and names of the regions to the DataArray<jupyter_code># extract the abbreviations and the names of the regions from regionmask
abbrevs = regionmask.defined_regions.giorgi[airtemps_all.region.values].abbrevs
names = regionmask.defined_regions.giorgi[airtemps_all.region.values].names

airtemps_all.coords['abbrevs'] = ('region', abbrevs)
airtemps_all.coords['names'] = ('region', names)
airtemps_all<jupyter_output><empty_output><jupyter_text>now we can select the regions in many ways<jupyter_code># as before, by the index of the region
r1 = airtemps_all.sel(region=6).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 == 'Eastern North America')).air.squeeze()

regs = [r1, r2, r3]<jupyter_output><empty_output><jupyter_text>Now, let's plot the three selected regions:<jupyter_code>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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11142857142857143, 'fraction_numerical': 0.013928571428571429, 'mean_word_length': 6.714876033057851, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '645950', 'license': 'MIT', 'n_tokens_mistral': 2038, 'n_tokens_neox': 1921, 'n_words': 508, 'path': 'docs/notebooks/mask_xarray.ipynb', 'repo_name': 'COVID-Weather/regionmask', 'size': '5600'}