Cartopy:删除国家边界的海岸线

pel*_*son 4 python shapely cartopy

我想用一种颜色绘制中国的轮廓,同时在另一种颜色中显示全球海岸线.我这样做的第一次尝试如下:

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()
Run Code Online (Sandbox Code Playgroud)

结果

我把海岸线做得很厚,这样你就可以看到它们与国界重叠.

我的问题是:有没有办法去除国家大纲旁边的海岸线,这样我就没有两条线在视觉上相互交互?

注意:这个问题是通过电子邮件直接向我询问的,我选择在此发布我的回复,以便其他人可以从解决方案中学习/受益.

pel*_*son 11

自然地球系列中没有名为"没有中国边界的海岸线"的数据集,因此我们将不得不自己制造它.为此,我们将要使用匀称操作,特别是它的difference方法.

差异方法如下图所示(取自Shapely的文档).下面突出显示两个示例圆(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.')
Run Code Online (Sandbox Code Playgroud)

接下来,我们掌握了海岸线几何形状:

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

coastlines = shapereader.Reader(coast).geometries()
Run Code Online (Sandbox Code Playgroud)

现在,把中国赶出海岸线:

coastlines_m_china = [geom.difference(china)
                      for geom in coastlines]
Run Code Online (Sandbox Code Playgroud)

当我们想象这一点时,我们发现差异并不完美:

不完美的差异

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

为了解决这个问题,可以在中国边境采用一个小的"黑客"来扩大边界,以达到这个交叉点的目的.该缓冲区的方法是非常适合这个用途.

# 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]
Run Code Online (Sandbox Code Playgroud)

有了这个"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()
Run Code Online (Sandbox Code Playgroud)

结果