我正在尝试使用 shapely
计算地理坐标多边形的大小,这似乎需要转换为合适的投影才能产生以平方米为单位的结果。我在网上找到了几个示例,但我无法让它适用于我的示例多边形。
因此,我尝试使用我找到的代码片段附带的相同示例多边形,并且我注意到它适用于某些整体而不适用于其他人。要重现结果,这里是最小的示例代码:
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
都只是从应该有效的代码片段中复制而来。但是,只有 coords2
适合我。我已经使用 https://geojsonlint.com 查看两个多边形之间是否存在差异,多边形的惯用手/方向似乎不是有效的 GeoJSON。我不知道 shapely
是否关心,但颠倒顺序——然后 https://geojsonlint.com 说它是有效的 GeoJSON,并且它在 map 上显示多边形——不会改变错误。
因此,它适用于 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
,工作片段如下所示:
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
在做什么。因此我不想提供答案。
最佳答案
我认为你做得很好,只是反向操作不会创建新数据集。
尝试使用这个函数来创建倒序列表:
def rev_slice(mylist):
'''
return a revered list
mylist: is a list
'''
a = mylist[::-1]
return a
像这样执行函数:
coords = rev_slice(coords1)
关于python - 为什么在某些情况下转换形状多边形不起作用?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/63034619/