[gdal-dev] OGR Geometry methods

[Marius Jigmond]

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.

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):
1. Select polygon(s) based on a SQL filter
2. Iterate over polygons from 1. and add geometries to a list
3. Iterate over a mesh of polygons, grab their centroids and verify
whether they are within the geometries from 2.

Script that fails (assumes I'm getting more than one geometry back and
adds it to a list):
#!/usr/bin/python

from osgeo import ogr
import sys, math, time, os

grid = 'model1.shp'
gridDS = ogr.Open(grid, 1)
gridLayer = gridDS.GetLayer(0)

aquifer = '/data/romania/judeteEPSG3844.shp'
aqDS = ogr.Open(aquifer)
aqLayer = aqDS.GetLayerByIndex(0)
aqLayer.SetAttributeFilter('DENJUD = "CLUJ"')
aqfeat = aqLayer.GetNextFeature()
aqgeom = aqfeat.GetGeometryRef()
geomList = [aqgeom]

while aqfeat is not None: #loop thru features in case of
islands/multiple polygons
  aqfeat = aqLayer.GetNextFeature()
  if aqfeat is not None:
    geomList.append(aqfeat.GetGeometryRef())

for i in range(gridLayer.GetFeatureCount()):
  feat = gridLayer.GetFeature(i)
  print feat.GetField('CELL_ID') #see how many features processed before
segfault, usually first two features
  geom = feat.GetGeometryRef()
  point = geom.Centroid()
  for polygeom in geomList:
    if polygeom.Contains(point):
      feat.SetField('IBOUND', 1)
      gridLayer.SetFeature(feat)
      break #break out of "for" if any poly contains the centroid

Script that works (assumes I'm only getting one geometry back, no need
for list):
#!/usr/bin/python

from osgeo import ogr
import sys, math, time, os

grid = 'model1.shp'
gridDS = ogr.Open(grid, 1)
gridLayer = gridDS.GetLayer(0)

aquifer = '/data/romania/judeteEPSG3844.shp'
aqDS = ogr.Open(aquifer)
aqLayer = aqDS.GetLayerByIndex(0)
aqLayer.SetAttributeFilter('DENJUD = "CLUJ"')
aqfeat = aqLayer.GetNextFeature()
aqgeom = aqfeat.GetGeometryRef()

for i in range(gridLayer.GetFeatureCount()):
  feat = gridLayer.GetFeature(i)
  geom = feat.GetGeometryRef()
  point = geom.Centroid()
  if aqgeom.Contains(point):
    feat.SetField('IBOUND', 1)
    gridLayer.SetFeature(feat)


-marius

On Mon, 2011-07-11 at 22:01 -0700, Frank Warmerdam wrote:

> Marius,
> 
> I think a bug is called for.
> 
> On Mon, Jul 11, 2011 at 8:00 PM, Marius Jigmond
> <mariusjigmond at hotmail.com> 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 at mobi:~/temp$ run_sql.py
> > 1
> > TELEORMAN
> > marius at 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
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> >
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> >
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
> >
> 
> 
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20110713/8503a378/attachment-0001.html