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

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

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)), round(ring.GetPoint(j)))
#also works as coordinates = (round(ring.GetX(j)), round(ring.GetY(j)))

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

print (poly.ExportToWkt())
``````

to get this:

``````POLYGON ((75695 452183,75571 451968,75323 451968,75199 452183,75323 452398,75571 452398,75695 452183))
``````
• Thanks @Taras, this is somewhat similar as to what I have tried. However, could you explain what is the benefit of using these rings? One of the benefits is that your solution is all `ogr` and mine uses `shapely`. Do you have any opinion on what would be more efficient? I want to apply it to large datasets, so every percent could potentially matter. Mar 21, 2022 at 16:39
• I have now tried and compared this and it seems to be a tiny bit faster than my solution: 31.5s vs 32.5s on my test dataset. Every little helps! Maybe because now I can do without `shapely`? Mar 22, 2022 at 8:17

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))
``````
• Are there ways of using this ST_SnapToGrid function in Python or would I need to call ogrinfo in a commandline? Mar 21, 2022 at 16:45
• I wrote `With Python you would use the ExecuteSQL function`. It is documented in gdal.org/python/osgeo.ogr.DataSource-class.html Mar 21, 2022 at 16:53
• I´m not sure how to use this one. When I use the same SELECT query (querystring) as in your example and use the ExecuteSQL method on the datasource (`driver_hex = ogr.GetDriverByName('ESRI Shapefile') datasource = driver_hex.Open('shapefile.shp', 0) datasource.ExecuteSQL(querystring)`), then I get: `ERROR 1: Undefined function 'ST_GeomFromText' used.` What am I doing wrong? Mar 22, 2022 at 9:09
• Did you set the SQL dialect SQLite? You can find usage examples from GDAL autotests github.com/OSGeo/gdal/blob/master/autotest/ogr/…. Mar 22, 2022 at 9:13
• Ok, changed to `result = ds_hex.ExecuteSQL(querystring, dialect = "SQLITE")` It returns a Layer object? `<osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x7f50f1987870> >` Mar 22, 2022 at 9:21

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_rounded = [[round(pair), round(pair)] for pair in coordinates]

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

for pair in coordinates_rounded:

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

print(poly.ExportToWkt())
``````

to achieve this:

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

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