GetPointCount() does not work

[dpinte@…]

Hi,

I've tested the following under Debian/Sid with python-gdal and gdal 1.3.1 and
under Windows XP with FWTools 1.0.a7 (gdal 1.3.1).

Using a shapefile containing 11 complicated polylines, the GetPointCount() on
myfeature.GetGeometryRef() does not work. It always give 0 ...

Here is my test script :

import os
import sys
sys.path.append('.')

import ogr

def test(network_file=None):
  if network_file == None:
	  network_file = 'D:\\ALERT\\data\\gis\\polygones\\managepoly_and.shp'
  
  lnet = ogr.Open(network_file)
 
  mylayer = lnet.GetLayer()
  for i in xrange(mylayer.GetFeatureCount()) :
    feature = mylayer.GetFeature(i)
    geom = feature.GetGeometryRef()
    print geom.GetPointCount()
 
if __name__ == '__main__':
  test()


I can provide my test file if needed.

Didrik

[warmerdam]

Yes, I will need the file. 

If it is big to attach here (say larger than 500KB), then
just stage it somewhere and I'll download it (or email it
directly to me if less than 10MB). 

[dpinte@…]

the testcase data file

[dpinte@…]

the file is uploaded. I've still not found any workaround.

Didrik

[warmerdam]

Didrik,

POLYGON objects don't own points directly.  Instead they own zero or more
rings.  The rings have points.  So you need to change your script around
a bit.  Also, polygon shapefiles can have MULTIPOLYGONs in some cases, 
so you also need to be ready for that.   An adjust form of the script might
look like this:

import os
import sys
sys.path.append('.')

import ogr

def report_geometry( geom, prefix = '  ' ):
  if geom.GetPointCount() > 0:
    print '%s%s (%d points)' % (prefix,geom.GetGeometryName(),geom.GetPointCount())
  else:
    print '%s%s' % (prefix,geom.GetGeometryName())
    
  if geom.GetGeometryCount() > 0:
    for i in xrange(geom.GetGeometryCount()):
      report_geometry( geom.GetGeometryRef(i), prefix + '  ' )

  
    
def test(network_file=None):
  if network_file == None:
	  network_file = 'managepoly_and.shp'
  
  lnet = ogr.Open(network_file)
 
  mylayer = lnet.GetLayer()
  for i in xrange(mylayer.GetFeatureCount()) :
    feature = mylayer.GetFeature(i)
    print 'Feature %d' % i
    geom = feature.GetGeometryRef()
    if geom is not None:
      report_geometry( geom )
 
if __name__ == '__main__':
  test()


This would produce:
Feature 0
  POLYGON
    LINEARRING (1851 points)
Feature 1
  POLYGON
    LINEARRING (1188 points)
Feature 2
  POLYGON
    LINEARRING (1264 points)
Feature 3
  POLYGON
    LINEARRING (779 points)
Feature 4
  POLYGON
    LINEARRING (816 points)
Feature 5
  POLYGON
    LINEARRING (794 points)
Feature 6
  POLYGON
    LINEARRING (560 points)
Feature 7
  POLYGON
    LINEARRING (616 points)
Feature 8
  POLYGON
    LINEARRING (268 points)
Feature 9
  POLYGON
    LINEARRING (1039 points)
Feature 10
  POLYGON
    LINEARRING (2608 points)