Lakes filled with land feature

[SUPScientist]

Hello—I'm trying to plot ocean/lake/water quality/color data (phytoplankton chlorophyll, for example) in the Great Lakes and shade in the land feature with facecolor gray. This shades the lakes also and water quality data won't show up in the lakes. A few questions: (1) is this the expected behavior? If not, any recommendations for how to fix this? I've tested up to 10m resolution. (2) If so, is there a way to do something like a Boolean NOT to do a land minus lakes feature? Thanks for any suggestions!

[greglucas]

Can you provide a small self contained example of code you've tried that doesn't work?

If I am understanding the question though, yes the land features have the lakes and rivers filled in, derived from the coastlines dataset. https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-land/

So, a few options come to mind. You could get the features (land and lakes) as shapely geometries and difference them to create a new geometry. Or perhaps simpler would be to plot the land, then plot the lakes in white, then plot whatever you need over that white area.

[SUPScientist]

Thanks for the help, @greglucas. Yes, it sounds like you are understanding the question.

Here's the relevant code:

    # Add coastline and fill continents with gray and use black for edge (coastline)
    # are prescribed resolution for the coast (mapres)...
    # ---
    feature = cfeature.NaturalEarthFeature(name='land', category='physical', scale=mapres, edgecolor='k', facecolor='gray')
    ax.add_feature(feature)

    lake_feature = cfeature.NaturalEarthFeature('physical', 'lakes', '10m', edgecolor='blue', facecolor='white')
    ax.add_feature(lake_feature)

    # NOTE: PJB added in 2023 code to try to show data in (Great) Lakes
    # add lakes and repeat same procedure if needed for rivers

    # Add ocean color or temperature data, stored in geophys_img
    if scale_type == 'LIN': sat_img= ax.imshow(geophys_img,  transform=ccrs.PlateCarree(), extent=map_data_extent, cmap=mycmap, vmin=low_limit, vmax=upper_limit)
    if scale_type == 'LOG': sat_img= ax.imshow(geophys_img,  transform=ccrs.PlateCarree(), extent=map_data_extent, cmap=mycmap, norm=matplotlib.colors.LogNorm(low_limit,upper_limit))

    plt.colorbar(sat_img,location='bottom')

Output:

When I do only the imshow (i.e., no land/lake features), there are data. Another workaround is to turn off the land feature altogether and use only edge color for the lake feature; then the geophys_img shows up, but this isn't the desired aesthetic.

One final note: I'd strongly prefer to plot the imshow of the geophys_img first, then layer the basemaps/features on top as there can be an abundance of flagged values outside of the bounds of the lakes which I would like to cover with a land feature, while leaving the flagged values inside the lakes visible (# specify color map with NaN's as black... mycmap.set_bad('k')).

I think your recommendation to start with shapely geometries and difference them may be the best way to go; if you have any example code you can recommend, I'd appreciate it, but understood if that's asking too much!

[SUPScientist]

Here's a better, self-contained toy example.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors

import scipy.misc

import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy
import cartopy.feature as cfeature
from cartopy.io import shapereader
from shapely.geometry import MultiPolygon, Polygon

# Create a Cartopy projection
projection = ccrs.PlateCarree()
mapres = '10m'

# Desired bounds
xmin = -83.96
xmax = -77.57
ymin = 40.18
ymax = 43.34

# Create a figure and axes for the plot
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=projection))

# Add land (automatically including lakes)
land_feature = cfeature.NaturalEarthFeature(name='land', category='physical', scale=mapres, edgecolor='k', facecolor='gray')
ax.add_feature(land_feature)

# Add lakes
lake_feature = cfeature.NaturalEarthFeature(name='lakes', category='physical', scale=mapres, edgecolor='blue', facecolor='none')
ax.add_feature(lake_feature)

# Generate your layer data (assuming it's stored in 'layer_data' variable)
layer_data = np.random.rand(100, 100)

# Plot the layer using imshow
im = ax.imshow(layer_data, extent=[xmin, xmax, ymin, ymax], transform=projection, origin='lower', cmap='viridis')

plt.colorbar(im,location='bottom')

plt.show()

Produces this:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors

import scipy.misc

import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy
import cartopy.feature as cfeature
from cartopy.io import shapereader
from shapely.geometry import MultiPolygon, Polygon

# Create a Cartopy projection
projection = ccrs.PlateCarree()
mapres = '10m'

# Desired bounds
xmin = -83.96
xmax = -77.57
ymin = 40.18
ymax = 43.34

# Create a figure and axes for the plot
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=projection))

# Add land (automatically including lakes)
land_feature = cfeature.NaturalEarthFeature(name='land', category='physical', scale=mapres, edgecolor='k', facecolor='gray')
# ax.add_feature(land_feature)

# Add lakes
lake_feature = cfeature.NaturalEarthFeature(name='lakes', category='physical', scale=mapres, edgecolor='blue', facecolor='none')
ax.add_feature(lake_feature)

# Generate your layer data (assuming it's stored in 'layer_data' variable)
layer_data = np.random.rand(100, 100)

# Plot the layer using imshow
im = ax.imshow(layer_data, extent=[xmin, xmax, ymin, ymax], transform=projection, origin='lower', cmap='viridis')

plt.colorbar(im,location='bottom')

plt.show()

Produces this (note land_feature is commented out, and you can see lake outline if you squint):

[greglucas]

I think you'll want to look into geopandas: https://geopandas.org/en/stable/docs/user_guide/set_operations.html which will handle a lot of the intersections and geometry manipulations for you. Here is a quick example to remove the lakes from the land.

lake_df = gpd.GeoDataFrame({'geometry': lake_feature.geometries()})
land_df = gpd.GeoDataFrame({'geometry': land_feature.geometries()})
land_df = land_df.overlay(lake_df, how='symmetric_difference')
land_df.plot(ax=ax, facecolor='gray')

[SUPScientist]

That worked! Required a tiny bit of tweaking of the .overlay() method, but this is a nice solution. Thanks very much. I will say it's a shame that lakes/rivers are filled in using the land feature by default, but I really appreciate the rapid and helpful suggestions. Here's how I got it:

land_feature = cfeature.NaturalEarthFeature(name='land', category='physical', scale='50m', edgecolor='k', facecolor='gray')
# ax.add_feature(feature) # don't add feature if using this script for lakes!

# add lakes and repeat same procedure if needed for rivers
lake_feature = cfeature.NaturalEarthFeature('physical', 'lakes', '50m', edgecolor='blue', facecolor='none')
# ax.add_feature(lake_feature) # don't add feature if using this script for lakes!

# Convert to geodataframe for subtraction
lake_df = gpd.GeoDataFrame({'geometry': lake_feature.geometries()})
land_df = gpd.GeoDataFrame({'geometry': land_feature.geometries()})

# Subtract lakes from land so that we don't fill over the lakes
land_minus_lake_df = gpd.overlay(land_df, lake_df, how='symmetric_difference')

# Plot land_minus_lake_df
land_minus_lake_df.plot(ax=ax, facecolor='gray')

[dopplershift]

Unfortunately, there's nothing we can do about the land feature since that's outside our control, we just plot what we get from Natural Earth.

Sounds like there's nothing left for us to do here so closing, but feel free to re-open if you think there's more for Cartopy to handle here.

[SUPScientist]

there's nothing we can do about the land feature since that's outside our control, we just plot what we get from Natural Earth

Fair enough! Thanks for the pointers, and good to know.