geopandas中地理区域之间边界段的返回长度

Question

我是使用geopandas 的新手，所以我有一个相当基本的问题。我想确定地理数据框中相邻地点之间发生了多少边界接触。

我将提供一个例子。以下代码读取预加载的地理框架，随机创建标记为“已处理”的国家/地区，定义一个为其邻国提供函数的函数，然后将结果绘制成与边界国家/地区稍浅的阴影。

import geopandas as gp
import numpy as np
import matplotlib.pyplot as plt
path = gp.datasets.get_path('naturalearth_lowres')
earth = gp.read_file(path)
africa = earth[earth.continent=='Africa']
africa['some_places'] = np.random.randint(0,2,size=africa.shape[0])*2

# Define and apply a function that determines which countries touch which others
def touches(x):
    result = 0
    if x in africa.loc[africa.some_places==2,'geometry']:
        result = 2
    else:
        for y in africa.loc[africa.some_places==2,'geometry']:
            if y.touches(x) :
                result = 1
                break
            else:
                continue
    return result
africa['touch'] = africa.geometry.apply(touches)

# Plot the main places which are 2, the ones that touch which are 1, and the non-touching 0
fig, ax = plt.subplots()
africa.plot(column='touch', cmap='Blues', linewidth=0.5, ax=ax, edgecolor='.2')
ax.axis('off')
plt.show()

对我来说，这给出了以下地图：

现在的问题是，实际上我不想不加选择地将所有区域都涂成浅蓝色。我 - 理想情况下 - 想要确定受处理国家/地区边界的长度，然后根据您与一个或多个受处理国家共享的边界量来计算您受到的影响程度。 p>

至少，我希望能够丢弃与另一个国家只有 1 或 2 英里边界的地方（或者可能在拐角处相遇）。欢迎任何建议或解决方案！

Answer 1

我认为您需要一个示例来演示正确的地理空间操作以获得结果。我下面的代码显示了如何在 N 和 S 美洲大陆之间执行intersection，获取相交线，然后计算其长度。最后，在地图上画线。我希望它对您的项目有用且适用。

import geopandas
import numpy as np
import matplotlib.pyplot as plt

# make use of the provided world dataset
world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

# drop some areas (they crash my program)
world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]

# create a geodataframe `wg`, taking a few columns of data
wg = world[['continent', 'geometry', 'pop_est']]
# create a geodataframe `continents`, as a result of dissolving countries into continents
continents = wg.dissolve(by='continent')

# epsg:3857, Spherical Mercator
# reproject to `Spherical Mercator` projection
continents3857 = continents.to_crs(epsg='3857')

# get the geometry of the places of interest
namerica = continents3857.geometry[3]   # north-america
samerica = continents3857.geometry[5]   # south-america

# get intersection between N and S America continents
na_intersect_sa = namerica.intersection(samerica)   # got multi-line

# show the length of the result (multi-line object)
blength = na_intersect_sa.length   # unit is meters on Spherical Mercator
print("Length in meters:", "%d" % blength)

# The output:
# Length in meters: 225030

ax = continents3857.plot(column="pop_est", cmap="Accent", figsize=(8,5))

for ea in na_intersect_sa.__geo_interface__['coordinates']:
    #print(ea)
    ax.plot(np.array(ea)[:,0], np.array(ea)[:,1], linewidth=3, color="red")

ax.set_xlim(-9500000,-7900000)
ax.set_ylim(700000, 1400000)

xmin,xmax = ax.get_xlim()
ymin,ymax = ax.get_ylim()

rect = plt.Rectangle((xmin,ymin), xmax-xmin, ymax-ymin, facecolor="lightblue", zorder=-10)
ax.add_artist(rect)

plt.show()   # sometimes not needed

结果图：