OGR PGeo driver's GetGeometryRef() returns NULL for points with Z coordinates

[jjr8]

The OGR shapefile driver cam read shapefiles that have Z coordinates. It recognizes the geometry of these features as wkbPoint25D and returns all three coordinates. Here is a Python script that demonstrates this.

First call ArcGIS APIs to create a point shapefile with Z coordinates and add one feature:

>>> import os
>>> os.mkdir('C:\\TestPGeo')
>>> from win32com.client import Dispatch
>>> gp = Dispatch('esriGeoprocessing.GPDispatch')
>>> gp.CreateFeatureClass_management('C:\\TestPGeo', 'TestPoints', 'POINT', None, 'DISABLED', 'ENABLED')    # Create point shapefile with Z coords
u'C:\\TestPGeo\\TestPoints.shp'
>>> cur = gp.InsertCursor('C:\\TestPGeo\\TestPoints.shp')
>>> row = cur.NewRow(); point = gp.CreateObject('Point'); point.X = 1; point.Y = 2; point.Z = 3; row.SetValue('Shape', point); cur.InsertRow(row)
>>> del cur

Now call ogr Python bindings to read the feature's geometry:

>>> from osgeo import ogr
>>> ds = ogr.Open('C:\\TestPGeo\\TestPoints.shp')
>>> layer = ds.GetLayer(0)
>>> feature = layer.GetNextFeature()
>>> geometry = feature.GetGeometryRef()
>>> geometry.GetGeometryType() == ogr.wkbPoint25D
True
>>> geometry.GetX(); geometry.GetY(); geometry.GetZ()
1.0
2.0
3.0

PGeo driver cannot do this. If the feature is a point with Z coords, GetGeometryRef() returns None from the Python bindings. Here's the demo:

First create a personal GDB with a point feature class having Z coords and add one feature:

>>> gp.CreatePersonalGDB_management('C:\\TestPGeo', 'MyGeoDB')
u'C:\\TestPGeo\\MyGeoDB.mdb'
>>> gp.CreateFeatureClass_management('C:\\TestPGeo\\MyGeoDB.mdb', 'TestPoints', 'POINT', None, 'DISABLED', 'ENABLED')    # Create point feature class with Z coords
u'C:\\TestPGeo\\MyGeoDB.mdb\\TestPoints'
>>> cur = gp.InsertCursor('C:\\TestPGeo\\MyGeoDB.mdb\\TestPoints')
>>> row = cur.NewRow(); point = gp.CreateObject('Point'); point.X = 4; point.Y = 5; point.Z = 6; row.SetValue('Shape', point); cur.InsertRow(row)
>>> del cur

Now try to read the feature's geometry:

>>> ds = ogr.Open('C:\\TestPGeo\\MyGeoDB.mdb')
>>> layer = ds.GetLayerByName('TestPoints')
>>> feature = layer.GetNextFeature()
>>> geometry = feature.GetGeometryRef()
>>> geometry.GetGeometryType() == ogr.wkbPoint25D
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
AttributeError: 'NoneType' object has no attribute 'GetGeometryType'
>>> geometry is None
True

I tested this by calling GDAL 1.5 from Python 2.5.2 as described in #2688. It was a Windows machine with ArcGIS 9.1 SP2. The version of Arc is only relevant if you care about how the personal geodb was created. I have not yet tried it with newer versions of ArcGIS, although I doubt they changed the binary format of point Z geometries in the newer 9.x versions.

I am attaching the shapefile and personal GDB created in this experiment.

For my application, it is not critical that the PGeo driver be able to read points with Z coords. I discovered this in testing (I have my own wrapper around OGR inside my app). For completeness and consistency inside OGR, it would be great if PGeo could read points with Zs. But my app can call the much slower ArcGIS cursor APIs if needed.

[jjr8]

Shapefile and geodatabase used in this example

[Even Rouault]

You can try the following patch that maps type 9 on SHPT_POINTZ. But I can't find any reference where 9 is a legitimate alias for SHPT_POINTZ...

Index: ogr/ogrsf_frmts/pgeo/ogrpgeolayer.cpp
===================================================================
--- ogr/ogrsf_frmts/pgeo/ogrpgeolayer.cpp	(révision 15790)
+++ ogr/ogrsf_frmts/pgeo/ogrpgeolayer.cpp	(copie de travail)
@@ -445,6 +445,8 @@
 /* -------------------------------------------------------------------- */
     if( nSHPType == 50 )
         nSHPType = SHPT_ARC;
+    else if (nSHPType == 9 )
+        nSHPType = SHPT_POINTZ;
 
 /* ==================================================================== */
 /*  Extract vertices for a Polygon or Arc.				*/

With that patch, I get:

Feature Count: 1
Extent: (4.000000, 5.000000) - (4.000000, 5.000000)
Layer SRS WKT:
(unknown)
OBJECTID: Integer (0.0)
OGRFeature(TestPoints):1
  OBJECTID (Integer) = 1
  POINT (4 5 6)

Is POINT(4 5 6) the expected value ? because the shapefile you've attached in the zip contains POINT(1 2 3) instead.

[jjr8]

Replying to rouault:

Is POINT(4 5 6) the expected value ? because the shapefile you've attached in the zip contains POINT(1 2 3) instead.

Yes, 4,5,6 are the correct coordinates in the attached geodatabase, while the shapefile contains 1,2,3. This is shown in the Python code above. Sorry for the confusion. It would have been clearer if I used the same coordinates in both.

[Even Rouault]

Patch applied in r15799

[jjr8]

I can confirm that this appears to be fixed in GDAL 1.6.0 with Python 2.5 on win32:

Python 2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> from osgeo import ogr
>>> ds = ogr.Open('C:\\TestPGeo\\MyGeoDB.mdb')
>>> layer = ds.GetLayerByName('TestPoints')
>>> feature = layer.GetNextFeature()
>>> geometry = feature.GetGeometryRef()
>>> geometry.GetGeometryType() == ogr.wkbPoint25D
True
>>> geometry.GetX(); geometry.GetY(); geometry.GetZ()
4.0
5.0
6.0
>>>