Merge pull request #10112 from gem/plotshakemap

ARISTOTLE: plot shakemap

openquake/calculators/getters.py

@@ -419,6 +419,11 @@ def get_ebrupture(dstore, rup_id):  # used in show rupture
     return get_ebr(rec, geom, trt)
 
 
+def get_rupture_from_dstore(dstore, rup_id=0):
+    ebr = get_ebrupture(dstore, rup_id)
+    return ebr.rupture
+
+
 # this is never called directly; get_rupture_getters is used instead
 class RuptureGetter(object):
     """

openquake/calculators/postproc/plots.py

@@ -1,27 +1,28 @@
 # -*- coding: utf-8 -*-
 # vim: tabstop=4 shiftwidth=4 softtabstop=4
-# 
+#
 # Copyright (C) 2024, GEM Foundation
-# 
+#
 # OpenQuake is free software: you can redistribute it and/or modify it
 # under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
-# 
+#
 # OpenQuake is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Affero General Public License for more details.
-# 
+#
 # You should have received a copy of the GNU Affero General Public License
 # along with OpenQuake.  If not, see <http://www.gnu.org/licenses/>.
 
+import io
 import os
+import base64
 import numpy
 from shapely.geometry import MultiPolygon
 from openquake.commonlib import readinput, datastore
 from openquake.hmtk.plotting.patch import PolygonPatch
-from openquake.calculators.getters import get_ebrupture
 
 
 def import_plt():
@@ -32,7 +33,7 @@ def import_plt():
     return plt
 
 
-def add_borders(ax, read_df=readinput.read_countries_df, buffer=0):
+def add_borders(ax, read_df=readinput.read_countries_df, buffer=0, alpha=0.1):
     plt = import_plt()
     polys = read_df(buffer)['geom']
     cm = plt.get_cmap('RdBu')
@@ -41,9 +42,27 @@ def add_borders(ax, read_df=readinput.read_countries_df, buffer=0):
         colour = cm(1. * idx / num_colours)
         if isinstance(poly, MultiPolygon):
             for onepoly in poly.geoms:
-                ax.add_patch(PolygonPatch(onepoly, fc=colour, alpha=0.1))
+                ax.add_patch(PolygonPatch(onepoly, fc=colour, alpha=alpha))
         else:
-            ax.add_patch(PolygonPatch(poly, fc=colour, alpha=0.1))
+            ax.add_patch(PolygonPatch(poly, fc=colour, alpha=alpha))
+    return ax
+
+
+def add_populated_places(ax, xlim, ylim, read_df=readinput.read_populated_places_df,
+                         lon_field='longitude', lat_field='latitude',
+                         label_field='name'):
+    data = read_df(lon_field, lat_field, label_field)
+    if data is None:
+        return ax
+    data = data[(data[lon_field] >= xlim[0]) & (data[lon_field] <= xlim[1])
+                & (data[lat_field] >= ylim[0]) & (data[lat_field] <= ylim[1])]
+    if len(data) == 0:
+        return ax
+    ax.scatter(data[lon_field], data[lat_field], label="Populated places",
+               s=2, color='black', alpha=0.5)
+    for _, row in data.iterrows():
+        ax.text(row[lon_field], row[lat_field], row[label_field], fontsize=7,
+                ha='right', alpha=0.5)
     return ax
 
 
@@ -59,6 +78,64 @@ def get_country_iso_codes(calc_id, assetcol):
     return id_0_str
 
 
+def plt_to_base64(plt):
+    """
+    The base64 string can be passed to a Django template and embedded
+    directly in HTML, without having to save the image to disk
+    """
+    bio = io.BytesIO()
+    plt.savefig(bio, format='png', bbox_inches='tight')
+    bio.seek(0)
+    img_base64 = base64.b64encode(bio.getvalue()).decode('utf-8')
+    return img_base64
+
+
+def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10),
+                  with_populated_places=False, return_base64=False,
+                  rupture=None):
+    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.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()
+    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,
+            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)
+    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:
+        return plt
+
+
 def plot_avg_gmf(ex, imt):
     plt = import_plt()
     _fig, ax = plt.subplots(figsize=(10, 10))
@@ -94,49 +171,68 @@ def plot_avg_gmf(ex, imt):
     return plt
 
 
-def add_surface(ax, surface, label):
-    ax.fill(*surface.get_surface_boundaries(), alpha=.5, edgecolor='grey',
-            label=label)
+def add_surface(ax, surface, label, alpha=0.5, facecolor=None, linestyle='-'):
+    fill_params = {
+        'alpha': alpha,
+        'edgecolor': 'grey',
+        'label': label
+    }
+    if facecolor is not None:
+        fill_params['facecolor'] = facecolor
+    ax.fill(*surface.get_surface_boundaries(), **fill_params)
     return surface.get_bounding_box()
 
 
-def add_rupture(ax, dstore, rup_id=0):
-    ebr = get_ebrupture(dstore, rup_id)
-    rup = ebr.rupture
+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'):
         min_x = 180
         max_x = -180
         min_y = 90
         max_y = -90
         for surf_idx, surface in enumerate(rup.surface.surfaces):
             min_x_, max_x_, max_y_, min_y_ = add_surface(
-                ax, surface, 'Surface %d' % surf_idx)
+                ax, surface, 'Surface %d' % surf_idx, alpha=surf_alpha,
+                facecolor=surf_facecolor, linestyle=surf_linestyle)
             min_x = min(min_x, min_x_)
             max_x = max(max_x, max_x_)
             min_y = min(min_y, min_y_)
             max_y = max(max_y, max_y_)
     else:
-        min_x, max_x, max_y, min_y = add_surface(ax, rup.surface, 'Surface')
+        min_x, max_x, max_y, min_y = add_surface(
+            ax, rup.surface, 'Surface', alpha=surf_alpha, facecolor=surf_facecolor,
+            linestyle=surf_linestyle)
     ax.plot(rup.hypocenter.x, rup.hypocenter.y, marker='*',
-            color='orange', label='Hypocenter', alpha=.5,
+            color='orange', label='Hypocenter', alpha=hypo_alpha,
             linestyle='', markersize=8)
     return ax, min_x, min_y, max_x, max_y
 
 
-def plot_rupture(dstore):
+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
     plt = import_plt()
-    _fig, ax = plt.subplots(figsize=(10, 10))
+    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)
-    # assuming there is only 1 rupture, so rup_id=0
-    ax, min_x, min_y, max_x, max_y = add_rupture(ax, dstore, rup_id=0)
+    ax, min_x, min_y, max_x, max_y = add_rupture(ax, rup)
     ax = add_borders(ax)
-    BUF_ANGLE = 4
-    ax.set_xlim(min_x - BUF_ANGLE, max_x + BUF_ANGLE)
-    ax.set_ylim(min_y - BUF_ANGLE, max_y + BUF_ANGLE)
+    BUF_ANGLE = 1
+    xlim = (min_x - BUF_ANGLE, max_x + BUF_ANGLE)
+    ylim = (min_y - BUF_ANGLE, max_y + BUF_ANGLE)
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
+    if with_populated_places:
+        ax = add_populated_places(ax, xlim, ylim)
     ax.legend()
-    return plt
+    if return_base64:
+        return plt_to_base64(plt)
+    else:
+        return plt
 
 
 def add_surface_3d(ax, surface, label):
@@ -147,13 +243,11 @@ def add_surface_3d(ax, surface, label):
     ax.plot_surface(lon_grid, lat_grid, depth_grid, alpha=0.5, label=label)
 
 
-def plot_rupture_3d(dstore):
+def plot_rupture_3d(rup):
     # NB: matplotlib is imported inside since it is a costly import
     plt = import_plt()
     fig = plt.figure()
     ax = fig.add_subplot(111, projection='3d')
-    ebr = get_ebrupture(dstore, rup_id=0)
-    rup = ebr.rupture
     if hasattr(rup.surface, 'surfaces'):
         for surf_idx, surface in enumerate(rup.surface.surfaces):
             add_surface_3d(ax, surface, 'Surface %d' % surf_idx)

openquake/commands/plot.py

@@ -29,6 +29,7 @@
 from openquake.hazardlib.geo.utils import PolygonPlotter
 from openquake.hazardlib.contexts import Effect, get_effect_by_mag
 from openquake.hazardlib.calc.filters import getdefault, IntegrationDistance
+from openquake.calculators.getters import get_rupture_from_dstore
 from openquake.calculators.extract import (
     Extractor, WebExtractor, clusterize)
 from openquake.calculators.postproc.plots import (
@@ -58,6 +59,7 @@ def make_figure_magdist(extractors, what):
         ax.legend()
     return plt
 
+
 def make_figure_hcurves(extractors, what):
     """
     $ oq plot "hcurves?kind=mean&imt=PGA&site_id=0"
