Indeed, SetAttributeFilter returns the right feature. Thanks Luke.
FYI the ExecuteSQL bug submission is http://trac.osgeo.org/gdal/ticket/4156 -marius Date: Tue, 12 Jul 2011 12:31:03 -0400 Subject: Re: [gdal-dev] OGR Geometry methods From: [email protected] To: [email protected] CC: [email protected] On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond <[email protected]> wrote: Strangely enough, ExecuteSQL works (it returns an osgeo.ogr.layer object with one feature) with 'select FID from judeteEPSG3844 where DENJUD = "CLUJ"' but later a GetField('FID') request returns "Illegal field requested in GetField()" -marius Another way to do this is to use the AttributeFilter. Run this by itself and see if you get the expected result: from osgeo import ogrfileName = '/data/romania/judeteEPSG3844.shp'shapeDS = ogr.Open(fileName)shapeLayer = shapeDS.GetLayerByIndex(0)shapeLayer.SetAttributeFilter('DENJUD="CLUJ"') targetFeature = shapeLayer.GetNextFeature()print targetFeature.GetField('DENJUD') # should print the DENJUD of the feature fitting the attribute filter. Date: Tue, 12 Jul 2011 09:13:00 -0400 Subject: Re: [gdal-dev] OGR Geometry methods From: [email protected] To: [email protected] ----- Luke Peterson - sent via mobile device - On Jul 11, 2011 11:00 PM, "Marius Jigmond" <[email protected]> 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? > Can you select the feature by its FID instead? (apologies, I'm banging this out on my phone ... probably won't run but you should get the idea): aqLayer = aqDS.GetLayerByIndex(0) #gets the whole shapefile aqFID = aqLayer.GetFIDColumn() # I forget the method name here... hope this is right but may have to pull in the layer def sql = 'select %s from judeteEPSG3844 where DENJUD = "CLUJ"' % aqFID # specifically asking for just the FID field in the resultset sqlResults = aqDS.ExecuteSQL(sql) # sql results same as before featFID = sqlResults.GetFeature(0).GetFieldAsInt(aqFID) #this should be the actual FID of the results (problems would occur if your sql returned more than one hit) feat = aqLayer.GetFeature(featFID) # using the fid to grab the feature from the 0th layer in the shapefile This may be a long way around ... > -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 >>> [email protected] >>> http://lists.osgeo.org/mailman/listinfo/gdal-dev >> >> _______________________________________________ >> gdal-dev mailing list >> [email protected] >> http://lists.osgeo.org/mailman/listinfo/gdal-dev > > > _______________________________________________ > gdal-dev mailing list > [email protected] > http://lists.osgeo.org/mailman/listinfo/gdal-dev
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
