ARISTOTLE: Plot shakemaps, ruptures and assets on top of basemaps obtained from external tiles services

openquake/calculators/base.py

@@ -37,7 +37,7 @@
 import configparser
 import collections
 
-from openquake.commands.plot_assets import main as plot_assets
+from openquake.commands.plot_assets_webmercator import main as plot_assets
 from openquake.baselib import general, hdf5, config
 from openquake.baselib import parallel
 from openquake.baselib.performance import Monitor
@@ -658,7 +658,7 @@ def import_perils(self):  # called in pre_execute
                 logging.warning('No sites were affected by %s' % name)
             data[name] = peril
             names.append(name)
-        self.datastore['events'] = numpy.zeros(1, rupture.events_dt)     
+        self.datastore['events'] = numpy.zeros(1, rupture.events_dt)
         create_gmf_data(self.datastore, [], names, data, N)
 
     def pre_execute(self):

openquake/calculators/postproc/plots.py

@@ -20,6 +20,8 @@
 import os
 import base64
 import numpy
+import contextily as ctx
+from pyproj import Transformer
 from shapely.geometry import MultiPolygon
 from openquake.commonlib import readinput, datastore
 from openquake.hmtk.plotting.patch import PolygonPatch
@@ -33,6 +35,36 @@ def import_plt():
     return plt
 
 
+def add_basemap(ax, min_x, min_y, max_x, max_y, source=ctx.providers.CartoDB.Positron):
+    # NOTE: another interesting option:
+    # source = ctx.providers.TopPlusOpen.Grey
+    img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source)
+    ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1)
+    ax.text(
+        0.01, 0.01,  # Position: Bottom-left corner (normalized coordinates)
+        source['attribution'],
+        transform=ax.transAxes,  # Place text relative to axes
+        fontsize=8, color="black", alpha=0.5,
+        ha="left",  # Horizontal alignment
+        va="bottom",  # Vertical alignment
+    )
+    return ax
+
+
+def adjust_limits(x_min, x_max, y_min, y_max, padding=0.5):
+    # Make the plot display all items with some margin, looking square
+    x_min, x_max = x_min - padding, x_max + padding
+    y_min, y_max = y_min - padding, y_max + padding
+    x_range = x_max - x_min
+    y_range = y_max - y_min
+    max_range = max(x_range, y_range)
+    x_center = (x_min + x_max) / 2
+    y_center = (y_min + y_max) / 2
+    xlim = x_center - max_range / 2, x_center + max_range / 2
+    ylim = y_center - max_range / 2, y_center + max_range / 2
+    return xlim, ylim
+
+
 def add_borders(ax, read_df=readinput.read_countries_df, buffer=0, alpha=0.1):
     plt = import_plt()
     polys = read_df(buffer)['geom']
@@ -100,36 +132,37 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10),
         matplotlib.use(backend)
     _fig, ax = plt.subplots(figsize=figsize)
     ax.set_aspect('equal')
-    ax.grid(True)
+    # ax.grid(True)
     ax.set_xlabel('Longitude')
     ax.set_ylabel('Latitude')
     title = 'Avg GMF for %s' % imt
     ax.set_title(title)
     gmf = shakemap_array['val'][imt]
-    markersize = 5
-    coll = ax.scatter(shakemap_array['lon'], shakemap_array['lat'], c=gmf,
-                      cmap='jet', s=markersize)
-    plt.colorbar(coll)
-    ax = add_borders(ax, alpha=0.2)
-    BUF_ANGLE = 1
-    min_x = shakemap_array['lon'].min()
-    max_x = shakemap_array['lon'].max()
-    min_y = shakemap_array['lat'].min()
-    max_y = shakemap_array['lat'].max()
+    markersize = 0.005
+    transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True)
+    x_webmercator, y_webmercator = transformer.transform(
+        shakemap_array['lon'], shakemap_array['lat'])
+    min_x, min_y, max_x, max_y = (min(x_webmercator),
+                                  min(y_webmercator),
+                                  max(x_webmercator),
+                                  max(y_webmercator))
     if rupture is not None:
-        ax, rup_min_x, rup_min_y, rup_max_x, rup_max_y = add_rupture(
-            ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.9,
+        ax, rup_min_x, rup_min_y, rup_max_x, rup_max_y = add_rupture_webmercator(
+            ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=1,
             surf_facecolor='none', surf_linestyle='--')
         min_x = min(min_x, rup_min_x)
         max_x = max(max_x, rup_max_x)
         min_y = min(min_y, rup_min_y)
         max_y = max(max_y, rup_max_y)
-    xlim = (min_x - BUF_ANGLE, max_x + BUF_ANGLE)
-    ylim = (min_y - BUF_ANGLE, max_y + BUF_ANGLE)
+    xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5)
+    min_x, max_x = xlim
+    min_y, max_y = ylim
+    add_basemap(ax, min_x, min_y, max_x, max_y)
+    coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize,
+                      alpha=0.4)
+    plt.colorbar(coll, ax=ax)
     ax.set_xlim(*xlim)
     ax.set_ylim(*ylim)
