Rounding all coordinates in shapely

Question

Following on from How can I add points to a LineString in shapely?, I'm hitting a problem where floating point representations mean objects aren't precisely where I expect them to be after they are moved using affinity.translate.

This would be fixed if there was a way of rounding all the coordinates to the nearest cm (or other insignificantly small value).

Does this exist?

For what it's worth, I'm doing all my intersection and union operations first before the translation and that seems to be working.

Answer 1

There are a few instances where @gene's answer does not work.

For example, the using the overprecise value -73.92391000000001

geojson = {'type': 'Polygon', 'coordinates': [[[-73.92391, 41.31064], [-73.92391, 41.31069], [-73.92388, 41.31069], [-73.92388, 41.31064], [-73.92391, 41.31064]]]}
print(str(shape(geojson)))
POLYGON ((-73.92391000000001 41.31064, -73.92391000000001 41.31069, -73.92388 41.31069, -73.92388 41.31064, -73.92391000000001 41.31064))

EDIT (MAY 2019): I recommend a roundtrip using shapely.wkt.dumps method instead, where rounding_precision=n is the n'th decimal point to which you want the coordinates rounded:

shapely.wkt.loads(shapely.wkt.dumps(geom, rounding_precision=n))

Answer 2

Shapely and GEOS cannot reduce precision (floating-point precision problem) but you can use other functions as numpy.round() or round(),with the GeoJSON format.

For polygons

from shapely.geometry import shape, mapping
import numpy as np
# one polygon
print poly.wkt
POLYGON ((101.2200527190607 -114.6493019170767, 146.6225142079163 -114.4488495484725, 185.0301801801801 -114.2792792792793, 184.6581081081081 -217.7153153153153, 14.99324324324321 -218.4594594594595, 16.4815315315315 -115.0234234234234, 101.2200527190607 -114.6493019170767))
# convert to GeoJSON format
geojson = mapping(poly)
print geojson
{'type': 'Polygon', 'coordinates': (((101.2200527190607, -114.6493019170767), (146.6225142079163, -114.4488495484725), (185.0301801801801, -114.2792792792793), (184.6581081081081, -217.7153153153153), (14.99324324324321, -218.4594594594595), (16.4815315315315, -115.0234234234234), (101.2200527190607, -114.6493019170767)),)}
 geojson['coordinates'] = np.round(np.array(geojson['coordinates']),2)
 print geojson
 {'type': 'Polygon', 'coordinates': array([[[ 101.22, -114.65],
    [ 146.62, -114.45],
    [ 185.03, -114.28],
    [ 184.66, -217.72],
    [  14.99, -218.46],
    [  16.48, -115.02],
    [ 101.22, -114.65]]])}
 print shape(geojson)
 POLYGON ((101.22 -114.65, 146.62 -114.45, 185.03 -114.28, 184.66 -217.72, 14.99 -218.46, 16.48 -115.02, 101.22 -114.65))

If you don't want to use Numpy, you can adapt the function def _set_precision(coords, precision) of the geojson-precision module

 def set_precision(coords, precision):
    result = []
    try:
        return round(coords, int(precision))
    except TypeError:
        for coord in coords:
            result.append(set_precision(coord, precision))
    return result
 geojson = mapping(poly)
 geojson['coordinates']= set_precision(geojson['coordinates'], 2)
 print geojson
 {'type': 'Polygon', 'coordinates': [[[101.22, -114.65], [146.62, -114.45], [185.03, -114.28], [184.66, -217.72], [14.99, -218.46], [16.48, -115.02], [101.22, -114.65]]]}
 print shape(geojson)
 POLYGON ((101.22 -114.65, 146.62 -114.45, 185.03 -114.28, 184.66 -217.72, 14.99 -218.46, 16.48 -115.02, 101.22 -114.65))

Answer 3

As of shapely 2.0 (see release notes), a precision model is now available: the set_precision function and the grid_size parameter for overlay functions.

However, it's not necessarily behaving as you might expect, see some unexpected results with contains here.

Answer 4

while there isn't a built in method or function in shapely you can use shapely.ops.transform to write a utility function to round coordinates for any geometry, including multipart geometries.

from shapley.ops import transform

def round_coordinates(geom, ndigits=2):
    
   def _round_coords(x, y, z=None):
      x = round(x, ndigits)
      y = round(y, ndigits)

      if z is not None:
          z = round(z, ndigits) # corrected from x to z

      return [c for c in (x, y, z) if c is not None]
   
   return transform(_round_coords, geom)

then running it like so

>>> _p = Polygon(((0.01, 1.234), (5.4321, 2.9348343943), (0.00000001, 9.23923), (0.01, 1.234)))
>>> round_coordinates(_p, 2).wkt
'POLYGON ((0.01 1.23, 5.43 2.93, 0 9.24, 0.01 1.23))'

Answer 5

In my case, the accepted answer did not work as my geojson was a multi polygon. Doing np.round(np.array(geosjon['coordinates']),2) throws a Python exception.

Instead, I wrote a short little recursive function that rounds all floats in a nested list. Meaning, that it also works on a list of list of list of list... for any depth.

def round_nested_list(l, precision):
    def round_element(e):
        if isinstance(e, list):
            return round_nested_list(e, precision)
        else:
            return round(e, precision)
    return [round_element(e) for e in l]

# test
l = [[[1.1123, 1.1123],[[1.1123, 1.1123],[[1.1123, 1.1123],[1.1123, 1.1123]]]], [[1.1123, 1.1123],[[1.1231, 1.1123],[[1.1231, 1.1231],[1.1231, 1.1231]]]]]

round_nested_list(l, 2)

# result is [[[1.11, 1.11], [[1.11, 1.11], [[1.11, 1.11], [1.11, 1.11]]]],
 [[1.11, 1.11], [[1.12, 1.11], [[1.12, 1.12], [1.12, 1.12]]]]]

I tried on a multipolygon with 20k+ points, and it took a few seconds. Good for me.

Edit : And if the coordinates are in tuples (at the last level), you can use :

def round_nested_list(l, precision):
    def round_element(e):
        if isinstance(e, list):
            return round_nested_list(e, precision)
        elif isinstance(e, tuple):
            return tuple(round(i, precision) for i in e)
        else:
            return round(e, precision)
    return [round_element(e) for e in l]

Answer 6

I don't see any way to reduce precision with Shapely even though GEOS, the base geometry engine, has precision reducer classes. You might try ogr2ogr. The geoJson driver has a "Coordinate Precision" option. You can export to geojson then create a shapefile from geojson. Just an idea.