3

I'm working with kml files using Python gdal.ogr module. When I extract feature's geometry using feature.GetGeometryRef(), I get:

POLYGON ((44.184069865521963 41.971859752507207 0,44.355248797385123 41.976261187945859 0, ...

However, I do not need this Z value - I need just a plain 2D WKT. After hours of searching I found FlattenTo2D() method, but feature.GetGeometryRef().FlattenTo2D() returns None

So how do I convert 3D WKT to 2D?

6

Using of FlattenTo2D() method is correct way, but the method does not returns new object (like you expected) - it changes current object. So you should use the method like this:

g = ogr.CreateGeometryFromWkt("POINT(1 1 1)")
g.flattenTo2D()
print(g.ExportToWkt())
>> 'POINT (1 1)'
3

WKT is only a text markup language for representing vector geometry, therefore you cannot convert 3D WKT to 2D (= text to text), you need to change the geometry

The problem with ogr is that it seems that all the geometries are 3D by defaut

1) Creation of a 2D point with ogr

point = ogr.Geometry(ogr.wkbPoint)
# 2D point
point.AddPoint(1198054.34, 648493.09)
print point.ExportToWkt()
POINT (1198054.34 648493.09 0)
point.GetGeometryType()==ogr.wkbPoint
False
point.GetGeometryType()==ogr.wkbPoint25D
True

2) But, with the creation of a 2D point with Shapely

from shapely.geometry import Point
pt = Point(1198054.34, 648493.09)
# conversion to an ogr geometry
point = ogr.CreateGeometryFromWkb(pt.wkb)
print(point.ExportToWkt())
POINT (1198054.34 648493.09)
point.GetGeometryType()==ogr.wkbPoint
True
point.GetGeometryType()==ogr.wkbPoint25D
False

3) Therefore it is possible to modify the geometry and an easy solution is to use the GeoJson format (feature.ExportToJson())

# The resulting WKT 
feature.ExportToWkt()
'Polygon ((250325.46051841106964275 142166.43077902274671942 0, 250775.29772301684715785 142019.67537893861299381 0, 250906.10144917882280424 141828.25529187230858952 0, 250555.1646228906174656 141550.69616562619921751 0, 250309.5088444888824597 141238.04335675123729743 0, 249933.04933992517180741 141777.20993532129796222 0, 250325.46051841106964275 142166.43077902274671942 0))'

json_geom = feature.ExportToJson()
print json_geom
'{ "type": "Polygon", "coordinates": [ [ [ 250325.460518411069643, 142166.430779022746719, 0.0 ], [ 250775.297723016847158, 142019.675378938612994, 0.0 ], [ 250906.101449178822804, 141828.25529187230859, 0.0 ], [ 250555.164622890617466, 141550.696165626199218, 0.0 ], [ 250309.50884448888246, 141238.043356751237297, 0.0 ], [ 249933.049339925171807, 141777.209935321297962, 0.0 ], [ 250325.460518411069643, 142166.430779022746719, 0.0 ] ] ] }'
# convert json format to a Python dictionary
import json
geom = json.loads(json_geom)
print geom['coordinates']
[[[250325.4605184111, 142166.4307790227, 100.0], [250775.2977230168, 142019.6753789386, 150.0], [250906.1014491788, 141828.2552918723, 200.0], [250555.1646228906, 141550.6961656262, 230.0], [250309.5088444889, 141238.0433567512, 340.0], [249933.0493399252, 141777.2099353213, 220.0], [250325.4605184111, 142166.4307790227, 100.0]]]
# slice the coordinates to eliminate 3D
new_coords =  [[i[:2] for i in geom['coordinates'][0]]]
print new_coords
[[[250325.4605184111, 142166.4307790227], [250775.2977230168, 142019.6753789386], [250906.1014491788, 141828.2552918723], [250555.1646228906, 141550.6961656262], [250309.5088444889, 141238.0433567512], [249933.0493399252, 141777.2099353213], [250325.4605184111, 142166.4307790227]]]
new = { "type": "Polygon",'coordinates': new_coords}
print new
{'type': 'Polygon', 'coordinates': [[[250325.4605184111, 142166.4307790227], [250775.2977230168, 142019.6753789386], [250906.1014491788, 141828.2552918723], [250555.1646228906, 141550.6961656262], [250309.5088444889, 141238.0433567512], [249933.0493399252, 141777.2099353213], [250325.4605184111, 142166.4307790227]]]}
# recreate the ogr geometry
poly = ogr.CreateGeometryFromJson(json.dumps(new))
poly.ExportToWkt() 
'POLYGON ((250325.460518411098747 142166.430779022688512,250775.29772301678895 142019.675378938612994,250906.1014491787937 141828.25529187230859,250555.164622890588362 141550.696165626199218,250309.508844488911564 141238.043356751208194,249933.049339925200911 141777.209935321297962,250325.460518411098747 142166.430779022688512))'

You can also use geomet to transform the correct GeoJson to WKT

from geomet import wkt
wkt.dumps(new, decimals=4)
'POLYGON ((250325.4605 142166.4308, 250775.2977 142019.6754, 250906.1014 141828.2553, 250555.1646 141550.6962, 250309.5088 141238.0434, 249933.0493 141777.2099, 250325.4605 142166.4308))'
1

Though quite an old post but I was having the same issue and this helped me:

ring = ogr.Geometry(ogr.wkbLinearRing)

ring.AddPoint(coordinates['xmin'], coordinates['ymin'])
ring.AddPoint(coordinates['xmin'], coordinates['ymax'])
ring.AddPoint(coordinates['xmax'], coordinates['ymax'])
ring.AddPoint(coordinates['xmax'], coordinates['ymin'])
ring.CloseRings()

geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(ring)

geom.FlattenTo2D()

Thanks to programcreek

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.