12

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

0

6 Answers 6

17

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))
2
  • isn't there a typo in your alternate code: result.append(_set_precision(coord, precision)) should be result.append(set_precision(coord, precision)) ???
    – Riccardo
    Commented Jan 31, 2022 at 21:48
  • Yes, corrected, thanks
    – gene
    Commented Feb 1, 2022 at 9:17
17

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))
2
  • 1
    And to remove sequential duplicates created by rounding: shapely.wkt.loads(...).simplify(0).
    – Ian
    Commented Apr 30, 2020 at 15:31
  • I like the simplicity of this answer. However what about performances ?
    – Astariul
    Commented Mar 24, 2021 at 2:48
2

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))'
1
  • imo this is the best answer since it uses the built in shapely transform method instead of converting to some other format. But it throws a weird error for me complaining about a generator. I fixed it by dropping the one line for loop. see here gis.stackexchange.com/a/445283
    – Shawn
    Commented Nov 16, 2022 at 14:12
1

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]
1

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.

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.