Cartopy: Drawing the coastlines with a country border removed

Question

I want to draw the outline of China in one color whilst also showing the global coastlines in another. My first attempt at doing this was as follows:

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

I've made the coastlines thick so that you can see the fact that they are overlapping with the country border.

My question is: Is there a way to remove the coastline next to the country outline so that I don't have the two lines interacting with one-another visually?

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.

Answer 1

There is no dataset in the Natural Earth collection called "coastlines without the Chinese border", so we are going to have to manufacture it ourselves. To do this we will want to use the shapely operations, and in particular it's difference method.

The difference method is shown below (taken from Shapely's docs) pictorially. The difference of two example circles (a and b) are highlighted below:

The aim then, is to get to the point of writing coastline.difference(china), and to visualise this result as our coastlines.

There are a number of ways of doing this. GeoPandas and Fiona are two technologies that would give very readable results. In this case though, let's use the tools that cartopy provides:

First, we get hold of the China border (see also: 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.')

Next, we get hold of the coastlines geometries:

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

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

Now, take China out of the coastlines:

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

When we visualise this, we see that the difference isn't quite perfect:

The reason for the black lines where we didn't want them is that the Natural Earth coastlines dataset is derived differently to the countries dataset, and they are therefore not perfectly coincident coordinates.

To get around this fact, a small "hack" can be applied to the China border to expand the boundary for the purposes of this intersection. The buffer method is perfect for this purpose.

# 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]

With this "hack" in place, I get the following result (full code included for completeness):

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