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

Question

(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.

Answer 1

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))

Answer 2

How about using the ST_SnapToGrid function that comes from SpatiaLite https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

return a new Geometry corresponding to the input Geometry; all points and vertices will be snapped to the grid defined by its origin and size(s).
Removes all consecutive points falling on the same cell.
All collapsed geometries will be stripped from the returned Geometry.
NULL is returned on failure.

Usage example with ogrinfo follows. With Python you would use the ExecuteSQL function. Don't pay attention to py.shp, it is in the command only because ogrinfo requires some datasource or it gives an error.

ogrinfo -dialect sqlite -sql "SELECT ST_SnapToGrid(ST_GeomFromText('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))'),1) AS geom" py.shp
INFO: Open of 'py.shp'
      using driver 'ESRI Shapefile' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Extent: (75199.000000, 451968.000000) - (75695.000000, 452398.000000)
Layer SRS WKT:
(unknown)
Geometry Column = geom
OGRFeature(SELECT):0
  POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 452398,75571 452398,75695 452183))

Answer 3

Use the following script (includes the application of json package)

import json
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)

coordinates = json.loads(geom.ExportToJson())['coordinates'][0]
coordinates_rounded = [[round(pair[0]), round(pair[1])] for pair in coordinates]

# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)

for pair in coordinates_rounded:
    ring.AddPoint_2D(pair[0], pair[1])
   
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

print(poly.ExportToWkt())

to achieve this:

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

Answer 4

Use simply the Regular expression module as in Find, round and replace numbers in a string with regexp for example

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))'
import re
# 0 decimals
re.sub(r'\d*\.\d+', lambda m: format(float(m.group(0)), '.0f'), testpolywkt)
'POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 452398,75571 452398,75695 452183))'
# 2 decimals
re.sub(r'\d*\.\d+', lambda m: format(float(m.group(0)), '.2f'), testpolywkt)
'POLYGON ((75694.86 452182.81,75570.78 451967.90,75322.61 451967.90,75198.53 452182.81,75322.61 452397.73,75570.78 452397.73,75694.86 452182.81))'
from osgeo import ogr
wkt= re.sub(r'\d*\.\d+', lambda m: format(float(m.group(0)), '.0f'), testpolywkt)
poly = ogr.CreateGeometryFromWkt(wkt)
poly.ExportToWkt()
'POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 452398,75571 452398,75695 452183))'