标签: shapely

我不知道如何使用 Shapely 在 Python 中添加两个多边形。

通过添加，我的意思是，例如，如果我添加两个高度为 4 和宽度为 2 的正方形，并且它们具有相同的坐标，则它应该返回一个高度为 8 和宽度为 2 的正方形。

我尝试过使用 MultiPolygons 并使用两个多边形之间的并集，但我无法获得所需的累积高度结果。

有人知道如何做我所描述的吗？或者还有其他 Python 模块可以让我做同样的事情吗？

我正在尝试验证平面上的一些多边形is_valid，但我得到 Too few points in geometry component at or near point的多边形的 z 不是恒定的。

有没有办法验证这些其他多边形？

这是一个例子：

from shapely.geometry import Polygon


poly1 = Polygon([(0,0), (1,1), (1,0)])
print(poly1.is_valid)
# True

# z=1
poly2 = Polygon([(0,0,1), (1,1,1), (1,0,1)])
print(poly2.is_valid)
# True

# x=1
poly3 = Polygon([(1,0,0), (1,1,1), (1,1,0)])
print(poly3.is_valid) 
# Too few points in geometry component at or near point 1 0 0
# False

我正在尝试在距其他坐标最近的点处拆分 Shapely LineString。我可以使用 获得线上最近的点project，interpolate但此时我无法分割线，因为它不是顶点。

我需要沿着边缘分割线，而不是捕捉到最近的顶点，以便最近的点成为线上的新顶点。

这是我到目前为止所做的：

from shapely.ops import split
from shapely.geometry import Point, LineString

line = LineString([(0, 0), (5,8)])
point = Point(2,3)

# Find coordinate of closest point on line to point
d = line.project(point)
p = line.interpolate(d)
print(p)
# >>> POINT (1.910112359550562 3.056179775280899)

# Split the line at the point
result = split(line, p)
print(result)
# >>> GEOMETRYCOLLECTION (LINESTRING (0 0, 5 8))

谢谢！

我的问题是：如何确定已分割的旋转矩形几何图形的哪一边Aside和哪一边是分割该几何图形的任意边的“左”和“右”边？BsideLineString

LineString出于此问题的目的，“左”和“右”被定义为按顺序从一个节点“行走”到另一个节点时分离器的左侧和右侧。

我创建了这个函数，用于将任意几何图形（非集合）分割为两侧shapely- “左”和“右”：

import shapely.geometry as geo
import shapely.ops as ops

def splitLR(geom, splitter):
    """Split a geometry into a 'left' and 'right' side using the shapely API"""
    if not isinstance(splitter, geo.LineString):
        raise TypeError("The splitter must be a LineString")
    if not splitter.is_simple:
        raise ValueError("Only simple splitter objects allowed")
    if hasattr(geom, "__iter__"):
        raise ValueError("Geometry collections not allowed")
    geom_extents = geo.GeometryCollection([geom, splitter]).minimum_rotated_rectangle
    sides = ops.split(geom_extents, splitter)
    try:
        Aside, Bside = sides …

我正在使用 Python 中的 Shapely 库。我找到了两条线的交点，返回值作为 MultiPoint 对象给出。

如何解构对象以获得交点中的各个点？

这是代码：

from shapely.geometry import LineString, MultiLineString
a = LineString([(0, 1), (0, 2), (1, 1), (2, 0)])
b = LineString([(0, 0), (1, 1), (2, 1), (2, 0)])
x = a.intersection(b)

输出：

print(x) 
MULTIPOINT (1 1, 2 0)

因此，在这种情况下，我会寻找一种方法来提取交点 (1,1) 和 (2,0)。

上一个问题“使用 cartopy 检查地理坐标点是陆地还是海洋”的答案请参阅https://stackoverflow.com/questions/asktitle=Checking%20if%20a%20geocooperative%20point%20is%20land%20or%20ocean% 20using%20cartopy%20and%20Natural%20Earth%2010m%20data建议使用以下代码来确定地理坐标是否为“is_land”：

    import cartopy.io.shapereader as shpreader
    import shapely.geometry as sgeom
    from shapely.ops import unary_union
    from shapely.prepared import prep

    land_shp_fname = shpreader.natural_earth(resolution='50m',
                                   category='physical', name='land')

    land_geom = 
    unary_union(list(shpreader.Reader(land_shp_fname).geometries()))
    land = prep(land_geom)

    def is_land(x, y):
       return land.contains(sgeom.Point(x, y))

