Search code examples
matplotlibcartopy

How to apply Earth Features and Land/Ocean masks in high resolution coastlines in Cartopy?

I am using the coastlines of the GSHHS dataset in Cartopy. This has a high resolution for coastlines. But I want to not only plot the high resolution coastline but also apply a mask for the ocean.

import matplotlib.pyplot as plt
import cartopy

fig = plt.figure(figsize=(20,12))
ax = plt.axes(projection=cartopy.crs.PlateCarree())
coast = cartopy.feature.GSHHSFeature(scale="full")
ax.add_feature(coast, linewidth=2)
ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m"))
ax.set_extent([-17, -16, 27.9, 28.7])

Executing the code there're differences in the images, since I guess that ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m")) is using the "10m" resolution, while GSHHS has a higher resolution.

How can I mask using GSHHS higher resolution? Thx.

Solution

  • Before one can answer the question how to apply a mask to hide features in the main plot, we need to investigate the available masks first.

    In our case, the main plot is Natural_Earth 10m resolution Physical Land features, and various resolutions of GSHHSFeature as the available masks.

    The code and the output plot below reveals the insight.

    # Code adapted from:-
# Src: https://ctroupin.github.io/posts/2019-09-02-fine-coast/
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

resolutions = {"c": "crude",
               "l": "low",
               "i": "intermediate",
               "h": "high",
               "f": "full"}
coordinates = (8.7, 8.81, 42.55, 42.60)
myproj = ccrs.PlateCarree()
fig = plt.figure(figsize=(8, 4))
for i, res in enumerate(resolutions):
    ax = plt.subplot(2, 3, i+1, projection=myproj)
    coast = cfeature.GSHHSFeature(scale=res)
    ax.add_feature(coast, facecolor="lightgray")
    ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)
    ax.set_xlim(coordinates[0], coordinates[1])
    ax.set_ylim(coordinates[2], coordinates[3])
    plt.title(resolutions[res])
plt.suptitle("GSHHS: gray Versus 10m_Physical_Land: yellow/red")
plt.show()

    Suppose we need a plot at this zoom level. It is clearly that the outlines from 2 data sources do not fit well enough to the eyes of the viewers. We may conclude that none of the available masks is fit for the target plot.

    But if the plot extents is wider, or smaller scale plots, coupled with some cartographic techniques, e.g. using thicker coastlines, one may get acceptable plots. The process is trial-and-error approach.

    gnshhh

    Edit1

    With (Global_land_mask) added, more choices can be plotted for comparison.

    from global_land_mask import globe
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

# Extent of map in degrees
minlon,maxlon,minlat,maxlat = (8.7, 8.81, 42.55, 42.60)

# Lat/lon points to get for `global_land_mask` uses
# Finer than 500x250 has no improvement 
lons = np.linspace(minlon,maxlon, 500)
lats = np.linspace(minlat,maxlat, 250)

# Make a grid 
lon_grid, lat_grid = np.meshgrid(lons,lats)

# Get whether the points are on land.
z = globe.is_land(lat_grid, lon_grid)

# GSHHS ... 
resolutions = {"c": "crude",
               "l": "low",
               "i": "intermediate",
               "h": "high",
               "f": "full"}

myproj = ccrs.PlateCarree()
fig = plt.figure(figsize=(8, 4))
for i, res in enumerate(resolutions):
    ax = plt.subplot(2, 3, i+1, projection=myproj)

    # GSHHSFeature
    coast = cfeature.GSHHSFeature(scale=res)
    ax.add_feature(coast, facecolor="brown", alpha=0.5)

    # 10m physical_land
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "land", "10m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)

    # Global_land_mask data is used to create fillcontour
    # The fillcontour with proper (colormap, zorder, alpha) can be used as land `mask`
    ax.contourf(lon_grid, lat_grid, z, cmap="Greys_r", alpha=0.4)

    ax.set_xlim(minlon, maxlon)
    ax.set_ylim(minlat, maxlat)
    plt.title(resolutions[res])

plt.suptitle("GSHHS:brown/black | 10m_Land:yellow/red | Global_land_mask:light_gray")
plt.show()

# The best resolutuion from `Global_land_mask` is plotted in `lightgray` covering the sea areas

    compare3a