query a polygon in an xarray object

[raybellwaves]

Let's say I have a global dataset

e.g.

import xarray as xr
url = "http://nomads.ncep.noaa.gov:80/dods/gfs_0p25/gfs20210215/gfs_0p25_00z" # YYYYMMDD
ds = xr.open_dataset(url)
da = ds.tmp2m.isel(time=0)
da = da.assign_coords(lon=(((ds.lon + 180) % 360) - 180))

and I would like to extract the data over my home state of Florida

[dcherian]

https://corteva.github.io/rioxarray/stable/examples/clip_geom.html ?

Also see regionmask: https://regionmask.readthedocs.io/en/stable/index.html but I don't see a crop or clip method.

[mathause]

In regionmask you cannot clip directly but you need to create a "mask" and then use the mask for the clipping.

import matplotlib.pyplot as plt
import numpy as np
import regionmask
import xarray as xr

%matplotlib

air = xr.tutorial.open_dataset("air_temperature")

# get the US states polygon
us_states = regionmask.defined_regions.natural_earth.us_states_50

# create a mask
mask_3D = us_states.mask_3D(air.air)

regionmask.plot_3D_mask(mask_3D)

# swap dim for indexing
mask_3D_a = mask_3D.swap_dims({"region": "abbrevs"})


plt.figure()
air.air.where(mask_3D_a.sel(abbrevs="FL")).isel(time=0).plot()


# now you can calculate the weighted average for all regions in one go
weights = np.cos(np.deg2rad(air.lat))

regional_means = air.air.weighted(weights * mask_3D).mean(("lat", "lon"))

plt.figure()
regional_means.plot()

[raybellwaves]

Thanks @dcherian and @mathause for your responses. And thanks @mathause for your work on region mask. As an FYI I got the rioxarray approach working here corteva/rioxarray#259