#2692 closed defect (fixed)
OGR PGeo driver's GetGeometryRef() returns NULL for points with Z coordinates
Reported by: | jjr8 | Owned by: | warmerdam |
---|---|---|---|
Priority: | normal | Milestone: | 1.6.0 |
Component: | OGR_SF | Version: | 1.5.2 |
Severity: | normal | Keywords: | pgeo type 9 pointz |
Cc: |
Description
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.
Attachments (1)
Change History (5)
by , 15 years ago
Attachment: | TestPGeo.zip added |
---|
follow-up: 2 comment:1 by , 15 years ago
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.
comment:2 by , 15 years ago
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.
comment:3 by , 15 years ago
Keywords: | pgeo type 9 pointz added |
---|---|
Resolution: | → fixed |
Status: | new → closed |
Patch applied in r15799
comment:4 by , 15 years ago
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 >>>
Shapefile and geodatabase used in this example