Cartopy:绘制移除国家边界的海岸线

Cartopy: Drawing the coastlines with a country border removed

我想用一种颜色勾勒出中国的轮廓,同时用另一种颜色展示全球的海岸线。我第一次这样做的尝试如下:

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

我把海岸线画得很厚,这样你就可以看到它们与国家边界重叠的事实。

我的问题是:有没有办法去掉国家轮廓旁边的海岸线,这样两条线就不会在视觉上相互影响?

Note: This question was asked to me directly via email, and I chose to post my response here so that others may learn/benefit from a solution.

Natural Earth 集合中没有名为 "coastlines without the Chinese border" 的数据集,因此我们将不得不自己制作它。要做到这一点,我们需要使用 shapely 操作,特别是 difference 方法。

差分法如下图所示(摘自Shapely's docs)。下面突出显示了两个示例圆圈(ab)的区别:

那么,目的就是达到写作的目的coastline.difference(china),并将这个结果可视化为我们的海岸线。

有多种方法可以做到这一点。 GeoPandas 和 Fiona 是两种可以提供非常可读的结果的技术。不过在这种情况下,让我们使用 cartopy 提供的工具:

首先,我们掌握了中国边界(另见:cartopy shapereader docs)。

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]

当我们将其可视化时,我们发现差异并不十分完美:

我们不想要黑线的原因是自然地球海岸线数据集与国家数据集的派生方式不同,因此它们不是完全重合的坐标。

为了解决这个问题,可以对中国边界应用一个小的 "hack" 以扩大边界以达到此交叉点的目的。 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]

有了这个 "hack",我得到了以下结果(包含完整代码以确保完整性):

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