Bathymetry contour using NaturalEarthFeature

[philippemiron]

Hi all,

I'm trying to plot bathymetry contour line and running into some issues. Depending on the selected region, I get straight lines that looks like artifacts resulting from the combination of multiple shape files.

Anyone has an idea how to fix this ?

You can see in this simple example.

import numpy as np
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

fig = plt.figure(figsize=(10, 4), dpi=300)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(), aspect='equal')
ax.set_extent([-100, 10, 0, 90], crs=ccrs.PlateCarree())

# ticks
ax.set_xticks(np.arange(-90, 30, 30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(40, 100, 20), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# add land and coastline and 1000m bathymetry line
ax.add_feature(cfeature.LAND, facecolor='grey', zorder=1)
ax.add_feature(cfeature.COASTLINE, linewidth=0.25, zorder=1)
bathym = cfeature.NaturalEarthFeature(name='bathymetry_J_1000', scale='10m', category='physical')
ax.add_feature(bathym, facecolor='none', edgecolor='black', linestyle='dashed', linewidth=1)

[greglucas]

If they are different geometries within the shapefile you may need to do a union on the two shapes to combine them into a single one. Shapely or geopandas can take care of that for you, see: https://shapely.readthedocs.io/en/latest/manual.html#object.union

[philippemiron]

Thanks for the suggestion.

I tried the union of the bathym object in the previous code:

bathym = cfeature.NaturalEarthFeature(name='bathymetry_J_1000', scale='10m', category='physical')

a = None
for geo in bathym.geometries():
    if a == None:
        a = geo
    else:
        a = a.union(geo)

which now return a shapely.geometry.multipolygon.MultiPolygon.

But then I get another error trying to add the new polygon to the figure.

ax.add_feature(a)
AttributeError: 'MultiPolygon' object has no attribute 'kwargs'

[greglucas]

You don't have a "feature" anymore, you have a bunch of geometries. Give this a try: ax.add_geometries(a, facecolor='none', edgecolor='black', linestyle='dashed', linewidth=1)

[philippemiron]

Thanks a lot.

This is actually a cleaner way to perform the union:

from shapely.ops import cascaded_union

bathym = cfeature.NaturalEarthFeature(name='bathymetry_J_1000', scale='10m', category='physical')
bathym = cascaded_union(list(bathym.geometries()))
ax.add_geometries(bathym, facecolor='none', edgecolor='black', linestyle='dashed', linewidth=1, crs=ccrs.PlateCarree())

Not sure if this is something that could be added to Cartopy?

[gabrilu]

Hi,

Thank you for the contribution. I came across this same problem today.

For those working with south ocean areas, performing a union may still preserve those polygon lines in some cases. Or at least preserve tiny single artifacts resulting from polygon borders not fitting perfectly one next to another (as for my case). Here is a quick increment for those who may encounter the same problem:

from shapely.geometry import JOIN_STYLE eps = 0.0001 #Increase eps depending on your contour resolution
bathym = bathym.buffer(eps, 1, join_style=JOIN_STYLE.mitre).buffer(-eps, 1, join_style=JOIN_STYLE.mitre)

From this topic