Merge pull request #67 from g-weatherill/trellis_test

Adds unit tests for trellis plots

Author: Graeme Weatherill

smtk/trellis/configure.py

@@ -558,13 +558,11 @@ def get_gsim_contexts(self):
         setattr(dctx, "rvolc", np.zeros_like(self.target_sites.mesh.lons))
         # Sites
         sctx = SitesContext()
-        key_list = ['_vs30', '_vs30measured', '_z1pt0', '_z2pt5', '_backarc']
+        key_list = ["lons", "lats", "vs30", "vs30measured", "z1pt0", "z2pt5",
+                    "backarc"]
         for key in key_list:
-            setattr(sctx, key[1:], getattr(self.target_sites, key))
-        for key in ['lons', 'lats']:
-            setattr(sctx, key, getattr(self.target_sites, key))
+            setattr(sctx, key, self.target_sites.array[key])
 
-        
         # Rupture
         rctx = RuptureContext()
         setattr(rctx, 'mag', self.magnitude)
@@ -727,7 +725,7 @@ def get_target_sites_point(self, distance, distance_type, vs30,
         :param floar z2pt5:
             Depth to 2.5 km/s interface
         """
-        if not distance_type in POINT_AT_MAPPING.keys():
+        if not distance_type in list(POINT_AT_MAPPING.keys()):
             raise ValueError("Distance type must be one of: Rupture ('rrup'), "
                              "Joyner-Boore ('rjb'), Epicentral ('repi') or "
                              "Hypocentral ('rhyp')")
@@ -799,7 +797,7 @@ def filter_hanging_wall(self, filter_type=None):
         if not filter_type:
             # Considers both footwall and hanging wall 
             return self.target_sites
-        elif not filter_type in ['HW', 'FW']:
+        elif not filter_type in ('HW', 'FW'):
             raise ValueError('Hanging wall filter must be either "HW" or "FW"')
         else:
             pass

smtk/trellis/trellis_plots.py

@@ -361,6 +361,7 @@ def from_rupture_model(cls, rupture, gsims, imts, stddevs='Total',
         kwargs.setdefault('distance_type', "rjb")
         kwargs.setdefault('xlim', None)
         kwargs.setdefault('ylim', None)
+        print(rupture.__class__)
         assert isinstance(rupture, GSIMRupture)
         magnitudes = [rupture.magnitude]
         sctx, rctx, dctx = rupture.get_gsim_contexts()
@@ -548,9 +549,9 @@ def _get_ylabel(self, i_m):
             units = PLOT_UNITS[i_m]
         return "Mean {:s} ({:s})".format(i_m, units)
 
-    def to_json(self):
+    def to_dict(self):
         """
-        Serializes the ground motion values to json
+        Parse the ground motion values to a dictionary
         """
         gmvs = self.get_ground_motion_values()
         gmv_dict = OrderedDict([
@@ -572,7 +573,13 @@ def to_json(self):
                         iml_to_list.append(val)
                 gmv_dict["figures"][imt].append((gsim, iml_to_list))
             gmv_dict["figures"][imt] = OrderedDict(gmv_dict["figures"][imt])
-        return json.dumps(gmv_dict)
+        return gmv_dict
+
+    def to_json(self):
+        """
+        Serializes the ground motion values to json
+        """
+        return json.dumps(self.to_dict())
 
     def get_ground_motion_values(self):
         """
@@ -773,7 +780,7 @@ class DistanceIMTTrellis(MagnitudeIMTTrellis):
     def __init__(self, magnitudes, distances, gsims, imts, params, 
             stddevs="Total", **kwargs):
         """
-        Instandi 
+        Instantiation 
         """
         if isinstance(magnitudes, float):
             magnitudes = [magnitudes]
@@ -881,9 +888,9 @@ def _get_ylabel(self, i_m):
             units = PLOT_UNITS[i_m]
         return "Mean {:s} ({:s})".format(i_m, units)
 
-    def to_json(self):
+    def to_dict(self):
         """
-        Exports ground motion values to json
+        Parses the ground motion values to a dictionary
         """
         gmvs = self.get_ground_motion_values()
         dist_label = "{:s} (km)".format(DISTANCE_LABEL_MAP[self.distance_type])
@@ -900,7 +907,13 @@ def to_json(self):
                 imt_dict["yvalues"].append((gsim, data))
             imt_dict["yvalues"] = OrderedDict(imt_dict["yvalues"])
             gmv_dict["figures"][imt] = imt_dict
-        return json.dumps(gmv_dict)
+        return gmv_dict
+
+    def to_json(self):
+        """
+        Exports ground motion values to json
+        """
+        return json.dumps(self.to_dict())
 
     def pretty_print(self, filename=None, sep=","):
         """
@@ -1370,9 +1383,9 @@ def _set_labels(self, i_m, ax):
         ax.set_xlabel("Period (s)", fontsize=14)
         ax.set_ylabel("Sa (g)", fontsize=14)
 
-    def to_json(self):
+    def to_dict(self):
         """
-        Export ground motion values dictionary to json
+        Export ground motion values to a dictionary
         """
         gmvs = self.get_ground_motion_values()
         periods = [float(val.split("SA(")[1].rstrip(")"))
@@ -1405,7 +1418,13 @@ def to_json(self):
                         else:
                             gmv_dict["figures"][pos_name]["yvalues"][gsim].\
                                 append(None)
-        return json.dumps(gmv_dict)
+        return gmv_dict
+
+    def to_json(self):
+        """
+        Exports the ground motion values to json
+        """
+        return json.dumps(self.to_dict())
 
     def _get_ylabel(self, i_m):
         """

tests/trellis/data/test_distance_imt_trellis.json

@@ -0,0 +1,244 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2014-2017 GEM Foundation and G. Weatherill
+#
+# 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/>.
+"""
+Tests for generation of data for trellis plots
+"""
+import unittest
+import os
+import json
+import numpy as np
+import smtk.trellis.trellis_plots as trpl
+import smtk.trellis.configure as rcfg
+
+
+BASE_DATA_PATH = os.path.join(os.path.dirname(__file__), "data")
+
+
+class BaseTrellisTest(unittest.TestCase):
+    """
+    This core test is designed to run a series of trellis plot calculations
+    and ensure compatibility with previously generated results
+    """
+
+    TEST_FILE = None
+
+    def setUp(self):
+        self.imts = ["PGA", "SA(0.2)", "SA(2.0)", "SA(3.0)"]
+        self.periods = [0.05, 0.075, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16,
+                        0.17, 0.18, 0.19, 0.20, 0.22, 0.24, 0.26, 0.28, 0.30,
+                        0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, 0.48,
+                        0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
+                        1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
+                        3.0, 4.001, 5.0, 7.5, 10.0]
+
+        self.gsims = ["AkkarBommer2010", "CauzziFaccioli2008",
+                      "ChiouYoungs2008", "ZhaoEtAl2006Asc", "AkkarEtAlRjb2014",
+                      "BindiEtAl2014Rjb", "CauzziEtAl2014", "DerrasEtAl2014",
+                      "AbrahamsonEtAl2014", "BooreEtAl2014", "ChiouYoungs2014",
+                      "CampbellBozorgnia2014", "KothaEtAl2016Italy",
+                      "KothaEtAl2016Other", "KothaEtAl2016Turkey",
+                      "ZhaoEtAl2016Asc", "BindiEtAl2017Rjb"]
+
+
+class DistanceTrellisTest(BaseTrellisTest):
+    TEST_FILE = "test_distance_imt_trellis.json"
+
+    def compare_jsons(self, old, new):
+        """
+        Compares the json data from file with the new data from trellis plot
+        """
+        # Check x-labels are the same
+        self.assertEqual(old["xlabel"], new["xlabel"])
+        # Check x-values are the same
+        np.testing.assert_array_almost_equal(old["xvalues"], new["xvalues"], 7)
+        for imt in old["figures"]:
+            self.assertEqual(old["figures"][imt]["ylabel"],
+                             new["figures"][imt]["ylabel"])
+            for gsim in old["figures"][imt]["yvalues"]:
+                np.testing.assert_array_almost_equal(
+                    old["figures"][imt]["yvalues"][gsim],
+                    new["figures"][imt]["yvalues"][gsim], 7)
+
+    def _run_trellis(self, rupture):
+        """
+        Executes the trellis plotting - for mean
+        """
+        return trpl.DistanceIMTTrellis.from_rupture_model(rupture,
+                                                          self.gsims,
+                                                          self.imts,
+                                                          distance_type="rrup")
+
+    def test_distance_imt_trellis(self):
+        """
+        Tests the DistanceIMT trellis data generation
+        """
+        reference = json.load(open(
+            os.path.join(BASE_DATA_PATH, self.TEST_FILE), "r"))
+        # Setup rupture
+        rupture = rcfg.GSIMRupture(6.5, 60., 1.5,
+                                   hypocentre_location=(0.5, 0.5))
+        _ = rupture.get_target_sites_line(250.0, 1.0, 800.0)
+        # Get trellis calculations
+        trl = self._run_trellis(rupture)
+        # Parse the json formatted string to a dictionary string
+        results = json.loads(trl.to_json())
+        # Compare the two dictionaries
+        self.compare_jsons(reference, results)
+
+
+class DistanceSigmaTrellisTest(DistanceTrellisTest):
+    TEST_FILE = "test_distance_sigma_imt_trellis.json"
+
+    def _run_trellis(self, rupture):
+        """
+        Executes the trellis plotting - for standard deviation
+        """
+        return trpl.DistanceSigmaIMTTrellis.from_rupture_model(
+            rupture,
+            self.gsims,
+            self.imts,
+            distance_type="rrup")
+
+
+class MagnitudeTrellisTest(BaseTrellisTest):
+    TEST_FILE = "test_magnitude_imt_trellis.json"
+
+    def _run_trellis(self, magnitudes, distance, properties):
+        """
+        Executes the trellis plotting - for mean
+        """
+        return trpl.MagnitudeIMTTrellis.from_rupture_model(properties,
+                                                           magnitudes,
+                                                           distance,
+                                                           self.gsims,
+                                                           self.imts)
+
+    def compare_jsons(self, old, new):
+        """
+        Compares the json data from file with the new data from trellis plot
+        """
+        self.assertEqual(old["xlabel"], new["xlabel"])
+        np.testing.assert_array_almost_equal(old["xvalues"], new["xvalues"], 7)
+        for imt in old["figures"]:
+            for gsim in old["figures"][imt]:
+                if gsim == u"ylabel":
+                    self.assertEqual(old["figures"][imt][gsim],
+                                     new["figures"][imt][gsim])
+                else:
+                    np.testing.assert_array_almost_equal(
+                        old["figures"][imt][gsim],
+                        new["figures"][imt][gsim], 7)
+
+    def test_magnitude_imt_trellis(self):
+        """
+        Tests the MagnitudeIMT trellis data generation
+        """
+        reference = json.load(open(
+            os.path.join(BASE_DATA_PATH, self.TEST_FILE), "r"))
+        magnitudes = np.arange(4., 8.1, 0.1)
+        distance = 20.
+        properties = {"dip": 60.0, "rake": -90.0, "aspect": 1.5, "ztor": 0.0,
+                      "vs30": 800.0, "backarc": False, "z1pt0": 50.0,
+                      "z2pt5": 1.0, "line_azimuth": 90.0}
+        trl = self._run_trellis(magnitudes, distance, properties)
+        results = json.loads(trl.to_json())
+        self.compare_jsons(reference, results)
+
+
+class MagnitudeSigmaTrellisTest(MagnitudeTrellisTest):
+    TEST_FILE = "test_magnitude_sigma_imt_trellis.json"
+
+    def _run_trellis(self, magnitudes, distance, properties):
+        """
+        Executes the trellis plotting - for standard deviation
+        """
+        return trpl.MagnitudeSigmaIMTTrellis.from_rupture_model(properties,
+                                                                magnitudes,
+                                                                distance,
+                                                                self.gsims,
+                                                                self.imts)
+
+
+class MagnitudeDistanceSpectraTrellisTest(BaseTrellisTest):
+    TEST_FILE = "test_magnitude_distance_spectra_trellis.json"
+
+    def compare_jsons(self, old, new):
+        """
+        Compares the MagnitudeDistanceSpectra jsons
+        """
+        self.assertEqual(old["xlabel"], new["xlabel"])
+        self.assertEqual(old["ylabel"], new["ylabel"])
+        np.testing.assert_array_almost_equal(old["xvalues"], new["xvalues"], 7)
+        for key in old["figures"]:
+            self.assertAlmostEqual(old["figures"][key]["magnitude"],
+                                   new["figures"][key]["magnitude"], 7)
+            self.assertAlmostEqual(old["figures"][key]["distance"],
+                                   new["figures"][key]["distance"], 7)
+            self.assertEqual(old["figures"][key]["row"],
+                             new["figures"][key]["row"])
+            self.assertEqual(old["figures"][key]["column"],
+                             new["figures"][key]["column"])
+            for gsim in old["figures"][key]["yvalues"]:
+                old_vals = np.array(old["figures"][key]["yvalues"][gsim])
+                new_vals = np.array(new["figures"][key]["yvalues"][gsim])
+                if old_vals.dtype == "O":
+                    # Has None Values - compare element by element
+                    for old_val, new_val in zip(old_vals, new_vals):
+                        if old_val and new_val:
+                            self.assertAlmostEqual(old_val, new_val, 7)
+                        else:
+                            self.assertEqual(old_val, new_val)
+                else:
+                    np.testing.assert_array_almost_equal(old_vals, new_vals, 7)
+
+    def _run_trellis(self, magnitudes, distances, properties):
+        """
+        Executes the trellis plotting - for mean
+        """
+        return trpl.MagnitudeDistanceSpectraTrellis.from_rupture_model(
+            properties, magnitudes, distances, self.gsims, self.periods,
+            distance_type="rrup")
+
+    def test_magnitude_distance_spectra_trellis(self):
+        """
+        Tests the MagnitudeDistanceSpectra Trellis data generation
+        """
+        reference = json.load(open(
+            os.path.join(BASE_DATA_PATH, self.TEST_FILE), "r"))
+        properties = {"dip": 60.0, "rake": -90.0, "aspect": 1.5, "ztor": 0.0,
+                      "vs30": 800.0, "backarc": False, "z1pt0": 50.0,
+                      "z2pt5": 1.0}
+        magnitudes = [4.0, 5.0, 6.0, 7.0]
+        distances = [5., 20., 50., 150.0]
+        trl = self._run_trellis(magnitudes, distances, properties)
+        results = json.loads(trl.to_json())
+        self.compare_jsons(reference, results)
+
+
+class MagnitudeDistanceSpectraSigmaTrellisTest(
+        MagnitudeDistanceSpectraTrellisTest):
+    TEST_FILE = "test_magnitude_distance_spectra_sigma_trellis.json"
+
+    def _run_trellis(self, magnitudes, distances, properties):
+        """
+        Executes the trellis plotting - for standard deviation
+        """
+        return trpl.MagnitudeDistanceSpectraSigmaTrellis.from_rupture_model(
+            properties, magnitudes, distances, self.gsims, self.periods,
+            distance_type="rrup")