On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond <span dir="ltr"><<a href="mailto:mariusjigmond@hotmail.com">mariusjigmond@hotmail.com</a>></span> wrote:<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div><div dir="ltr">
Strangely enough, ExecuteSQL works (it returns an osgeo.ogr.layer object with one feature) with<div class="im"><br>
<br>
'select FID from judeteEPSG3844 where DENJUD = "CLUJ"'<br>
<br></div>
but later a GetField('FID') request returns "Illegal field requested in GetField()"<br></div></div></blockquote><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div><div dir="ltr"><br>
-marius<br>
<br></div></div></blockquote><div><br></div><div><span class="Apple-style-span" style="border-collapse: collapse; font-family: arial, sans-serif; font-size: 13px; "><div>Another way to do this is to use the AttributeFilter. Run this by itself and see if you get the expected result:</div>
<div><br></div><div>from osgeo import ogr</div><div>fileName = '/data/romania/judeteEPSG3844.shp'</div><div>shapeDS = ogr.Open(fileName)</div><div>shapeLayer = shapeDS.GetLayerByIndex(0)</div><div>shapeLayer.SetAttributeFilter('DENJUD="CLUJ"')</div>
<div>targetFeature = shapeLayer.GetNextFeature()</div><div>print targetFeature.GetField('DENJUD') # should print the DENJUD of the feature fitting the attribute filter.</div></span></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div><div dir="ltr"><div><div><div class="h5"><div><div><blockquote style="border-left:1px #ccc solid;padding-left:1ex"><div><div dir="ltr"><div><div>Date: Tue, 12 Jul 2011 09:13:00 -0400<br>Subject: Re: [gdal-dev] OGR Geometry methods<br>
From: <a href="mailto:luke.peterson@gmail.com" target="_blank">luke.peterson@gmail.com</a><br>
To: <a href="mailto:mariusjigmond@hotmail.com" target="_blank">mariusjigmond@hotmail.com</a></div><div><div></div><div><br><br><br>
-----<br>
Luke Peterson<br>
- sent via mobile device -<br>
On Jul 11, 2011 11:00 PM, "Marius Jigmond" <<a href="mailto:mariusjigmond@hotmail.com" target="_blank">mariusjigmond@hotmail.com</a>> wrote:<br>
><br>
> 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:<br>
><br>
> #!/usr/bin/python<br>
><br>
> from osgeo import ogr<br>
> import sys, math, time, os<br>
><br>
> aquifer = '/data/romania/judeteEPSG3844.shp'<br>
> aqDS = ogr.Open(aquifer)<br>
> sql = 'select * from judeteEPSG3844 where DENJUD = "CLUJ"'<br>
> aqLayer = aqDS.ExecuteSQL(sql)<br>
> feat = aqLayer.GetFeature(0)<br>
> print aqLayer.GetFeatureCount()<br>
> print feat.GetField('DENJUD')<br>
> aqDS.ReleaseResultSet(aqLayer)<br>
><br>
> The result:<br>
> marius@mobi:~/temp$ run_sql.py <br>
> 1<br>
> TELEORMAN<br>
> marius@mobi:~/temp$ <br>
><br>
> This explains why none of my centroids were remotely close to the polygon.<br>
><br>
> Is there something wrong with my query or should I file a bug?<br>
><br>
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):<br>
aqLayer = aqDS.GetLayerByIndex(0) #gets the whole shapefile<br>
aqFID = aqLayer.GetFIDColumn() # I forget the method name here... hope this is right but may have to pull in the layer def<br>
sql = 'select %s from judeteEPSG3844 where DENJUD = "CLUJ"' % aqFID # specifically asking for just the FID field in the resultset<br>
sqlResults = aqDS.ExecuteSQL(sql) # sql results same as before<br>
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)<br>
feat = aqLayer.GetFeature(featFID) # using the fid to grab the feature from the 0th layer in the shapefile<br>
This may be a long way around ... <br><br>
> -marius<br>
><br>
><br>
> On Mon, 2011-07-11 at 20:46 -0500, Marius Jigmond wrote:<br>
>><br>
>> I suppose a piece of code speaks louder :):<br>
>><br>
>> xsect = False<br>
>> for i in range(gridLayer.GetFeatureCount()):<br>
>> feat = gridLayer.GetFeature(i)<br>
>> geom = feat.GetGeometryRef()<br>
>> point = geom.Centroid()<br>
>> for j in range(aqLayer.GetFeatureCount()):<br>
>> aqfeat = aqLayer.GetFeature(j)<br>
>> aqgeom = aqfeat.GetGeometryRef()<br>
>> if point.Intersect(aqgeom):<br>
>> xsect = True<br>
>> print 'ibound = 1'<br>
>> break<br>
>> if xsect:<br>
>> feat.SetField('IBOUND', 1)<br>
>> gridLayer.SetFeature(feat)<br>
>> xsect = False<br>
>><br>
>> -marius<br>
>><br>
>> On Mon, 2011-07-11 at 20:34 -0500, Marius Jigmond wrote:<br>
>>><br>
>>> Hi everyone,<br>
>>><br>
>>> 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.<br>
>>><br>
>>> -marius<br>
>>><br>
>>> _______________________________________________<br>
>>> gdal-dev mailing list<br>
>>> <a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
>>> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
>><br>
>> _______________________________________________<br>
>> gdal-dev mailing list<br>
>> <a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
>> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
><br>
><br>
> _______________________________________________<br>
> gdal-dev mailing list<br>
> <a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
<br></div></div></div> </div></div>
</blockquote></div><br></div></div></div></div> </div></div>
</blockquote></div><br>