Shapely 2.0 doesn't accept Decimal Fields in the Polygon constructor.

[ryanneilyoung]

Expected behavior and actual behavior.

Before the Shapely 2.0 update our code worked:

    def to_geo(self):
        p1 = (self.Longitude1, self.Latitude1)
        p2 = (self.Longitude2, self.Latitude2)
        p3 = (self.Longitude3, self.Latitude3)
        p4 = (self.Longitude4, self.Latitude4)
        p5 = (self.Longitude1, self.Latitude1)
        poly_coords = (p1, p2, p3, p4, p5)

        poly = Polygon(poly_coords)
        poly = poly.convex_hull

        return poly

now it returns a:
ValueError: Inconsistent coordinate dimensionality

I've run this also with the latest on master and gotten the same result. We are able to mitigate this by casting the Decimals to floats, but I'd rather not if this is truly a regression and not an intended change.

now it

Steps to reproduce the problem.

 # the following are latitude/longitude values normally on a Django model, but I've just changed them to the value for reproducibility. 
 from shapely.geometry import Polygon
 from decimal import Decimal
 p1 = (Decimal('-76.572088'), Decimal('39.256796'))
 p2 = (Decimal('-76.569245'), Decimal('39.258841'))
 p3 = (Decimal('-76.567123'), Decimal('39.257046'))
 p4 = (Decimal('-76.571833'), Decimal('39.255359'))
 p5 = (Decimal('-76.572088'), Decimal('39.256796'))
 poly_coords = (p1, p2, p3, p4, p5)
 
 poly = Polygon(poly_coords)
 poly = poly.convex_hull

Operating system

Ubuntu 22.04.1 LTS, but also recreated in a python docker image

Shapely version and provenance

shapely==2.0.0 installed from PyPI using pip

[idanmiara]

The root of the issue that is that NumPy doesn't recognize decimal.Decimal as a specific type. The closest it can get is the most general dtype, object. (https://stackoverflow.com/a/7772386)
Thus the following line:
https://github.com/shapely/shapely/blob/main/shapely/geometry/polygon.py#L93
creates an object numpy array which the Polygon constructor doesn't handle.

Not sure if this is a known limitation.

[ryanneilyoung]

I'm very comfortable with the answer being "we've dropped Decimal support for Polygons, you'll have to cast", but I do want to make sure that it is a known regression from this new release and should either be fixed or documented somewhere.

[caspervdw]

I think it is a low effort to keep supporting decimal.Decimal, for the single geometry constructors. For the vectorized versions (points, linestrings, linearrings) we can't do that as numpy doesn't support it.

So I am 👍 of fixing the regression.

[idanmiara]

I think it is a low effort to keep supporting decimal.Decimal, for the single geometry constructors. For the vectorized versions (points, linestrings, linearrings) we can't do that as numpy doesn't support it.

So I am 👍 of fixing the regression.

I could make a PR to fix this, but I need someone to review my PRs...
Would you be able to review my opened PRs please?

[jorisvandenbossche]

Yeah, it seems that ctypes (or cython, for functions that had a speedups version) automatically converted Decimals to floats where needed. Also +1 on keeping a general conversion to float for individual coordinates in the class constructors.

We currently have this function in several places (for linestrings/linearring):

shapely/shapely/geometry/polygon.py

Lines 86 to 91 in 5a3d023

	# check coordinates on points
	def _coords(o):
	if isinstance(o, Point):
	return o.coords[0]
	else:
	return o

We can probably update the else: return o to else: return [float(c) for c in o]

Would you be able to review my opened PRs please?

I am prioritizing issues/PRs for 2.0.1, so yes you can be ensured to get a review.

[idanmiara]

I am prioritizing issues/PRs for 2.0.1, so yes you can be ensured to get a review.

OK, I'll be happy to fix this 👍
Just wondering if you are going to merge #1681 and #1682 for 2.0.1 as it would be more convenient for me to fix this issue (and add the required tests) once those are merged.

[dotysan]

I landed here because my app reads billions of GeoJSON items and converts the feature items to shapes for geoprocessing. I use ijson for memory efficiency. And after upgrading shapely from 1.8 to 2.0, I started getting the ValueError: Inconsistent coordinate dimensionality exception.

It appears that ijson protects floats from floating point drift with Decimal() by default. The trick is to convert this:

    with open(geojson) as f:
        for item in ijson.items(f, 'features.item'):
            geom = item['geometry']
            shape = geometry.shape(geom)

Into this:

    with open(geojson) as f:
        for item in ijson.items(f, 'features.item', use_float=True):
            geom = item['geometry']
            shape = geometry.shape(geom)

Since my GeoJSON is lat/lon decimal degrees with at least 8-digits of precision, we don't care about eliminating the floating point drift as its buried well below the perceivable accuracy on our maps.

And now use_float=True apparently runs even faster without trying to wrap every value in Decimal()!

Just passing this info on as an FYI...

[jorisvandenbossche]

@dotysan thanks for the heads up about that! Apart from speed, it also makes sense to directly keep floats (and not convert to Decimal), since shapely (GEOS) would otherwise convert the coordinates to floating point values anyway.