【问题标题】:Why is transforming a shapely polygon not working in some cases?为什么在某些情况下转换形状多边形不起作用?
【发布时间】:2020-11-12 00:29:25
【问题描述】:

我正在尝试使用shapely 计算地理坐标多边形的大小,这似乎需要转换为合适的投影才能产生以平方米为单位的结果。我在网上找到了几个示例,但我无法让它适用于我的示例多边形。

因此,我尝试使用我发现的代码 sn-ps 附带的相同示例多边形,我注意到它适用于某些整体而不适用于其他部分。为了重现结果,这是最小的示例代码:

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)

coords1 和 coords2 都只是从据说可以工作的代码 sn-ps 复制而来。但是,只有coords2 适合我。我已经使用https://geojsonlint.com 来查看两个多边形之间是否存在差异,并且多边形的惯用手/方向似乎不是有效的 GeoJSON。我不知道shapely 是否在乎,但颠倒顺序——然后https://geojsonlint.com 说它是有效的 GeoJSON,它会在地图上显示多边形——不会改变错误。

所以,它适用于 coords2,但是当我使用 coords1 时,我收到以下错误:

~/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)

我认为coords1(以及我自己数据中的示例多边形)有一些不同的地方会导致问题,但我不知道与coords2相比有什么不同。

简而言之,coords1 和coords2 有什么区别,一个工作一个不工作?

更新:我通过将always_xy=True 添加到预测的定义中使其工作。加上shapely 提供的新语法,避免partial,工作的sn-p 看起来像这样:

    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)

说实话,即使在阅读了文档之后,我也不知道always_xy 在做什么。因此,我不想提供答案。

【问题讨论】:

    标签: python gis polygon geojson shapely


    【解决方案1】:

    我认为你做得很好,只是反转不会创建新数据集。

    尝试使用此功能创建倒序列表:

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

    像这样执行函数:

    coords = rev_slice(coords1)
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2021-11-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2021-01-18
      • 2021-09-29
      • 2021-11-20
      相关资源
      最近更新 更多