Merge pull request #6660 from swiss-seismological-service/feature/mmi…

…-shakemaps

Feature/mmi shakemaps

debian/changelog

@@ -1,6 +1,7 @@
   [Nicolas Schmid]
   * Improve performance for ShakeMap calculations when spatialcorr and crosscorr
     are both set to 'no'
+  * Add feature to do ShakeMap calculations for vulnerability models using MMI.
 
   [Michele Simionato]
   * Supported exposures with generic CSV fields thanks to the `exposureFields`

openquake/calculators/base.py

@@ -1183,15 +1183,33 @@ def read_shakemap(calc, haz_sitecol, assetcol):
         logging.info('Extracted %d assets', len(assetcol))
 
     # assemble dictionary to decide on the calculation method for the gmfs
-    if oq.spatial_correlation != 'no' or oq.cross_correlation != 'no':
-        # cross correlation and/or spatial correlation after S&H
-        gmf_dict = {'kind': 'Silva&Horspool',
-                    'spatialcorr': oq.spatial_correlation,
-                    'crosscorr': oq.cross_correlation,
-                    'cholesky_limit': oq.cholesky_limit}
+    if 'MMI' in oq.imtls:
+        # calculations with MMI should be executed
+        if len(oq.imtls) == 1:
+            # only MMI intensities
+            if oq.spatial_correlation != 'no' or oq.cross_correlation != 'no':
+                logging.warning('Calculations with MMI intensities do not '
+                                'support correlation. No correlations '
+                                'are applied.')
+
+            gmf_dict = {'kind': 'mmi'}
+        else:
+            # there are also other intensities than MMI
+            raise RuntimeError(
+                'There are the following intensities in your model: %s '
+                'Models mixing MMI and other intensities are not supported. '
+                % ', '.join(oq.imtls.keys()))
     else:
-        # no correlation required, basic calculation is faster
-        gmf_dict = {'kind': 'basic'}
+        # no MMI intensities, calculation with or without correlation
+        if oq.spatial_correlation != 'no' or oq.cross_correlation != 'no':
+            # cross correlation and/or spatial correlation after S&H
+            gmf_dict = {'kind': 'Silva&Horspool',
+                        'spatialcorr': oq.spatial_correlation,
+                        'crosscorr': oq.cross_correlation,
+                        'cholesky_limit': oq.cholesky_limit}
+        else:
+            # no correlation required, basic calculation is faster
+            gmf_dict = {'kind': 'basic'}
 
     logging.info('Building GMFs')
     with calc.monitor('building/saving GMFs'):

openquake/hazardlib/shakemap/gmfs.py

@@ -173,6 +173,8 @@ def calculate_gmfs_sh(kind, shakemap, imts, Z, mu, spatialcorr,
     :param crosscorr: 'no', 'yes' or 'full'
     :returns: F(Z, mu) to calculate gmfs
     """
+    # make sure all imts used have a period, needed for correlation
+    imts = [im for im in imts if hasattr(im, 'period')]
     # checks
     N = len(shakemap)
     M = len(imts)
@@ -217,6 +219,26 @@ def calculate_gmfs_basic(kind, shakemap, imts, Z, mu):
     return numpy.exp(sig @ Z + numpy.log(mu)) / PCTG
 
 
+@ calculate_gmfs.add('mmi')
+def calculate_gmfs_mmi(kind, shakemap, imts, Z, mu):
+    """
+    Basic calculation method to sample data from shakemap values
+    given mmi intensities.
+
+    :param shakemap: site coordinates with shakemap values
+    :param imts: list of required imts
+    :returns: F(Z, mu) to calculate gmfs
+    """
+    try:
+        # create diag matrix with std values
+        std = numpy.array(shakemap['std']['MMI'])
+        sig = numpy.diag(std.flatten())  # shape (M*N, M*N)
+    except ValueError as e:
+        raise ValueError('No stds for MMI intensities supplied, only %s' %
+                         ', '.join(shakemap['std'].dtype.names)) from e
+    return sig @ Z + mu
+
+
 def to_gmfs(shakemap, gmf_dict, site_effects, trunclevel,
             num_gmfs, seed, imts=None):
     """

openquake/hazardlib/shakemap/parsers.py

@@ -36,16 +36,19 @@
 SHAKEMAP_URL = US_GOV + '/fdsnws/event/1/query?eventid={}&format=geojson'
 F32 = numpy.float32
 SHAKEMAP_FIELDS = set(
-    'LON LAT SVEL PGA PSA03 PSA10 PSA30 STDPGA STDPSA03 STDPSA10 STDPSA30'
+    'LON LAT SVEL MMI PGA PSA03 PSA10 PSA30 '
+    'STDMMI STDPGA STDPSA03 STDPSA10 STDPSA30'
     .split())
 FIELDMAP = {
     'LON': 'lon',
     'LAT': 'lat',
     'SVEL': 'vs30',
+    'MMI': ('val', 'MMI'),
     'PGA': ('val', 'PGA'),
     'PSA03': ('val', 'SA(0.3)'),
     'PSA10': ('val', 'SA(1.0)'),
     'PSA30': ('val', 'SA(3.0)'),
+    'STDMMI': ('std', 'MMI'),
     'STDPGA': ('std', 'PGA'),
     'STDPSA03': ('std', 'SA(0.3)'),
     'STDPSA10': ('std', 'SA(1.0)'),

openquake/hazardlib/tests/shakemap/shakemap_test.py

@@ -28,7 +28,7 @@ def mean_std(shakemap, site_effects):
                      'spatialcorr': 'yes', 'crosscorr': 'yes'})
     _, gmfs = to_gmfs(
         shakemap, gmf_dict, site_effects, trunclevel=3,
-        num_gmfs=1000, seed=42)
+        num_gmfs=1000, seed=42, imts=['PGA', 'SA(0.3)', 'SA(1.0)', 'SA(3.0)'])
     return gmfs.mean(axis=1), numpy.log(gmfs).std(axis=1)
 
 
@@ -159,3 +159,26 @@ def test_from_files(self):
                          [0.6077153, 0.6661571, 0.6296381, 0.668559],
                          [0.6146356, 0.6748830, 0.6714424, 0.6613612],
                          [0.5815353, 0.6460007, 0.6491335, 0.6603457]])
+
+    def test_for_mmi(self):
+        f1 = 'file://' + os.path.join(CDIR, 'ghorka_grid.xml')
+        f2 = 'file://' + os.path.join(CDIR, 'ghorka_uncertainty.xml')
+        uridict = dict(kind='usgs_xml', grid_url=f1, uncertainty_url=f2)
+        sitecol, shakemap, *_ = get_sitecol_shakemap(uridict.pop('kind'),
+                                                     uridict, imt_dt.names)
+        n = 4  # number of sites
+        self.assertEqual(len(sitecol), n)
+
+        _, gmfs = to_gmfs(shakemap, {'kind': 'mmi'}, False,
+                          trunclevel=3, num_gmfs=1000, seed=42, imts=['MMI'])
+
+        gmf_by_imt = gmfs.mean(axis=1)
+        std_by_imt = gmfs.std(axis=1)
+        aae(gmf_by_imt, [[3.80704848],
+                         [3.89791949],
+                         [3.88040454],
+                         [3.93584243]])
+        aae(std_by_imt, [[0.71068558],
+                         [0.72233552],
+                         [0.72033749],
+                         [0.69906908]])