Map scalebar has incorrect length

[tsonne]

The scalebars, e.g. in Shakemap, are shown at incorrect lengths on the map if the location is not at the Equator. For a simple scenario in Alaska at latitude 66 deg, the scalebar shows 150 km when it really is about 60 km long based on the map extent and grid lines. For an event in Equador at latitude 0 deg, the scalebar has the correct length.

This seems to affect also the example images in the Shakemap Tutorial etc.

To reproduce it, I created some test Shakemaps with simple event.xml input files:

<earthquake id="scenario1" netid="TEST" network="TEST" lat="66.0" lon="-161.5" depth="4"
 mag="4.5" time="2020-08-27T12:00:00.000000Z" locstring="Alaska"/>

<earthquake id="scenario1" netid="TEST" network="TEST" lat="0.0" lon="-80.0" depth="4"
 mag="4.5" time="2020-08-27T12:00:00.000000Z" locstring="Equador"/>

See the images below. I've also drawn thick red lines to indicate the actually expected length for the given scalebar using GIMP.

Versions

• OS: Ubuntu 18.04
• Python 3.7
• Impact Utils 0.8.23
• Numpy 1.18.5
• Cartopy 0.18.0

[tsonne]

I found a workaround to fix the scalebar, after line 181 in impactutils/mapping/scalebar.py insert this code:

    # Get factor to convert data units to meters in a different way:
    # get axis limits in lon, lat
    lonlat_ext = ax.get_extent(cartopy.crs.PlateCarree())
    # get axis limits in meters (TransverseMercator based on map center)
    m_ext = ax.get_extent(cartopy.crs.TransverseMercator(
                central_longitude=((lonlat_ext[0] + lonlat_ext[1]) / 2.),
                central_latitude=((lonlat_ext[2] + lonlat_ext[3]) / 2.),
                approx=False))
    # get ratio meters / data units and inverse
    dataunits_to_meters = (m_ext[1] - m_ext[0])/(xmax - xmin)
    meters_to_dataunits = 1./dataunits_to_meters

Now the length is correct at all latitudes, see the updated Alaska example:

Incidentally the extra town comes up after also fixing the cities issue #135

A possible alternative to fix the coordinates was proposed by a coworker. For line 171 in impactutils/mapping/scalebar.py:

            #dataunits_to_meters = UNITS[proj_units]
            lonlat_ext = ax.get_extent(cartopy.crs.PlateCarree())
            clat = (lonlat_ext[2] + lonlat_ext[3]) / 2.
            dataunits_to_meters = UNITS[proj_units] * \
                (np.pi*6378137*np.cos(np.deg2rad(clat)) / \
                (180*np.sqrt(1-((6378137**2-6356752**2) / \
                6378137**2)*np.sin(np.deg2rad(clat))**2)))/(111.32*10**3)

It gets the ratio of length of degree of longitude at the Equator versus the map's central latitude. Somehow, the two solutions are slightly different though (~3%), we're not sure why.

[emthompson-usgs]

@mhearne-usgs Have you looked at this?

[mhearne-usgs]

We could look at the solutions outlined here as well:

SciTools/cartopy#490

[emthompson-usgs]

Thanks @tsonne. Sorry this took a very long time for us to address. I added your solution in PR #140.

[emthompson-usgs]

@tsonne, also you seem to be adept at both finding bugs and also fixes. In the future, we would welcome pull requests with your fixes.