1

I found something unexpected when trying to create a simple polygon, following the "cookbook" example in:

https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html

from osgeo import ogr

ring = ogr.Geometry(ogr.wkbLinearRing)
print(ring.GetGeometryType())
ring.AddPoint(1179091.1646903288, 712782.8838459781)
ring.AddPoint(1161053.0218226474, 667456.2684348812)
ring.AddPoint(1214704.933941905, 641092.8288590391)
ring.AddPoint(1228580.428455506, 682719.3123998424)
ring.AddPoint(1218405.0658121984, 721108.1805541387)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
print(ring.GetGeometryType())

# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
print(poly)
print(poly.GetGeometryType())

This gives as output:

2
-2147483646
POLYGON ((1179091.16469033 712782.883845978 0,1161053.02182265 667456.268434881 0,1214704.9339419 641092.828859039 0,1228580.42845551 682719.312399842 0,1218405.0658122 721108.180554139 0,1179091.16469033 712782.883845978 0))
-2147483645

i.e. the geometryType changes from 2 to some weird negative integer after creation, and similarly for the resulting polygon.

This matters, because later on, I want to check geometryType and the polygon is obviously no longer of type ogr.wkbPolygon (3)

As a workaround I check for GetGeometryName() instead, but this is definitely odd.

Guido Lemoine

NB. Python 3.6.9 with the latest GDAL source build (3.1.0) but also different older versions.

6
1

OGR geometry type codes of 25D geometries are small negative numbers https://gist.github.com/walkermatt/7121427

-2147483647 Point25D
-2147483646 LineString25D
-2147483645 Polygon25D
-2147483644 MultiPoint25D
-2147483643 MultiLineString25D
-2147483642 MultiPolygon25D
          0 Geometry
          1 Point
          2 Line
          3 Polygon
          4 MultiPoint
          5 MultiLineString
          6 MultiPolygon
        100 No Geometry

Function AddPoint adds a point with x, y, and z coordinates https://gdal.org/python/osgeo.ogr-pysrc.html

 def AddPoint(self, *args, **kwargs): 
 """AddPoint(Geometry self, double x, double y, double z=0)""" 
 return _ogr.Geometry_AddPoint(self, *args, **kwargs) 

There is another function AddPoint_2D that is using only x and y.

def AddPoint_2D(self, *args): 
"""AddPoint_2D(Geometry self, double x, double y)""" 
return _ogr.Geometry_AddPoint_2D(self, *args) 

You have two alternatives for getting answer "2" in your first example.

  1. Build the geometry as you did but ask the geometry type as flattened into 2D

    print(ogr.GT_Flatten(ring.GetGeometryType()))
    2
    
  2. Build the geometry as 2D from the beginning

    from osgeo import ogr
    
    ring = ogr.Geometry(ogr.wkbLinearRing)
    print(ring.GetGeometryType())
    ring.AddPoint_2D(1179091.1646903288, 712782.8838459781)
    ring.AddPoint_2D(1161053.0218226474, 667456.2684348812)
    ring.AddPoint_2D(1214704.933941905, 641092.8288590391)
    ring.AddPoint_2D(1228580.428455506, 682719.3123998424)
    ring.AddPoint_2D(1218405.0658121984, 721108.1805541387)
    ring.AddPoint_2D(1179091.1646903288, 712782.8838459781)
    print(ring.GetGeometryType())
    2
    
1
  • Thanks! I did not know about the 25D geometries. Trouble is that GetGeometryName does not mention the 25D part, otherwise it would have given some clue about this. May 20 '20 at 7:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.