<div dir="ltr">I reply to myself, to see if someone has any insights.<br><br><div class="gmail_quote">2008/9/22 Jose Gomez-Dans <span dir="ltr"><<a href="mailto:jgomezdans@gmail.com">jgomezdans@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><div dir="ltr">Hi,<br>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)</div></blockquote><div><br>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.<br>
<br>here's some relevant snippet of code, just in case someone's interested:<br># I am assuming I have a raster_map, it's X and Y size, and its geotransform<br># L is my ogr layer<br>#[...]<br> for row in range(g.RasterXSize):<br>
for col in range(g.RasterYSize):<br> if (raster_map[col,row] > 0): #Loop over each pixel, and if the value of the pixel >0, try to find it's polygon.<br> x = geotransform[0] + (row+0.5)*geotransform[1]<br>
y = geotransform[3] + (col+0.5)*geotransform[5]<br> (lon,lat,z) = t.TransformPoint ( x,y ) # I need to transform the data first.<br> geop = ogr.CreateGeometryFromWkt("POINT (%f %f)"%(lon,lat)) # This is the centroid of my pixel<br>
feat = L.GetNextFeature()<br> 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<br>
while feat is not None:<br> geo_grid = feat.GetGeometryRef()<br> if geo_grid.Contains( geop ): # Found gridcell<br> #Found a match. Do something with it, and stop searching geometries. <br>
break<br> feat.Destroy()<br> feat = L.GetNextFeature()<br> <br> L.SetSpatialFilter(None) #Reset the spatial filter<br> L.ResetReading() #Go back to the first feature of the layer<br>
</div></div><br>Cheers,<br>Jose<br clear="all"><br>-- <br>Centre for Terrestrial Carbon Dynamics<br>Department of Geography, University College London<br>Gower Street, London WC1E 6BT, UK<br>
</div>