GSHHS coastlines

[emlynmjones]

The GSHHS coastlines appear to be performing poorly for some regions. There is a 5km offset between the GSHHS coastline and those from OpenStreet Map. When plotted with a google earth backdrop, the OpenStreetmap coastlines appear correct. See example below (RED=GSHHS, GREEN = OpenStreetmap):

In the Python environment, it looks like the "Natural Earth" shape files faithfully represent the coastlines.

I suggest changing from GSHHS to Natural Earth as the default coastlines used by Cartopy in EMSARRAY.

Cheers, EJ

[sharon-tickell]

That sort of offset looks like the GSHHS layer may use a different CRS projection?

According to google-fu, Google Earth and the OpenStreetMap database are EPSG:4326 (units in degrees), Google Maps and the OpenStreetMap tiles and WMS are EPSG:3857 (spherical mercator a.k.a. web mercator, units in metres).

What is the GSHHS coastline in?

[emlynmjones]

I think there may be something odd with GSHHS, if we look a little to the south, there is an island that is correctly represented, but other island not:

[mx-moth]

I've mostly tested plots around Tasmania, where the GSHHS coastlines are adequate, although not great to be honest. With your example, where some islands are aligned and some are not, I suspect the issue is entirely within the GSHHS coastlines dataset.

Fixing GSHHS is out of scope for emsarray itself, but if you have a suggestion for a better coastline source I'd be happy to switch to that instead of GSHHS

[sharon-tickell]

I suspect GSHHS may just need reprojecting, but I usually defer to Erin for expertise there!

If you need them, the OzCoasts GeoServer has the GeoNetwork Natural Earth 50m Coastline (gn:ne_50m_coastline), The Bright Earth World eAtlas Coastline (ea-be:World_e-Atlas_VMap0-AU-Coast100k_Coast-split) and the GBR Features, Reefs, Coastline layer published by GBRMPA (ea-be:GBR_GBRMPA_GBR-features).

The CMAR GeoServer also has the Bright Earth eAtlas Coastline and the Reef Features layer, and adds the a caab:coastlines layer.

[emlynmjones]

We've had issues with GSHHS elsewhere. @sharon-tickell - you're likely right, in that there could be a projection issue, but we've noticed GSHHS to have some undesirable errors around the Aus coastline.

My preference would be to use the Natural Earth coastlines if possible. To date, they have given good results.

[mx-moth]

I added the Natural Earth coastlines to a plot of the Hobart region, and the resolution of the features is fairly poor at this scale. This is using the highest resolution 10m scale Natural Earth dataset. Which is worse: low resolution coastlines with Natural Earth, or some islands in the GBR being misaligned?

This is the Natural Earth coastlines for the region above, and below is the GSHHS coastlines for the same region

[mx-moth]

If you are making custom plots, you can add the natural earth coastlines instead of the GSHHS coastlines by replacing the call to emsarray.plot.add_coastlines(axes) with something like:

import cartopy.feature

axes.add_feature(cartopy.feature.COASTLINE)
axes.add_feature(cartopy.feature.LAND)

A larger example is available in the cartopy docs

[mx-moth]

The errors noticed above seem to be localised to specific areas, while most of the GSHHS coastlines are correctly aligned. This issue should be raised against GSHHS.

The Natural Earth coastlines are not a suitable replacement as they are too low resolution for many applications.

emsarray now allows users to disable adding coastlines during plotting if the GSHHS coastlines are not suitable for a specific location. Users can then add their own choice of coast line to the plot. For example the following will render a plot using Natural Earth coastlines:

import emsarray
import matplotlib.pyplot as plt
from cartopy.feature import NaturalEarthFeature

figure = plt.figure()
bran2020 = emsarray.tutorial.open_dataset('bran2020')
bran2020.ems.plot_on_figure(
    figure, bran2020['temp'].isel(Time=0, st_ocean=0),
    coast=False,
)
figure.axes[0].add_feature(
    NaturalEarthFeature('physical', 'coastline', '50m'),
    facecolor=(0.7, 0.7, 0.7, 0.5),
    edgecolor='darkgrey',
    linewidth=0.5,
)

[mx-moth]

A further example of adding custom coastlines from a Shapefile

import emsarray
import matplotlib.pyplot as plt
import shapefile
import shapely.geometry

from pyproj import CRS, Transformer
from shapely.ops import transform
from zipfile import ZipFile

# Open the dataset
bran2020 = emsarray.tutorial.open_dataset('bran2020')

# Open the coastline shapefile
# This coastfile can be downloaded from
# https://www.thelist.tas.gov.au/app/content/data/geo-meta-data-record?detailRecordUID=c10ccd03-dd2d-492d-8c2f-fca25baed194
shp_path = "LIST_COASTLINE_STATEWIDE.zip"
shp = shapefile.Reader(shp_path)
shp_zip = ZipFile(open(shp_path, 'rb'))
shp_crs = CRS(shp_zip.read(next(name for name in shp_zip.namelist() if name.endswith('.prj'))).decode('utf-8'))

# Extract the coastline as a shapely MultiLineString
shp_lines = shp.shapes()
coastline = shapely.geometry.MultiLineString([
    shapely.geometry.shape(s) for s in shp_lines
])

# Reproject in to the dataset CRS
transformer = Transformer.from_crs(shp_crs, bran2020.ems.data_crs)
coastline = transform(transformer.transform, coastline)

# Plot the data
figure = plt.figure()
bran2020.ems.plot_on_figure(
    figure, bran2020['temp'].isel(Time=0, st_ocean=0),
    coast=False,
)

# Plot the coastline
axes = figure.axes[0]
for line in coastline.geoms:
    axes.plot(*line.xy, '-k', linewidth=0.5)