突出显示共享边/边的多边形

Question

我如何创建... ---我不知道；

list, dict

在 GeoDataFrame 中，共享边/边的所有多边形。多边形将相交但永远不会交叉。

import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt

polys = gpd.GeoSeries([Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
                     Polygon([(0,2), (2,2), (2,4), (0,4)]),
                     Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
                     Polygon([(3,3), (5,3), (5,5), (3,5)])])

fp = gpd.GeoDataFrame({'geometry': polys, 'name': ['a', 'b', 'c', 'd'],
                      'grnd': [25, 25, 25, 25],
                      'rf': [29, 35, 26, 31]})

fig, ax = plt.subplots(figsize=(5, 5))
fp.plot(ax=ax, alpha=0.3, cmap='tab10', edgecolor='k',)
fp.apply(lambda x: ax.annotate(text=x['name'], xy=x.geometry.centroid.coords[0], ha='center'), axis=1)
plt.show()

for i, row in fp.iterrows():
    oring = list(row.geometry.exterior.coords)#, row['ground_height']
    if row.geometry.exterior.is_ccw == False:
        #-- to get proper orientation of the normals
        oring.reverse() 
    for (j, v) in enumerate(oring[:-1]):
        print(oring[j][0], oring[j][1], oring[j+1][0], oring[j+1][1], row['name'])

预期结果：

0.0 0.0 2.0 0.0 a
2.0 0.0 2.0 1.5 a c
2.0 1.5 2.0 2.0 a 
2.0 2.0 0.0 2.0 a b
0.0 2.0 0.0 0.0 a
0.0 2.0 2.0 2.0 b a
2.0 2.0 2.0 4.0 b 
2.0 4.0 0.0 4.0 b
0.0 4.0 0.0 2.0 b... and so on

Answer 1

再看看预期的结果，我会先创建一个包含所有线段的系列，然后检查每个线段是否都接触多边形并且相交的长度是否大于0：

import geopandas as gpd
from shapely.geometry import Polygon, LineString
import matplotlib.pyplot as plt

polys = gpd.GeoSeries([
    Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
    Polygon([(0,2), (2,2), (2,4), (0,4)]),
    Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
    Polygon([(3,3), (5,3), (5,5), (3,5)])
])

fp = gpd.GeoDataFrame({
    'geometry': polys, 
    'name': ['a', 'b', 'c', 'd'],
    'grnd': [25, 25, 25, 25],
    'rf': [29, 35, 26, 31]
})

# Create series of all the line segments
lines = fp.geometry.apply(lambda x: list(map(
    LineString, 
    zip(x.boundary.coords[:-1], x.boundary.coords[1:]))
)).explode()

result = {
    str(line): list(fp.loc[
        (fp.geometry.touches(line)) # line touches the polygon
        & (fp.geometry.intersection(line).length > 0), # And the intersection is more than just a point
        'name'
    ].values)
    for line in lines
}

输出：

{'LINESTRING (0 0, 2 0)': ['a'],
 'LINESTRING (2 0, 2 1.5)': ['a', 'c'],
 'LINESTRING (2 1.5, 2 2)': ['a'],
 'LINESTRING (2 2, 0 2)': ['a', 'b'],
 'LINESTRING (0 2, 0 0)': ['a'],
 'LINESTRING (0 2, 2 2)': ['a', 'b'],
 'LINESTRING (2 2, 2 4)': ['b'],
 'LINESTRING (2 4, 0 4)': ['b'],
 'LINESTRING (0 4, 0 2)': ['b'],
 'LINESTRING (2 0, 5 0)': ['c'],
 'LINESTRING (5 0, 5 1.5)': ['c'],
 'LINESTRING (5 1.5, 2 1.5)': ['c'],
 'LINESTRING (2 1.5, 2 0)': ['a', 'c'],
 'LINESTRING (3 3, 5 3)': ['d'],
 'LINESTRING (5 3, 5 5)': ['d'],
 'LINESTRING (5 5, 3 5)': ['d'],
 'LINESTRING (3 5, 3 3)': ['d']}

编辑： 要处理

MultiPolygon

s 你可以尝试：

import json
from itertools import chain

import geopandas as gpd
from shapely.geometry import Polygon, LineString, MultiPolygon
import matplotlib.pyplot as plt

