Adds logging of shakemap statistis in shakemap_from_quakeml workflow

shakyground2/workflows.py

@@ -2,6 +2,8 @@
 Defines complete end-to-end shakemap calculation workflows
 """
 import os
+import numpy as np
+from scipy import stats
 from typing import Optional, Dict, List, Tuple, Union
 from shakyground2.earthquake import Earthquake
 from shakyground2.site_model import SiteModel
@@ -46,6 +48,7 @@ class ShakemapWorkflowResult(object):
         "ground_motion_models",
         "tectonic_region",
         "earthquake",
+        "statistics",
     ]
 
     def __getitem__(self, key):
@@ -140,6 +143,7 @@ def shakemaps_from_quakeml(
     # Export to file (if an export directory is given) or to a dictionary of byte arrays
     results.shakemaps = {"mean": {}, "stddevs": {}}
     results.contours = {"mean": {}, "stddevs": {}}
+    results.statistics = {"mean": {}}
     for imt in imts:
         results.shakemaps["mean"][imt] = shakemap.to_geotiff(mean_shakemaps, imt)
         results.shakemaps["stddevs"][imt] = shakemap.to_geotiff(
@@ -151,6 +155,13 @@ def shakemaps_from_quakeml(
         results.contours["stddevs"][imt] = shakemap.get_contours(
             imt, stddev_shakemaps, contour_levels_stddev, is_stddev=True
         )
+        # Retrieve summary statistics
+        results.statistics["mean"][imt] = {
+            "maximum": np.exp(np.max(mean_shakemaps[imt])),
+            "median": np.exp(np.mean(mean_shakemaps[imt])),
+            "minimum": np.exp(np.min(mean_shakemaps[imt])),
+            "IQR": stats.iqr(mean_shakemaps[imt]),
+        }
         if export_folder:
             filestem = "{:s}_{:s}".format(results.earthquake.id, imt)
             # Export the bytes to raster files

tests/workflows_test.py

@@ -3,6 +3,7 @@ Test cases for end-to-end workflows
 """
 import os
 import unittest
+import numpy as np
 from geopandas import GeoDataFrame
 from shakyground2.workflows import shakemaps_from_quakeml, ShakemapWorkflowResult
 
@@ -41,6 +42,9 @@ class GeofonWorkflowTestCase(unittest.TestCase):
             for maptype in results.shakemaps:
                 self.assertIsInstance(results.shakemaps[maptype][imt], bytes)
             self.assertIsInstance(results.contours["mean"][imt], GeoDataFrame)
+            for stat_type in results.statistics["mean"][imt]:
+                # Verify that none of the statistics are nans
+                self.assertFalse(np.isnan(results.statistics["mean"][imt][stat_type]))
 
     def test_workflow_with_exports(self):
         test_folder = os.path.join(DATA_PATH, "tmp_shakemaps_from_quakeml")