Outline the cells in DataArray that match a simple, returning a Shapely Polygon

[Brewster Malevich]

Hello everyone,

I need to do something and I feel like there should be an existing way to do this in Python. Hopefully you all can point me in the right direction.

Let's say I have a 2D DataArray (lat, lon) that represents a climate field spanning the globe. The values in the DataArray are integers or strings that classify all points into a couple of categories. For example, we might have a nice sized blob of "12" values over the Indian Ocean. To the north of this blob, we might have a cluster of points with value "13", and so forth... All points are classified.

Here is the part where I could use  some help.

I'm trying to find a way of turning these segments DataArray into Shapely Polygons. I'm looking for a Shapely Polygon (or MultiPolygon) that contains all and only the "12" values. Another Polygon that contains only the "13" values, and so forth...

I can get a sloppy Polygon by by turning the DataArray cells into Shapely Points and then creating a convex hull. This is no good because the shape is too simple. I can also kill a mosquito with a cannon by using Delaunay triangulation to create a tighter polygon, but this is a better method for irregular points that are not cells on a regular grid. This seems like it should be a fairly common raster query -> vector polygon operation. Just outline the cells in the DataArray that match a dumb query, returning a Shapely Polygon.

Is there an existing way to easily do this in Python? Any Python GIS experts here? I'd appreciate it if anyone can point me to a library, book, method, function, or anything, really.

If need be, I can draw up a simple code example to show what I'm talking about.

Thanks much for your help!

[Fabien]

Hi,

Le jeudi 13 octobre 2016 22:23:49 UTC+2, Brewster Malevich a écrit :

I'm trying to find a way of turning these segments DataArray into Shapely Polygons. I'm looking for a Shapely Polygon (or MultiPolygon) that contains all and only the "12" values. Another Polygon that contains only the "13" values, and so forth...

 

I know two ways to do this. One is based on rasterio: https://mapbox.github.io/rasterio/topics/features.html (recommended), the other is based on matplotib's undocumented _cntr module. I've used it in one of our projects here (but my use case is narrower than yours).

Note that it is more a general problem than an "xarray" one, so I'm not sure if the xarray devs will be willing to implement it in xarray directly.

Cheers,

Fabien

[Brewster Malevich]

Hey Fab,

Thanks for the reply! This is very helpful.

You're right. I don't expect xarray devs to implement anything for this. I was just hoping that someone had dealt with this problem before. In retrospect I should have posted this to stack exchange.

Thanks again for your help,