polys = gpd.GeoSeries([
    Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
    MultiPolygon([Polygon([(0,2), (2,2), (2,4), (0,4)])]),
    Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
    Polygon([(3,3), (5,3), (5,5), (3,5)])
])

fp = gpd.GeoDataFrame({
    'geometry': polys, 
    'name': ['a', 'b', 'c', 'd'],
    'grnd': [25, 25, 25, 25],
    'rf': [29, 35, 26, 31]
})

# Create series of all the line segments
lines = fp.geometry.apply(lambda x: (
    list(map(LineString, zip(x.boundary.coords[:-1], x.boundary.coords[1:])))
    if isinstance(x, Polygon)
    else list(chain(*list(list(map(
        LineString, zip(poly.boundary.coords[:-1], poly.boundary.coords[1:])
    )) for poly in x.geoms)))
)).explode()

result = {
    str(line): list(fp.loc[
        (fp.geometry.touches(line)) # line touches the polygon
        & (fp.geometry.intersection(line).length > 0), # And the intersection is more than just a point
        'name'
    ].values)
    for line in lines
}

Answer 2

你在找

geopandas.GeoSeries.touches

吗？

您可以执行以下操作：

import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt

polys = gpd.GeoSeries([
    Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
    Polygon([(0,2), (2,2), (2,4), (0,4)]),
    Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
    Polygon([(3,3), (5,3), (5,5), (3,5)])
])

fp = gpd.GeoDataFrame({
    'geometry': polys, 
    'name': ['a', 'b', 'c', 'd'],
    'grnd': [25, 25, 25, 25],
    'rf': [29, 35, 26, 31]
})

result = {
    row['name']: {
        'touching polygons': list(fp.loc[fp.geometry.touches(row.geometry), 'name'].values)
    }
    for i, row in fp.iterrows()
}

输出：

{
    'a': {'touching polygons': ['b', 'c']}, 
    'b': {'touching polygons': ['a']}, 
    'c': {'touching polygons': ['a']}, 
    'd': {'touching polygons': []}
}

Answer 3

看看

Expected result

，我可以说，没有可用的空间谓词来操作并获得那种结果。

Intersects

或

touches

在两个几何图形之间会得到一些假阳性结果。

在这里，我实施了一个简单的检查，

same_lineQ(x1y1, x2y2)

，作为在需要缺少空间操作的地方使用的过程。

# PART 1
import geopandas as gpd
from shapely.geometry import Polygon, LineString, Point
import matplotlib.pyplot as plt
import pandas as pd

polys = gpd.GeoSeries([
    Polygon([(0,0), (2,0), (2, 1.5), (2,2), (0,2)]),
    Polygon([(0,2), (2,2), (2,4), (0,4)]),
    Polygon([(2,0), (5,0), (5,1.5), (2,1.5)]),
    Polygon([(3,3), (5,3), (5,5), (3,5)])
])

fp = gpd.GeoDataFrame({
    'geometry': polys, 
    'name': ['a', 'b', 'c', 'd'],
    'grnd': [25, 25, 25, 25],
    'rf': [29, 35, 26, 31]
})

fig, ax = plt.subplots(figsize=(5/2, 5/2))
fp.plot(ax=ax, alpha=0.3, cmap='tab10', edgecolor='k',)
fp.apply(lambda x: ax.annotate(text=x['name'], xy=x.geometry.centroid.coords[0], ha='center'), axis=1)

# Part 2
# Collect all the line segments from all polygons
def get_all_xy0xy1(p1, attrib="none"):
    """
    p1: a Polygon object, has single `exterior`
    returns: list of [[x0,y0],[x1,y1]] ready for LineString creation
    eg.: LineString([(0, 0), (9, 9)])
    """
    xy0_xy1_list = []
    attribs = []
    for ix,xy in enumerate(zip(p1.exterior.xy[0], p1.exterior.xy[1])):
        # 3 or more items
        #print(ix,xy)  #either x, or y separately
        if ix>0:
            #print([prev, xy]) #list of x,y; from-to
            xy0_xy1_list.append([prev, xy])
            attribs.append(attrib)
        prev = xy
    return xy0_xy1_list, attribs