当自然地球“物理”“土地”形状文件的分辨率更改为“10m”时，此代码返回地理坐标 (0,0) 的意外结果“True”。

    >>> print(is_land(0, 0))
    True

这是 Natural Earth shapefile 数据或 shapely 实用程序代码的问题吗？

    print shapely.__version__
    1.6.4.post1

如何仅从数组中提取值而不使用文本“array”和“typecode”？

数组是shapely.linestring.centroid.xy：

a = LineString.centroid.xy
print(a)
>> (array('d', [-1.72937...45182697]), array('d', [2.144161...64685937]))
print(a[0])
>> array('d', [-1.7293720645182697])

我只需要作为-1.7293...浮点数而不是整个数组业务。

我有一个 geopandas 数据框 gdf

gdf
    ID  longitude   latitude    geometry
0   80  103.619501  1.2810      POINT (103.619500987 1.281)
1   81  103.619501  1.2855      POINT (103.619500987 1.2855)

按照这个建议，我在它周围创建了一个方形缓冲区，其距离bd定义为：

bd = abs((gdf['latitude'][0]-gdf['latitude'][1])/2)

最后我能够得到以下信息：

buffer = gdf.buffer(bd)
envelope = buffer.envelope

f, ax = plt.subplots(figsize=(7.5, 7.5))
envelope.plot(color='white', edgecolor='gray',ax=ax)
gdf.plot(ax=ax)

如何设置bd对应于 500 米的距离？

我在这里找到了荷兰 4 位邮政编码的 shapefile：

https://public.opendatasoft.com/explore/dataset/openpostcodevlakkenpc4/export/?location=9,52.16172,5.56595&basemap=jawg.streets

我想做的是将共享前两位数字的邮政编码组合起来，然后将它们绘制在地图上。

不幸的是，数据似乎有问题——荷兰的一些地区似乎没有被覆盖。结果，当我组合多边形以获取 2 位数多边形（从 4 位数多边形）时，生成的多边形中存在孔。

我需要填补这些漏洞。我看过类似问题的帖子，但似乎没有什么能完全满足我的需要。特别是，我看到一篇关于使用凹壳的帖子，但这在这里似乎有点矫枉过正。

我已经设法修复了几个工件（例如“条子”），但漏洞仍然存在。

这是我到目前为止所拥有的：

from shapely.ops import cascaded_union
from shapely.geometry import JOIN_STYLE, Polygon, MultiPolygon

import os
import folium
import pandas as pd
import geopandas as gpd

GEOM_LOC = r"PATH_TO_FILE_ABOVE\openpostcodevlakkenpc4.shx"

# Get the data
geom = gpd.read_file(GEOM_LOC)

# Remove empty or nan records
is_empty = geom.geometry.is_empty
is_nan = geom.geometry.isna()
geom = geom[~(is_empty | is_nan)]

# Add field for 2 digit postcode
geom["digit"] = geom.pc4.apply(lambda x: x[0:2])
geom = geom[["digit", "geometry"]]

# …

这是一个非常简单的案例，但到目前为止我还没有找到任何简单的方法。这个想法是获得 a 中定义的所有点GeoDataFrame和 another 中定义的点之间的一组距离GeoDataFrame。

import geopandas as gpd
import pandas as pd

# random coordinates
gdf_1 = gpd.GeoDataFrame(geometry=gpd.points_from_xy([0, 0, 0], [0, 90, 120]))
gdf_2 = gpd.GeoDataFrame(geometry=gpd.points_from_xy([0, 0], [0, -90]))
print(gdf_1)
print(gdf_2)

#  distances are calculated elementwise
print(gdf_1.distance(gdf_2))

这会产生 ingdf_1和gdf_2共享相同索引的点之间的元素距离（还有一个警告，因为两个 GeoSeries 没有相同的索引，这就是我的情况）。

                geometry
0    POINT (0.000 0.000)
1   POINT (0.000 90.000)
2  POINT (0.000 120.000)
                    geometry
0    POINT (0.00000 0.00000)
1  POINT (0.00000 -90.00000)
/home/seydoux/anaconda3/envs/chelyabinsk/lib/python3.8/site-packages/geopandas/base.py:39: UserWarning: The indices of the two GeoSeries …