4

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.

11

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)),)}
 geosjson['coordinates'] = np.round(np.array(geosjon['coordinates']),2)
 print geosjson
 {'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(geosjson)
 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 geosjson-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)
 geosjson['coordinates']= set_precision(geosjson['coordinates'], 2)
 print geosjson
 {'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(geosjson)
 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))
6

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))
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]
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, privacy policy and cookie policy

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