# Line segments are collected in `all_line_segs`
all_line_segs = []
names = []
for ix, row in fp.iterrows():
    name = row['name']
    geom = row.geometry
    all_line_segs += get_all_xy0xy1(geom, name)[0]
    names += get_all_xy0xy1(geom, name)[1]

# Create a dataframe using the line segments
line_segs = pd.DataFrame({
    'xy1_xy2': all_line_segs, 
    'name': names
})

# Part 3
def same_lineQ(x1y1, x2y2):
    """
    Input: x1y1, x2y2; two list of (x,y).
    Returns:
     True if they represent the same LineString
       ignoring the direction
     else returns False.
    """
    return (x1y1[0] in x2y2) and (x1y1[1] in x2y2)

for ir, irow in line_segs.iterrows():
    iname = irow['name']
    ixys = irow['xy1_xy2']
    targets = set()
    for kr, krow in line_segs.iterrows():
        if ir != kr:
            kname = krow['name']
            if iname != kname:
                kxys = krow['xy1_xy2']

                if same_lineQ(ixys, kxys)==True:
                    #print(ixys, kxys, same_lineQ(ixys, kxys))
                    targets.update(kname)

        else:
            pass
    if len(targets)==0:
        print(ixys[0][0],ixys[0][1], ixys[1][0],ixys[1][1], iname)
    else:    
        print(ixys[0][0],ixys[0][1], ixys[1][0],ixys[1][1], iname, targets.pop())

输出：

0.0 0.0 2.0 0.0 一个
2.0 0.0 2.0 1.5 交流
2.0 1.5 2.0 2.0 一
2.0 2.0 0.0 2.0 一 b
0.0 2.0 0.0 0.0 一个
0.0 2.0 2.0 2.0 乙
2.0 2.0 2.0 4.0 b
2.0 4.0 0.0 4.0 b
0.0 4.0 0.0 2.0 乙
2.0 0.0 5.0 0.0 c
5.0 0.0 5.0 1.5℃
5.0 1.5 2.0 1.5℃
2.0 1.5 2.0 0.0
3.0 3.0 5.0 3.0 天
5.0 3.0 5.0 5.0 天
5.0 5.0 3.0 5.0 天
3.0 5.0 3.0 3.0 天

编辑

如果geodataframe比较复杂，比如有些多边形有洞，或者有些行有MultiPolygon而不是Polygon，上面的代码就不行了。它仅适用于具有仅无孔的多边形的地理数据框。

遇到这种情况怎么办？

一种方法是分解没有简单多边形的行，即具有带孔的多边形或多边形，并获得具有所有简单多边形的结果地理数据框。

以这个修改后的geodataframe为例：

（注意

e

是2个多边形的MultiPolygon，其中一个是单孔。而

f

是简单的多边形。）

如果与上面的代码一起使用，输出将是：

0.0 0.0 2.0 0.0 一个
2.0 0.0 2.0 1.5 交流
2.0 1.5 2.0 2.0 一
2.0 2.0 0.0 2.0 一 b
0.0 2.0 0.0 0.0 一个
0.0 2.0 2.0 2.0 乙
2.0 2.0 2.0 4.0 b
2.0 4.0 0.0 4.0 b
0.0 4.0 0.0 2.0 乙
2.0 0.0 5.0 0.0 c
5.0 0.0 5.0 1.5℃
5.0 1.5 2.0 1.5℃
2.0 1.5 2.0 0.0
3.0 3.0 5.0 3.0 天
5.0 3.0 5.0 5.0 d e
5.0 5.0 3.0 5.0 天
3.0 5.0 3.0 3.0 天
6.5 0.0 7.5 0.0 f
7.5 0.0 7.5 1.5 楼
7.5 1.5 6.5 1.5 铁
6.5 1.5 6.5 0.0 f
6.0 1.0 8.0 1.0 电子
8.0 1.0 8.0 3.0 电子
8.0 3.0 6.0 3.0 电子
6.0 3.0 6.0 1.0 电子
5.0 3.0 8.0 3.0 电子
8.0 3.0 8.0 5.0 电子
8.0 5.0 5.0 5.0 电子
5.0 5.0 5.0 3.0 日
6.5 2.5 7.5 2.5 电子
7.5 2.5 7.5 1.5 电子
7.5 1.5 6.5 1.5 英法
6.5 1.5 6.5 2.5 电子