[gdal-dev] Re: Best way to test wheter point is contained in Polygon

[Jose Gomez-Dans]

I reply to myself, to see if someone has any insights.

2008/9/22 Jose Gomez-Dans <jgomezdans at gmail.com>

> Hi,
> I have a raster file and a vector file. I want to assign to each pixel in
> my raster size a value derived from the vector file (which are polygons). To
> do this, I want to test that the centroid of my pixel lies within a given
> polygon. In the past, I have done this gdal_rasterize-ing the polygon file,
> and going through numpy arrays. It is a bit cumbersome, and the
> rasterization is not really needed with OGR. However, it seems that if my
> pixel centroid is (x,y), I will need to go through all features in my layer,
> extract the geometry and do a Contains() test for my (x,y) point. I can cut
> the number of iterations (I guess, not tried) by using a SpatialFilter on my
> vector layer, so only one (ish) geometry is returned, so I would only need
> to do one Contains test (I do have a _lot_ of polygons in this file)
>

The code I have loops over each pixel, calculates its centroid from the
GeoTransform, and reprojects this to match the OGR's dataset projection. I
then produce a geometry for this centroid, and apply a spatial filter to my
OGR layer so that it encompasses said centroid. It is then a matter of
looping through the geometries until I get a match. The code works, but it
is quite slow. I have speeded it up substantially by creating an spatial
index on my shapefile, but it still takes a long time to go through this
(there are a lot of pixels to go through, and I need to do this for many
raster files!), so I was wondering if anyone has some suggestions on how to
improve this.

here's some relevant snippet of code, just in case someone's interested:
# I am assuming I have a raster_map, it's X and Y size, and its geotransform
# L is my ogr layer
#[...]
 for row in range(g.RasterXSize):
    for col in range(g.RasterYSize):
    if (raster_map[col,row] > 0): #Loop over each pixel, and if the value of
the pixel >0, try to find it's polygon.
        x = geotransform[0] + (row+0.5)*geotransform[1]
        y = geotransform[3] + (col+0.5)*geotransform[5]
        (lon,lat,z) = t.TransformPoint ( x,y ) # I need to transform the
data first.
        geop = ogr.CreateGeometryFromWkt("POINT (%f %f)"%(lon,lat)) # This
is the centroid of my pixel
        feat = L.GetNextFeature()
        filter_layer = L.SetSpatialFilterRect(lon-.5, lat-.5, lon+.5,lat+.5)
# To  minimise searches, do a spatial filter on the layer, so we only have a
handful of geometries to worry about
        while feat is not None:
            geo_grid = feat.GetGeometryRef()
            if geo_grid.Contains( geop ): # Found gridcell
                #Found a match. Do something with it, and stop searching
geometries.
                break
            feat.Destroy()
            feat = L.GetNextFeature()

        L.SetSpatialFilter(None) #Reset the spatial filter
        L.ResetReading() #Go back to the first feature of the layer

Cheers,
Jose

-- 
Centre for Terrestrial Carbon Dynamics
Department of Geography, University College London
Gower Street, London WC1E 6BT, UK
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080923/9a8d9dfb/attachment.html