Marius,<br><br>Please do file a ticket with a small shapefile that shows this error.<br><a href="http://trac.osgeo.org/gdal/newticket">http://trac.osgeo.org/gdal/newticket</a><br><br><div class="gmail_quote">On Tue, Jul 12, 2011 at 8:30 AM, Marius Jigmond <span dir="ltr"><<a href="mailto:mariusjigmond@hotmail.com">mariusjigmond@hotmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><u></u>
<div>
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>
<a href="mailto:marius@mobi" target="_blank">marius@mobi</a>:~/temp$ run_sql.py <br>
1<br>
TELEORMAN<br>
<a href="mailto:marius@mobi" target="_blank">marius@mobi</a>:~/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><font color="#888888">
<br>
-marius</font><div><div></div><div class="h5"><br>
<br>
On Mon, 2011-07-11 at 20:46 -0500, Marius Jigmond wrote:<br>
<blockquote type="CITE">
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>
<blockquote type="CITE">
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
<pre>_______________________________________________
gdal-dev mailing list
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
<pre>_______________________________________________
gdal-dev mailing list
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
</blockquote>
</div></div></div>
<br>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">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></blockquote></div><br><br clear="all"><br>-- <br>Best regards,<br>Chaitanya kumar CH.<br>
<br>+91-9494447584<br>17.2416N 80.1426E<br>