-    if with_populated_places:
-        ax = add_populated_places(ax, xlim, ylim)
     if return_base64:
         return plt_to_base64(plt)
     else:
@@ -140,34 +173,34 @@ def plot_avg_gmf(ex, imt):
     plt = import_plt()
     _fig, ax = plt.subplots(figsize=(10, 10))
     ax.set_aspect('equal')
-    ax.grid(True)
+    # ax.grid(True)
     ax.set_xlabel('Lon')
     ax.set_ylabel('Lat')
-
     title = 'Avg GMF for %s' % imt
     assetcol = get_assetcol(ex.calc_id)
     if assetcol is not None:
         country_iso_codes = get_country_iso_codes(ex.calc_id, assetcol)
         if country_iso_codes is not None:
             title += ' (Countries: %s)' % country_iso_codes
     ax.set_title(title)
-
     avg_gmf = ex.get('avg_gmf?imt=%s' % imt)
     gmf = avg_gmf[imt]
     markersize = 5
-    coll = ax.scatter(avg_gmf['lons'], avg_gmf['lats'], c=gmf, cmap='jet',
-                      s=markersize)
-    plt.colorbar(coll)
-
-    ax = add_borders(ax)
-
-    minx = avg_gmf['lons'].min()
-    maxx = avg_gmf['lons'].max()
-    miny = avg_gmf['lats'].min()
-    maxy = avg_gmf['lats'].max()
-    w, h = maxx - minx, maxy - miny
-    ax.set_xlim(minx - 0.2 * w, maxx + 0.2 * w)
-    ax.set_ylim(miny - 0.2 * h, maxy + 0.2 * h)
+    transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True)
+    x_webmercator, y_webmercator = transformer.transform(
+        avg_gmf['lons'], avg_gmf['lats'])
+    min_x, min_y, max_x, max_y = (min(x_webmercator),
+                                  min(y_webmercator),
+                                  max(x_webmercator),
+                                  max(y_webmercator))
+    xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5)
+    min_x, max_x = xlim
+    min_y, max_y = ylim
+    add_basemap(ax, min_x, min_y, max_x, max_y)
+    coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize)
+    plt.colorbar(coll, ax=ax)
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
     return plt
 
 
@@ -183,6 +216,25 @@ def add_surface(ax, surface, label, alpha=0.5, facecolor=None, linestyle='-'):
     return surface.get_bounding_box()
 
 
+def add_surface_webmercator(
+        ax, surface, label, alpha=0.5, facecolor=None, linestyle='-'):
+    transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
+    lon, lat = surface.get_surface_boundaries()
+    x, y = transformer.transform(lon, lat)  # Transform lat/lon to Web Mercator
+    fill_params = {
+        'alpha': alpha,
+        'edgecolor': 'grey',
+        'label': label
+    }
+    if facecolor is not None:
+        fill_params['facecolor'] = facecolor
+    ax.fill(x, y, **fill_params)
+    lon_min, lon_max, lat_max, lat_min = surface.get_bounding_box()
+    x_min, y_max = transformer.transform(lon_min, lat_max)  # Top-left corner
+    x_max, y_min = transformer.transform(lon_max, lat_min)  # Bottom-right corner
+    return x_min, x_max, y_min, y_max
+
+
 def add_rupture(ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5,
                 surf_facecolor=None, surf_linestyle='-'):
     if hasattr(rup.surface, 'surfaces'):
