[gdal-dev] OGR Geometry methods

Luke Peterson luke.peterson at gmail.com
Tue Jul 12 12:31:03 EDT 2011


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 ogr
fileName = '/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/db788cb8/attachment-0001.html


More information about the gdal-dev mailing list