Re: [gdal-dev] OGR Geometry methods
Marius,
Please do file a ticket with a small shapefile that shows this error.
http://trac.osgeo.org/gdal/newticket
On Tue, Jul 12, 2011 at 8:30 AM, Marius Jigmond
mariusjigm...@hotmail.comwrote:
**
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?
-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
listgdal-dev@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing
listgdal-dev@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
--
Best regards,
Chaitanya kumar CH.
+91-9494447584
17.2416N 80.1426E
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR Geometry methods
On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond
mariusjigm...@hotmail.comwrote:
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.peter...@gmail.com
To: mariusjigm...@hotmail.com
-
Luke Peterson
- sent via mobile device -
On Jul 11, 2011 11:00 PM, Marius Jigmond mariusjigm...@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@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
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
RE: [gdal-dev] OGR Geometry methods
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.peter...@gmail.com
To: mariusjigm...@hotmail.com
CC: gdal-dev@lists.osgeo.org
On Tue, Jul 12, 2011 at 11:56 AM, Marius Jigmond mariusjigm...@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.peter...@gmail.com
To: mariusjigm...@hotmail.com
-
Luke Peterson
- sent via mobile device -
On Jul 11, 2011 11:00 PM, Marius Jigmond mariusjigm...@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@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
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR Geometry methods
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@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] OGR Geometry methods
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?
-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@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
