在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()

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

现在的问题是，实际上我不想不加选择地将所有区域都遮住浅蓝色。我 - 理想情况下 - 想要确定受治疗国家的边界​​长度，然后根据您与一个或多个受治疗国家共享的边界数量来缩小您的受影响程度。

至少，我希望能够扔掉那些只与另一个国家分享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

由此产生的情节：