<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 TRANSITIONAL//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; CHARSET=UTF-8">
<META NAME="GENERATOR" CONTENT="GtkHTML/3.32.2">
</HEAD>
<BODY>
Frank, thanks for the clarification on iterating over the selection. I was able to run my script and in the simple case of returning one polygon feature it works great. However, I ran into some problems when I tried to assume that the filter would return multiple features.<BR>
<BR>
Is it possible to iterate over a selection and add every geometry object to a list, in other words can geometries be handled like other Python objects? I would later like to iterate over the list of geometries and verify whether they contain the centroids I'm interested in. My procedure can be summarized as follows (sample scripts at the end):<BR>
1. Select polygon(s) based on a SQL filter<BR>
2. Iterate over polygons from 1. and add geometries to a list<BR>
3. Iterate over a mesh of polygons, grab their centroids and verify whether they are within the geometries from 2.<BR>
<BR>
Script that fails (assumes I'm getting more than one geometry back and adds it to a list):<BR>
#!/usr/bin/python<BR>
<BR>
from osgeo import ogr<BR>
import sys, math, time, os<BR>
<BR>
grid = 'model1.shp'<BR>
gridDS = ogr.Open(grid, 1)<BR>
gridLayer = gridDS.GetLayer(0)<BR>
<BR>
aquifer = '/data/romania/judeteEPSG3844.shp'<BR>
aqDS = ogr.Open(aquifer)<BR>
aqLayer = aqDS.GetLayerByIndex(0)<BR>
aqLayer.SetAttributeFilter('DENJUD = "CLUJ"')<BR>
aqfeat = aqLayer.GetNextFeature()<BR>
aqgeom = aqfeat.GetGeometryRef()<BR>
geomList = [aqgeom]<BR>
<BR>
while aqfeat is not None: #loop thru features in case of islands/multiple polygons<BR>
aqfeat = aqLayer.GetNextFeature()<BR>
if aqfeat is not None:<BR>
geomList.append(aqfeat.GetGeometryRef())<BR>
<BR>
for i in range(gridLayer.GetFeatureCount()):<BR>
feat = gridLayer.GetFeature(i)<BR>
print feat.GetField('CELL_ID') #see how many features processed before segfault, usually first two features<BR>
geom = feat.GetGeometryRef()<BR>
point = geom.Centroid()<BR>
for polygeom in geomList:<BR>
if polygeom.Contains(point):<BR>
feat.SetField('IBOUND', 1)<BR>
gridLayer.SetFeature(feat)<BR>
break #break out of "for" if any poly contains the centroid<BR>
<BR>
Script that works (assumes I'm only getting one geometry back, no need for list):<BR>
#!/usr/bin/python<BR>
<BR>
from osgeo import ogr<BR>
import sys, math, time, os<BR>
<BR>
grid = 'model1.shp'<BR>
gridDS = ogr.Open(grid, 1)<BR>
gridLayer = gridDS.GetLayer(0)<BR>
<BR>
aquifer = '/data/romania/judeteEPSG3844.shp'<BR>
aqDS = ogr.Open(aquifer)<BR>
aqLayer = aqDS.GetLayerByIndex(0)<BR>
aqLayer.SetAttributeFilter('DENJUD = "CLUJ"')<BR>
aqfeat = aqLayer.GetNextFeature()<BR>
aqgeom = aqfeat.GetGeometryRef()<BR>
<BR>
for i in range(gridLayer.GetFeatureCount()):<BR>
feat = gridLayer.GetFeature(i)<BR>
geom = feat.GetGeometryRef()<BR>
point = geom.Centroid()<BR>
if aqgeom.Contains(point):<BR>
feat.SetField('IBOUND', 1)<BR>
gridLayer.SetFeature(feat)<BR>
<BR>
<BR>
-marius<BR>
<BR>
On Mon, 2011-07-11 at 22:01 -0700, Frank Warmerdam wrote:
<BLOCKQUOTE TYPE=CITE>
<PRE>
Marius,
I think a bug is called for.
On Mon, Jul 11, 2011 at 8:00 PM, Marius Jigmond
<<A HREF="mailto:mariusjigmond@hotmail.com">mariusjigmond@hotmail.com</A>> wrote:
> After some more investigation that is likely NOT the issue. I have an
> ExecuteSQL statement which selects a certain polygon based on an attribute
> value. Unfortunately it seems to return the wrong feature. The feature I
> query for is unique so a duplicate is out of the question. Here's the code:
>
> #!/usr/bin/python
>
> from osgeo import ogr
> import sys, math, time, os
>
> aquifer = '/data/romania/judeteEPSG3844.shp'
> aqDS = ogr.Open(aquifer)
> sql = 'select * from judeteEPSG3844 where DENJUD = "CLUJ"'
> aqLayer = aqDS.ExecuteSQL(sql)
> feat = aqLayer.GetFeature(0)
> print aqLayer.GetFeatureCount()
> print feat.GetField('DENJUD')
> aqDS.ReleaseResultSet(aqLayer)
>
> The result:
> marius@mobi:~/temp$ run_sql.py
> 1
> TELEORMAN
> marius@mobi:~/temp$
>
> This explains why none of my centroids were remotely close to the polygon.
>
> Is there something wrong with my query or should I file a bug?
>
> -marius
>
> On Mon, 2011-07-11 at 20:46 -0500, Marius Jigmond wrote:
>
> I suppose a piece of code speaks louder :):
>
> xsect = False
> for i in range(gridLayer.GetFeatureCount()):
> feat = gridLayer.GetFeature(i)
> geom = feat.GetGeometryRef()
> point = geom.Centroid()
> for j in range(aqLayer.GetFeatureCount()):
> aqfeat = aqLayer.GetFeature(j)
> aqgeom = aqfeat.GetGeometryRef()
> if point.Intersect(aqgeom):
> xsect = True
> print 'ibound = 1'
> break
> if xsect:
> feat.SetField('IBOUND', 1)
> gridLayer.SetFeature(feat)
> xsect = False
>
> -marius
>
> On Mon, 2011-07-11 at 20:34 -0500, Marius Jigmond wrote:
>
> Hi everyone,
>
> I am trying to test whether centroids of polygons lie/intersect within
> another polygon. I have tried Intersect, Within, and Contains but they
> always return false. Should these methods work for my intended purpose or do
> I need to implement a point in polygon function? Thanks.
>
> -marius
>
> _______________________________________________
> gdal-dev mailing list
> <A HREF="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</A>
> <A HREF="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</A>
>
> _______________________________________________
> gdal-dev mailing list
> <A HREF="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</A>
> <A HREF="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</A>
>
> _______________________________________________
> gdal-dev mailing list
> <A HREF="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</A>
> <A HREF="http://lists.osgeo.org/mailman/listinfo/gdal-dev">http://lists.osgeo.org/mailman/listinfo/gdal-dev</A>
>
</PRE>
</BLOCKQUOTE>
</BODY>
</HTML>