Картопия: рисование береговых линий без границы страны.

Я хочу нарисовать контур Китая одним цветом, а также показать береговые линии мира другим. Моя первая попытка сделать это была следующая:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feature
import cartopy.io.shapereader as shapereader


countries = shapereader.natural_earth(resolution='110m',
                                      category='cultural',
                                      name='admin_0_countries')

# Find the China boundary polygon.
for country in shapereader.Reader(countries).records():
    if country.attributes['su_a3'] == 'CHN':
        china = country.geometry
        break
else:
    raise ValueError('Unable to find the CHN boundary.')


plt.figure(figsize=(8, 4))
ax = plt.axes(projection=ccrs.PlateCarree())

ax.set_extent([50, 164, 5, 60], ccrs.PlateCarree())

ax.add_feature(feature.LAND)
ax.add_feature(feature.OCEAN)
ax.add_feature(feature.COASTLINE, linewidth=4)

ax.add_geometries([china], ccrs.Geodetic(), edgecolor='red',
                  facecolor='none')

plt.show()

Результат

Я сделал береговые линии толстыми, чтобы вы могли видеть тот факт, что они пересекаются с границей страны.

У меня вопрос: есть ли способ удалить береговую линию рядом с контуром страны, чтобы две линии не взаимодействовали друг с другом визуально?

Примечание. Этот вопрос был задан мне непосредственно по электронной почте, и я решил опубликовать свой ответ здесь, чтобы другие могли узнать / извлечь выгоду из решения.


person pelson    schedule 14.07.2017    source источник


Ответы (1)


В коллекции Natural Earth нет набора данных под названием «Береговые линии без китайской границы», поэтому нам придется изготовить его самостоятельно. Для этого мы захотим использовать фигурные операции и, в частности, метод difference.

Метод различия показан ниже (взят из документации Shapely) наглядно. Разница в двух примерах кругов (a и b) выделена ниже:

метод различия

Итак, цель состоит в том, чтобы перейти к написанию coastline.difference(china) и визуализировать этот результат как наши береговые линии.

Есть несколько способов сделать это. GeoPandas и Fiona - две технологии, которые дадут очень удобочитаемые результаты. Но в этом случае давайте воспользуемся инструментами, которые предоставляет cartopy:

Сначала мы захватываем границу с Китаем (см. Также: cartopy shapereader документы).

import cartopy.io.shapereader as shapereader


countries = shapereader.natural_earth(resolution='110m',
                                      category='cultural',
                                      name='admin_0_countries')

# Find the China boundary polygon.
for country in shapereader.Reader(countries).records():
    if country.attributes['su_a3'] == 'CHN':
        china = country.geometry
        break
else:
    raise ValueError('Unable to find the CHN boundary.')

Далее мы получаем геометрию береговых линий:

coast = shapereader.natural_earth(resolution='110m',
                                  category='physical',
                                  name='coastline')

coastlines = shapereader.Reader(coast).geometries()

А теперь уберите Китай с берегов:

coastlines_m_china = [geom.difference(china)
                      for geom in coastlines]

Когда мы визуализируем это, мы видим, что разница не совсем идеальная:

«Несовершенная

Причина появления черных линий там, где они нам не нужны, заключается в том, что набор данных береговых линий Natural Earth получен иначе, чем набор данных стран, и поэтому они не являются полностью совпадающими координатами.

Чтобы обойти этот факт, можно применить небольшой «взлом» к границе Китая, чтобы расширить границу для целей этого пересечения. Для этой цели идеально подходит метод buffer.

# Buffer the Chinese border by a tiny amount to help the coordinate
# alignment with the coastlines.
coastlines_m_china = [geom.difference(china.buffer(0.001))
                      for geom in coastlines]

С этим «хаком» я получаю следующий результат (полный код включен для полноты):

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as feature
import cartopy.io.shapereader as shapereader


coast = shapereader.natural_earth(resolution='110m',
                                  category='physical',
                                  name='coastline')

countries = shapereader.natural_earth(resolution='110m',
                                      category='cultural',
                                      name='admin_0_countries')

# Find the China boundary polygon.
for country in shapereader.Reader(countries).records():
    if country.attributes['su_a3'] == 'CHN':
        china = country.geometry
        break
else:
    raise ValueError('Unable to find the CHN boundary.')

coastlines = shapereader.Reader(coast).geometries() 

# Hack to get the border to intersect cleanly with the coastline.
coastlines_m_china = [geom.difference(china.buffer(0.001))
                      for geom in coastlines]

ax = plt.axes(projection=ccrs.PlateCarree())

ax.set_extent([50, 164, 5, 60], ccrs.PlateCarree())
ax.add_feature(feature.LAND)
ax.add_feature(feature.OCEAN)

ax.add_geometries(coastlines_m_china, ccrs.Geodetic(), edgecolor='black', facecolor='none', lw=4)
ax.add_geometries([china], ccrs.Geodetic(), edgecolor='red', facecolor='none')

plt.show()

«Результат»

person pelson    schedule 14.07.2017