Converting OGR geometry as WKT of polygon from float to integer in the WKT string

coordinatesogrpolygonpythonwell-known-text

(This is a repost into this community of my original question posted on SO, because there has been no answer so far. Hopefully in this community more luck. I'll ask for the other one to be closed.)

I am working with a hexagonal grid that needs to be exported to text and each feature needs to be represented by a WKT-string so it can be later imported in a next step.

I use OGR in Python (3.9.7).

Problem: The geometry data for the points that make up the polygon in the WKT string are all as floats. I want them to be as integers instead, because that is accurate enough for my purposes as the decimal places represent micrometers in the coordinate system used (SRID=28992). When all points in the polygons will be rounded the same way, my grid will remain a fully covering grid.

Example:

from osgeo import ogr

testpolywkt='POLYGON ((75694.8564184426 452182.812141684,75570.7757705624 451967.898155319,75322.6144748018 451967.898155319,75198.5338269215 452182.812141684,75322.6144748018 452397.726128049,75570.7757705624 452397.726128049,75694.8564184426 452182.812141684))'

geom = ogr.CreateGeometryFromWkt(testpolywkt)

I want to get a WKT string like:

'POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 52398,75571 452398,75695 452183))'

What could be the best way to do that?

–Update:
Let me add what I have tried myself, which is to go point by point in the geometry and then round those numbers, like:

from shapely import wkt

for feat_hex in layer_hex:
    strShape = str(feat_hex.GetField('receptor_i')) + '\tSRID=28992;POLYGON(('

    geom_hex = feat_hex.GetGeometryRef()
    strWkt = geom_hex.ExportToWkt()

    # Get points from polygon
    p1 = wkt.loads(strWkt)
    x, y = p1.exterior.coords.xy
    
    # do this for each point
    for i in range(0,len(x)):
        if i > 0:
            strShape += ' '
        
        strShape += str(round(x[i])) + ' ' + str(round(y[i]))
               
        if i < len(x)-1:
            strShape += ',' # no comma at end
            
    strShape += '))\n'
    
    # write line to file 
    outputfile.write(strShape) # write line

It works, though this feels like a workaround rather than a more elegant solution. The question is whether there are some settings in the ogr package that somehow instruct the underlying C-library to round it or set a lower precision.

Best Answer

Otherwise try this script:

from osgeo import ogr

testpolywkt = 'POLYGON ((75694.8564184426 452182.812141684,75570.7757705624 451967.898155319,75322.6144748018 451967.898155319,75198.5338269215 452182.812141684,75322.6144748018 452397.726128049,75570.7757705624 452397.726128049,75694.8564184426 452182.812141684))'

geom = ogr.CreateGeometryFromWkt(testpolywkt)

numRings = geom.GetGeometryCount()
# Create ring
newRing = ogr.Geometry(ogr.wkbLinearRing)

for i in range(numRings):
    ring = geom.GetGeometryRef(i)
    numPoints = ring.GetPointCount()
    for j in range(numPoints):
        coordinates = (round(ring.GetPoint(j)[0]), round(ring.GetPoint(j)[1]))
        #also works as coordinates = (round(ring.GetX(j)), round(ring.GetY(j)))
        newRing.AddPoint_2D(coordinates[0], coordinates[1])

# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(newRing)

print (poly.ExportToWkt())

to get this:

POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 452398,75571 452398,75695 452183))
Related Question