@@ -1042,7 +1044,8 @@ def make_figure_rupture(extractors, what):
     """
     [ex] = extractors
     dstore = ex.dstore
-    return plot_rupture(dstore)
+    rup = get_rupture_from_dstore(dstore, rup_id=0)
+    return plot_rupture(rup)
 
 
 def make_figure_rupture_3d(extractors, what):
@@ -1051,7 +1054,8 @@ def make_figure_rupture_3d(extractors, what):
     """
     [ex] = extractors
     dstore = ex.dstore
-    return plot_rupture_3d(dstore)
+    rup = get_rupture_from_dstore(dstore, rup_id=0)
+    return plot_rupture_3d(rup)
 
 
 def plot_wkt(wkt_string):

openquake/commands/plot_assets.py

@@ -22,6 +22,7 @@
 import logging
 from openquake.commonlib import datastore
 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
@@ -83,8 +84,8 @@ def main(calc_id: int = -1, site_model=False,
         print('rupture(%s, %s), dist=%s' % (lon, lat, dist))
         if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE':
             # assuming there is only 1 rupture, so rup_id=0
-            ax, _min_x, _min_y, _max_x, _max_y = add_rupture(
-                ax, dstore, 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)
         else:
             p.scatter(xlon, xlat, marker='*', color='orange',
                       label='hypocenter', alpha=.5)

openquake/commands/tests/independence_test.py

@@ -35,6 +35,9 @@
 
 
 class IndependenceTestCase(unittest.TestCase):
+    def test_hazardlib(self):
+        assert_independent('openquake.hazardlib', 'openquake.calculators')
+
     def test_risklib(self):
         assert_independent('openquake.risklib', 'openquake.commonlib')
         assert_independent('openquake.risklib', 'openquake.calculators')

openquake/commonlib/readinput.py

@@ -1662,6 +1662,20 @@ def read_geometries(fname, code, buffer=0):
     return pandas.DataFrame(dict(code=codes, geom=geoms))
 
 
+@functools.lru_cache()
+def read_populated_places(fname, lon_name, lat_name, label_name):
+    """
+    Reading coordinates and names of populated places from a CSV file
+
+    :returns: a Pandas DataFrame
+    """
+    data = pandas.read_csv(fname)
+    expected_colnames_set = {lon_name, lat_name, label_name}
+    if not expected_colnames_set.issubset(data.columns):
+        raise ValueError(f"CSV file must contain {expected_colnames_set} columns.")
+    return data
+
+
 def read_mosaic_df(buffer):
     """
     :returns: a DataFrame of geometries for the mosaic models
@@ -1680,6 +1694,21 @@ def read_countries_df(buffer=0.1):
     return read_geometries(fname, 'shapeGroup', buffer)
 
 
+def read_populated_places_df(lon_field='longitude', lat_field='latitude',
+                             label_field='name'):
+    """
+    Reading from a 'worldcities.csv' file in the mosaic_dir, if present, or returning
+    None otherwise
+
+    :returns: a DataFrame of coordinates and names of populated places
+    """
+    mosaic_dir = config.directory.mosaic_dir
+    fname = os.path.join(mosaic_dir, 'worldcities.csv')
+    if not os.path.isfile(fname):
+        return
+    return read_populated_places(fname, lon_field, lat_field, label_field)
+
+
 def read_source_models(fnames, hdf5path='', **converterparams):
     """
     :param fnames: a list of source model files

openquake/engine/aristotle.py

@@ -121,6 +121,8 @@ def get_aristotle_params(arist):
     inputs = {'exposure': [arist.exposure_hdf5],
               'job_ini': '<in-memory>'}
     rupdic = get_rupture_dict(arist.rupture_dict, arist.ignore_shakemap)
+    if 'shakemap_array' in rupdic:
+        del rupdic['shakemap_array']
     if arist.station_data_file is None:
         # NOTE: giving precedence to the station_data_file uploaded via form
         try: