Missing region border when specifying a projection; In addition: a (faster) alternative to using projections

[mherrmann3]

1. A bug: the region border is not projected

When using 'projection': cartopy.crs.Mercator() in plot_args to csep.utils.plots.plot_spatial_dataset (i.e., via csep.core.forecasts.GriddedForecast.plot, the region boundary disappears. This doesn't happen if not specifying this argument (i.e., when PlateCaree is used).

Suggested fix: add transform=cartopy.crs.PlateCarree() to the .plot() command here.

2. A (faster) alternative to using projections

Instead of using a fancy projections (i.e. everything else than PlateCaree), you may consider to simply scale the axes according to the currently shown region:

import cartopy.feature as cfeature

fig = plt.figure()

# Set projection
ax = plt.axes(projection=cartopy.crs.PlateCarree())

# Re-aspect plot (only for plain matplotlib plots, or when using PlateCarree)
LATKM = 110.574
lat_region = <approximate latitude at the center of the plot>
ax.set_aspect(1 / (LATKM / 111.320 * np.cos(np.deg2rad(lat_region)))) 

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE, ...)
ax.coastlines(resolution='50m', ...)
ax.add_feature(cfeature.BORDERS, ...)

ax.scatter(...)

(Note that you don't need to specify the transform=cartopy.crs.PlateCarree() argument in the plot/scatter/pcolor functions, but you can.)

I found this approach completely sufficient at country-level zoom levels (or closer). By setting the correct aspect ratio at the current latitude, you'll hardly notice a difference to projection-based plots at such zoom levels (unless you plot larger countries like Russia, etc. 😉). Plotting is also faster (~2-3x), as it avoids using more expensive transform objects (such as Mercator(), Orthographic(), you name it); note that this speed increase is independent of the number of plotted events / grid cells. (Btw: In case you want to benchmark the plotting speed, be aware that cartopy uses some sort of memoizing, so you need to reload the module every time.)

Suggested "fix": Maybe we can implement an intelligent switch (and/or an arg) to use plain PlateCaree + set_aspect. Essentially, this would need to happen in csep.utils.plots.plot_basemap.

[wsavran]

thank you for the very well thought out report. i opened a pull request #110 that fixes this issue. it also adds an optional argument called 'fast' to plot_args that performs the approximate Mercator projection you mentioned.

[wsavran]

# Import statements
import timeit
from importlib import reload

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

from csep import load_catalog_forecast
from csep.utils.datasets import ucerf3_ascii_format_landers_fname
from csep.core.regions import (
    california_relm_region, 
    magnitude_bins, 
    create_space_magnitude_region
)
from csep.utils.plots import plot_spatial_dataset

fore = load_catalog_forecast(
    ucerf3_ascii_format_landers_fname,
    region = create_space_magnitude_region(
        magnitudes=magnitude_bins(4.95, 8.05, 0.2),
        region=california_relm_region())
)
rates = fore.get_expected_rates()

Benchmark Cartopy transform: PlateCarree

test_statements = '''
reload(cartopy)
ax = plot_spatial_dataset(
    rates.spatial_counts(cartesian=True), 
    rates.region, 
    plot_args = {
        'figsize': (10, 10),
        'projection': cartopy.crs.PlateCarree()
    }
)
plt.close(ax.get_figure())
'''

r = timeit.repeat(globals=globals(),
    stmt = test_statements,
    repeat = 10,
    number = 1
)

print(f'mean: {np.array(r).mean():.2f} s, '
      f'min: {np.array(r).min():.2f} s, ',
      f'max: {np.array(r).max():.2f} s')

mean: 0.22 s, min: 0.21 s,  max: 0.22 s

Benchmark Cartopy transform: Mercator

test_statements = '''
reload(cartopy)
ax = plot_spatial_dataset(
    rates.spatial_counts(cartesian=True), 
    rates.region, 
    plot_args = {
        'figsize': (10, 10),
        'projection': cartopy.crs.Mercator()
    }
)
plt.close(ax.get_figure())
'''

r = timeit.repeat(globals=globals(),
    stmt = test_statements,
    repeat = 10,
    number = 1
)

print(f'mean: {np.array(r).mean():.2f} s, '
      f'min: {np.array(r).min():.2f} s, ',
      f'max: {np.array(r).max():.2f} s')

mean: 5.68 s, min: 5.50 s,  max: 5.94 s

Benchmark Fast transform

test_statements = '''
reload(cartopy)
ax = plot_spatial_dataset(
    rates.spatial_counts(cartesian=True), 
    rates.region, 
    plot_args = {
        'figsize': (10, 10),
        'projection': 'fast'
    }
)
plt.close(ax.get_figure())
'''

r = timeit.repeat(globals=globals(),
    stmt = test_statements,
    repeat = 10,
    number = 1
)
print(f'mean: {np.array(r).mean():.2f} s, '
      f'min: {np.array(r).min():.2f} s, ',
      f'max: {np.array(r).max():.2f} s')

mean: 0.23 s, min: 0.23 s,  max: 0.24 s

[mherrmann3]

Great work, Bill! It's impressive how costly those transform operations are.
I believe this faster plotting allows users to play more with plot_args and find their preferred ones. (btw: kudos for having so many!) 😀

[pabloitu]

I love the fast projection option. It is super useful for exploratory visualizations. And believe it is also good to leave the cartopy.projection options, if one needs to be accurate or use weird cut projections like in NZ

[mherrmann3]

Damn, I just noticed that the aspect ratio simplification was wrong, see #110 (comment)

It did exactly the opposite axis scaling. You possibly didn't notice that yet due to testing it only on rather global parts. Much better noticeable at more regional level.

Sorry for the confusion.

[mherrmann3]

bump: @wsavran: what to do about: #125 (comment)?

[wsavran]

@mherrmann3 sorry I missed this! the issue was closed and I didn't see your last comment. I apologize about that. You will have to open a new PR to propose the changes to the aspect ratio. Once a pull request is closed and merged a new pull request must be made for additional changes.