# [SciPy-User] ROI, the sequel

Joe Kington jkington@wisc....
Wed Aug 4 12:56:34 CDT 2010

```For whatever it's worth, matplotlib has a fairly fast function
matplotlib.nxutils.points_inside_poly that will do this for polygons with no
"holes".

To use it to rasterize a polygon region of interest, just build a 2xN array
of the x,y locations of the verticies of your grid, pass it to to
points_inside_poly along with the polygon verticies, and then reshape the
output boolean mask to the size of your original grid.

Hope that helps a bit,
-Joe

2010/8/4 Thøger Emil Juul Thorsen <thoeger@fys.ku.dk>

> Thanks to all of your for your suggestions - I think I'm going to try
> out the gdal solution first.
>
> But where can I request a similar function for includion in NumPy? I
> cannot be the only one needing to work with non-trivial Regions Of
> Interest on large arrays?
>
> Best;
>
> Emil
>
>
>
> On Wed, 2010-08-04 at 13:07 +1000, Pinner, Luke wrote:
> > GDAL/OGR perhaps? www.gdal.org
> >
> > Some semi-pseudo code below.
> >
> > from osgeo import gdal,gdal_array,ogr
> > import numpy
> >
> > # Open a polygon layer to rasterize from.
> > poly_ds = ogr.Open( some file dataset ) #Can also create an in-memory
> poly layer from scratch
> > poly_lyr=poly_ds.GetLayer(0)
> >
> > # Create an in-memory raster to rasterize into
> > target_ds = gdal.GetDriverByName('MEM').Create( approriate args...)
> >
> > #Rasterize...
> > gdal.RasterizeLayer( appropriate args... )
> > numpy_array = gdal_aray.DatasetReadAsArray(target_ds)
> >
> > Then use the numpy_array as you lie.
> >
> > Luke
> >
> > -----Original Message-----
> > From: scipy-user-bounces@scipy.org [mailto:scipy-user-bounces@scipy.org]
> On Behalf Of Thøger Emil Juul Thorsen
> > Sent: Wednesday, 4 August 2010 10:05 AM
> > To: scipy-user@scipy.org
> > Subject: [SciPy-User] ROI, the sequel
> >
> > Hello list;
> >
> > I need a tool for polygon clipping of 2d NumPy arrays, but I cannot seem
> > to find anything. Shapely will only export the vertices of the polygon.
> > NXutils and patches in matplotlib seem to be purely visual, in any case
> > I cannot see how they would be saved as arrays in data coordinates.
> >
> > There is a very neat IDL routine, POYFILLV which returns a 1d array of
> > the (1d) indices of the pixels that are inside the polygon, sort of like
> > a "where" function. Is there anything like it out there for numpy
> > arrays?
> >
> > Best;
> >
> > Emil
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> >
> > ------
> > If you have received this transmission in error please notify us
> immediately by return e-mail and delete all copies. If this e-mail or any
> attachments have been sent to you in error, that error does not constitute
> waiver of any confidentiality, privilege or copyright in respect of
> information in the e-mail or attachments.
> >
> >
> >
> > Please consider the environment before printing this email.
> >
> > ------
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100804/7d9e5edb/attachment.html
```