Возвращает длину отрезка границы между географическими областями в геопандах

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

Я думаю, вам нужен пример, который демонстрирует правильные геопространственные операции для получения результата.Мой код ниже показывает, как сделать 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

Полученный сюжет: