[gdal-dev] OGR Geometry methods

Tue Jul 12 12:45:13 EDT 2011

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: luke.peterson at gmail.com
To: mariusjigmond at hotmail.com
CC: gdal-dev at lists.osgeo.org

On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond <mariusjigmond at hotmail.com> 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: luke.peterson at gmail.com

To: mariusjigmond at hotmail.com

-----

Luke Peterson

- sent via mobile device -

On Jul 11, 2011 11: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?

>

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

>>> 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/20110712/eb9eb126/attachment.html