rasterio and geopandas can be combined with xray to make converting shapefiles into raster masks pretty easy. Here's a quick demo:

``` python import geopandas from rasterio import features from affine import Affine

def transform_from_latlon(lat, lon): lat = np.asarray(lat) lon = np.asarray(lon) trans = Affine.translation(lon[0], lat[0]) scale = Affine.scale(lon[1] - lon[0], lat[1] - lat[0]) return trans * scale

def rasterize(shapes, coords, fill=np.nan, kwargs): """Rasterize a list of (geometry, fill_value) tuples onto the given xray coordinates. This only works for 1d latitude and longitude arrays. """ transform = transform_from_latlon(coords['latitude'], coords['longitude']) out_shape = (len(coords['latitude']), len(coords['longitude'])) raster = features.rasterize(shapes, out_shape=out_shape, fill=fill, transform=transform, dtype=float, kwargs) return xray.DataArray(raster, coords=coords, dims=('latitude', 'longitude'))

this shapefile is from natural earth data

http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-1-states-provinces/

states = geopandas.read_file('/Users/shoyer/Downloads/ne_10m_admin_1_states_provinces_lakes') us_states = states.query("admin == 'United States of America'").reset_index(drop=True) state_ids = {k: i for i, k in enumerate(us_states.woe_name)} shapes = [(shape, n) for n, shape in enumerate(us_states.geometry)]

ds = xray.Dataset(coords={'longitude': np.linspace(-125, -65, num=5000), 'latitude': np.linspace(50, 25, num=3000)}) ds['states'] = rasterize(shapes, ds.coords)

example of applying a mask

ds.states.where(ds.states == state_ids['California']).plot() ```

Once you have the rasterized geometries, you can use them as arrays to do arithmetic: https://github.com/xray/xray/issues/503

When we figure out how to represent coordinate reference systems properly in xray we might add in a direct wrapper for some of these rasterio functions.