@liyaojun have you tried the website http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-1-states-provinces/ ? It still works for me.

Here's a notebook that includes a working example: https://gist.github.com/shoyer/0eb96fa8ab683ef078eb

To get the facetted plot to work, you'll need to be running the development version of xray.

Hmm. My code does seem to have a bug when given a dataset with more than 2 dimensions. I'll see if I can fix this up and add it to the examples page on the docs.

On Mon, Sep 28, 2015 at 4:20 PM, slharris notifications@github.com wrote:

Thank you for fixing the last line, I can get your example to run with my own data without any errors being raised however the output does not appear to have any masked applied! I can use where() to plot a selected area but when it comes to being used with another dataset it doesn't appear to be working for me. Is there an example that shows something like temperature data over time for different states, extracted using the method above? This method will be so useful for me - if I can get it to work!

Any feedback will be greatly appreciated.

Reply to this email directly or view it on GitHub: https://github.com/xray/xray/issues/501#issuecomment-143899899

@slharris I think you need to modify your second to last line like so:

monthlymean_california = monthlymean.where(ds.states == state_ids['California']).groupby('time').mean()

@slharris I just updated the earlier example to show how to use where for masking.

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.

We'll probably release v0.6 in several weeks, once we have the key functionality we want for new the plotting module.

On Wed, Aug 5, 2015 at 7:24 PM, slharris notifications@github.com wrote:

This example is very helpful. Thank you.

I think the where() method you refer to will be very useful - When will the version 0.5.3 be released?

Thanks Sarah

On 31 July 2015 at 06:00, Stephan Hoyer notifications@github.com wrote:

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, *_kwargs) return xray.DataArray(raster, coords=coords, dims=('latitude', 'longitude'))

states =

geopandas.read_file('/Users/shoyer/Downloads/ne_10m_admin_1_states_provinces_lakes') geometries = states.query("admin == 'United States of America'").geometry shapes = [(shape, n) for n, shape in enumerate(geometries)]

ds = xray.Dataset(coords={'longitude': np.linspace(-125, -65, num=2000), 'latitude': np.linspace(50, 25, num=1000)}) ds['states'] = rasterize(shapes, ds.coords) plotting requires the dev version of xray

ds.states.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.

— Reply to this email directly or view it on GitHub https://github.com/xray/xray/issues/501#issuecomment-126461466.

Dr Sarah Harris Research Fellow School of Earth, Atmosphere and Environment Room 227, 9 Rainforest Walk Monash University Ph: 03 9902 4243 Email: sarah.harris@monash.edu Email%3Asarah.harris@monash.edu

— Reply to this email directly or view it on GitHub https://github.com/xray/xray/issues/501#issuecomment-128216372.