2

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

0

4 Answers 4

3

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))
2
  • 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.
    – romepa
    Mar 21 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?
    – romepa
    Mar 22 at 8:17
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))
7
  • Are there ways of using this ST_SnapToGrid function in Python or would I need to call ogrinfo in a commandline?
    – romepa
    Mar 21 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
    – user30184
    Mar 21 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?
    – romepa
    Mar 22 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/….
    – user30184
    Mar 22 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> >
    – romepa
    Mar 22 at 9:21
1

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

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

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.