Ticket #2692 (closed defect: fixed)
OGR PGeo driver's GetGeometryRef() returns NULL for points with Z coordinates
|Reported by:||jjr8||Owned by:||warmerdam|
|Severity:||normal||Keywords:||pgeo type 9 pointz|
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.