@@ -208,6 +260,30 @@ def add_rupture(ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5,
     return ax, min_x, min_y, max_x, max_y
 
 
+def add_rupture_webmercator(
+        ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5,
+        surf_facecolor=None, surf_linestyle='-'):
+    transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
+    min_x, max_x = float('inf'), float('-inf')
+    min_y, max_y = float('inf'), float('-inf')
+    if hasattr(rup.surface, 'surfaces'):
+        for surf_idx, surface in enumerate(rup.surface.surfaces):
+            min_x_, max_x_, min_y_, max_y_ = add_surface_webmercator(
+                ax, surface, 'Surface %d' % surf_idx, alpha=surf_alpha,
+                facecolor=surf_facecolor, linestyle=surf_linestyle)
+            min_x, max_x = min(min_x, min_x_), max(max_x, max_x_)
+            min_y, max_y = min(min_y, min_y_), max(max_y, max_y_)
+    else:
+        min_x, max_x, min_y, max_y = add_surface_webmercator(
+            ax, rup.surface, 'Surface', alpha=surf_alpha,
+            facecolor=surf_facecolor, linestyle=surf_linestyle)
+    hypo_x, hypo_y = transformer.transform(rup.hypocenter.longitude,
+                                           rup.hypocenter.latitude)
+    ax.plot(hypo_x, hypo_y, marker='*', color='orange', label='Hypocenter',
+            alpha=hypo_alpha, linestyle='', markersize=hypo_markersize)
+    return ax, min_x, min_y, max_x, max_y
+
+
 def plot_rupture(rup, backend=None, figsize=(10, 10),
                  with_populated_places=False, return_base64=False):
     # NB: matplotlib is imported inside since it is a costly import
@@ -221,9 +297,7 @@ def plot_rupture(rup, backend=None, figsize=(10, 10),
     ax.grid(True)
     ax, min_x, min_y, max_x, max_y = add_rupture(ax, rup)
     ax = add_borders(ax)
-    BUF_ANGLE = 1
-    xlim = (min_x - BUF_ANGLE, max_x + BUF_ANGLE)
-    ylim = (min_y - BUF_ANGLE, max_y + BUF_ANGLE)
+    xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y)
     ax.set_xlim(*xlim)
     ax.set_ylim(*ylim)
     if with_populated_places:
@@ -235,6 +309,32 @@ def plot_rupture(rup, backend=None, figsize=(10, 10),
         return plt
 
 
+def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), return_base64=False):
+    # NB: matplotlib is imported inside since it is a costly import
+    plt = import_plt()
+    if backend is not None:
+        # we may need to use a non-interactive backend
+        import matplotlib
+        matplotlib.use(backend)
+    _fig, ax = plt.subplots(figsize=figsize)
+    ax.set_aspect('equal')
+    # ax.grid(True)
+    ax, min_x, min_y, max_x, max_y = add_rupture_webmercator(
+        ax, rup, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.3,
+        surf_linestyle='--')
+    xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5)
+    min_x, max_x = xlim
+    min_y, max_y = ylim
+    add_basemap(ax, min_x, min_y, max_x, max_y, source=ctx.providers.TopPlusOpen.Color)
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
+    ax.legend()
+    if return_base64:
+        return plt_to_base64(plt)
+    else:
+        return plt
+
+
 def add_surface_3d(ax, surface, label):
     lon, lat, depth = surface.get_surface_boundaries_3d()
     lon_grid = numpy.array([[lon[0], lon[1]], [lon[3], lon[2]]])

openquake/commands/plot.py

@@ -33,7 +33,8 @@
 from openquake.calculators.extract import (
     Extractor, WebExtractor, clusterize)
 from openquake.calculators.postproc.plots import (
-    plot_avg_gmf, import_plt, add_borders, plot_rupture, plot_rupture_3d)
+    plot_avg_gmf, import_plt, add_borders, plot_rupture, plot_rupture_webmercator,
+    plot_rupture_3d)
 from openquake.calculators.postproc.aelo_plots import (
     plot_mean_hcurves_rtgm, plot_disagg_by_src, plot_governing_mce)
 
@@ -1048,6 +1049,16 @@ def make_figure_rupture(extractors, what):
     return plot_rupture(rup)
 
 
+def make_figure_rupture_webmercator(extractors, what):
+    """
+    $ oq plot "rupture_webmercator?"
+    """
+    [ex] = extractors
+    dstore = ex.dstore
+    rup = get_rupture_from_dstore(dstore, rup_id=0)
+    return plot_rupture_webmercator(rup)
+
+
 def make_figure_rupture_3d(extractors, what):
     """
     $ oq plot "rupture_3d?"

openquake/commands/plot_assets.py

@@ -24,8 +24,7 @@
 from openquake.hazardlib.geo.utils import cross_idl, get_bbox
 from openquake.calculators.getters import get_rupture_from_dstore
 from openquake.calculators.postproc.plots import (
-    add_borders, get_assetcol, get_country_iso_codes)
-from openquake.calculators.postproc.plots import add_rupture
+    add_borders, get_assetcol, get_country_iso_codes, add_rupture, adjust_limits)
 
 
 def main(calc_id: int = -1, site_model=False,
@@ -85,7 +84,7 @@ def main(calc_id: int = -1, site_model=False,
         if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE':
             # assuming there is only 1 rupture, so rup_id=0
             rup = get_rupture_from_dstore(dstore, rup_id=0)
-            ax, _min_x, _min_y, _max_x, _max_y = add_rupture(ax, rup)
+            ax, _minx, _miny, _maxx, _maxy = add_rupture(ax, rup)
         else:
             p.scatter(xlon, xlat, marker='*', color='orange',
                       label='hypocenter', alpha=.5)
@@ -99,9 +98,13 @@ def main(calc_id: int = -1, site_model=False,
     else:
         minx, miny, maxx, maxy = get_bbox(
             assetcol['lon'], assetcol['lat'], xlon, xlat)
-    BUF_ANGLE = 1
-    ax.set_xlim(minx - BUF_ANGLE, maxx + BUF_ANGLE)
-    ax.set_ylim(miny - BUF_ANGLE, maxy + BUF_ANGLE)
+    minx = min(minx, _minx)
+    maxx = max(maxx, _maxx)
+    miny = min(miny, _miny)
+    maxy = max(maxy, _maxy)
+    xlim, ylim = adjust_limits(minx, maxx, miny, maxy)
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
 
     country_iso_codes = get_country_iso_codes(calc_id, assetcol)
     if country_iso_codes is not None: