Why is transforming a shapely polygon not working in some cases?

Viewed 2247

I'm trying to calculate the size of a polygon of geographic coordinates using shapely, which seems to require a transformation into a suitable projection to yield a results in square meter. I found a couple of examples online, but I couldn't get it working for my example polygon.

I therefore tried to use the same example polygons that came with the code snippets I found, and I noticed that it works for some whole not for others. To reproduce the results, here's the minimal example code:

import json
import pyproj

from shapely.ops import transform
from shapely.geometry import Polygon, mapping
from functools import partial

coords1 = [(-97.59238135821987, 43.47456565304017),
           (-97.59244690469288, 43.47962399877412),
           (-97.59191951546768, 43.47962728271748),
           (-97.59185396090983, 43.47456565304017),
           (-97.59238135821987, 43.47456565304017)]

coords1 = reversed(coords1) # Not sure if important, but https://geojsonlint.com says it's wrong handedness
                            # Doesn't seem to affect the error message though

coords2 = [(13.65374516425911, 52.38533382814119),
           (13.65239769133293, 52.38675829106993),
           (13.64970274383571, 52.38675829106993),
           (13.64835527090953, 52.38533382814119),
           (13.64970274383571, 52.38390931824483),
           (13.65239769133293, 52.38390931824483),
           (13.65374516425911, 52.38533382814119)]

coords = coords1 # DOES NOT WORK
#coords = coords2 # WORKS

polygon = Polygon(coords)

# Print GeoJON to check on https://geojsonlint.com
print(json.dumps(mapping(polygon)))

projection = partial(pyproj.transform, 
                     pyproj.Proj('epsg:4326'), 
                     pyproj.Proj('esri:54009'))

transform(projection, polygon)

Both coords1 and coords2 are just copied from code snippets that supposedly work. However, only coords2 works for me. I've used https://geojsonlint.com to see if there's a difference between the two polygons, and it seems that the handedness/orientation of the polygon is not valid GeoJSON. I don't know if shapely even cares, but reversing the order -- and https://geojsonlint.com says it's valid GeoJSON then, and it shows the polygon on the map -- does not change the error.

So, it works with coords2, but when I use coords1 I get the following error:

~/env/anaconda3/envs/py36/lib/python3.6/site-packages/shapely/geometry/base.py in _repr_svg_(self)
    398             if xmin == xmax and ymin == ymax:
    399                 # This is a point; buffer using an arbitrary size
--> 400                 xmin, ymin, xmax, ymax = self.buffer(1).bounds
    401             else:
    402                 # Expand bounds by a fraction of the data ranges

ValueError: not enough values to unpack (expected 4, got 0)

I assume there's something different about coords1 (and the example polygon from my own data) that causes the problem, but I cannot tell what could be different compared to coords2.

In short, what's the difference between coords1 and coords2, with one working and the other not?

UPDATE: I got it working by adding always_xy=True to the definition of the projections. Together with the newer syntax provided by shapely, avoiding partial, the working snippet looks like this:

    project = pyproj.Transformer.from_proj(
        pyproj.Proj('epsg:4326'), # source coordinate system
        pyproj.Proj('epsg:3857'),
        always_xy=True
    ) # destination coordinate system
 
    transform(project.transform, polygon)

To be honest, even after reading the docs, I don't really know what always_xy is doing. Hence I don't want to provide is an answer.

1 Answers

i think you did good, only that the reversed does not create new dataset.

try to use this function to create reversed order list:

def rev_slice(mylist):
    '''
    return a revered list
    mylist: is a list
    '''
    a = mylist[::-1]
    return a

execute the function like so:

coords = rev_slice(coords1)
Related