Merge remote-tracking branch 'upstream/master'

shakyground2/site_model.py

+"""
+Core class to define a set of target sites for the shakemap calculation
+"""
+from __future__ import annotations
+import warnings
+import numpy as np
+import pandas as pd
+from typing import Tuple, Optional, Dict
+from openquake.hazardlib.site import SiteCollection
+
+
+# Note that OpenQuake unit conventions are being adopted here for consistency
+# with the OpenQuake software itself
+SITE_PROPERTIES = {
+    "sids": np.int64,  # Unique site identifiers
+    "lon": np.float64,  # Longitudes (in decimal degrees)
+    "lat": np.float64,  # Latitudes (in decimal degrees)
+    "depth": np.float64,  # Site depth (in km, negative indicates above ground)
+    "vs30": np.float64,  # Vs30 of the sites (in m/s)
+    "vs30measured": np.bool,  # Vs30 is measured (True) or inferred (False)
+    "z1pt0": np.float64,  # Depth (m) to the Vs 1.0 km/s layer
+    "z2pt5": np.float64,  # Depth (km) to the Vs 2.5 km/s layer
+    "xvf": np.float64,  # Distance (km) to the volcanic front
+    "backarc": np.bool,  # Site is in the subduction backarc (True) or else
+    "region": np.int32,  # Region to which the site belongs
+    "geology": (np.string_, 20)  # Geological classification for the site
+}  
+
+
+# In some cases reasonable default values can be used for relevant ground motion models
+# unless specifically defined by a user
+SITE_DEFAULTS = {
+    "vs30measured": True,
+    "xvf": 150.0,
+    "region": 0,
+    "geology": "UNCLASSIFIED",
+    "backarc": False,
+}
+
+
+class SiteModel(object):
+    """
+    Class to manage the definition of the site properties for application in a Shakemap
+
+    Attributes:
+        site_array: Array of sites and their corresponding properties as numpy.ndarray
+        bbox_properties: Properties of the bounding box used to create the site array
+                         (where relevant)
+    """
+
+    def __init__(self, site_array: np.ndarray, bbox_properties: Dict = {}):
+        self.site_array = site_array
+        self.bbox_properties = bbox_properties
+
+    def __len__(self):
+        # Returns the number of sites in the site array
+        return self.site_array.shape[0]
+
+    def __repr__(self):
+        # Returns an informative string regarding the composition of the sites,
+        # including the bounding box information when relevant
+        info_string = "SiteModel(%g sites)" % len(self)
+        if self.bbox_properties:
+            bbox_string = " BBOX: %.3f/%.3f/%.3f/%.3f DX: %.3e DY: %.3e" % (
+                self.bbox_properties["bbox"][0],
+                self.bbox_properties["bbox"][1],
+                self.bbox_properties["bbox"][2],
+                self.bbox_properties["bbox"][3],
+                self.bbox_properties["spcx"],
+                self.bbox_properties["spcy"],
+            )
+            info_string += bbox_string
+        return info_string
+
+    @property
+    def dataframe(self):
+        # Return the site array as a dataframe
+        return pd.DataFrame(self.site_array)
+
+    def get_site_collection(self, idx: Optional[np.ndarray] = None):
+        """
+        Convert the site model to an OpenQuake SiteCollection, filtering the site
+        model where defined
+        Args:
+            idx: Index of rows in the site model to be returned in the SiteCollection, input
+                 as either a logical vector (numpy.ndarray) or a vector of indices
+                 (numpy.ndarray)
+        Returns:
+            Site collection as instance of :class:`openquake.hazardlib.site.SiteCollection`
+        """
+        # The site collection is usually instantiated with a list of single Site objects
+        # and then passed to the site array. This step is bypassed and the SiteCollection's
+        # array is build directly from the site arrays
+        site_col = object.__new__(SiteCollection)
+        site_col.complete = site_col
+        if idx is not None:
+            # Filter the site array
+            site_col.array = self.site_array[idx]
+        else:
+            site_col.array = self.site_array.copy()
+        site_col.array.flags.writeable = False
+        return site_col
+
+    @classmethod
+    def from_bbox(
+        cls,
+        bbox: Tuple[float, float, float, float],
+        spcx: float,
+        spcy: float,
+        vs30: float,
+        vs30measured: bool = True,
+        z1pt0: Optional[float] = None,
+        z2pt5: Optional[float] = None,
+        xvf: float = 150.0,
+        backarc: bool = False,
+        region: int = 0,
+        geology: np.string_ = "UNCLASSIFIED",
+    ) -> SiteModel:
+        """
+        Builds a site model from a bounding box and scalar site properties to be
+        fixed for all sites
+
+        Args:
+            bbox:
+                Bounding box of a region as [llon, llat, ulon, ulat]
+            spcx:
+                Fixed longitude spacing of successive points in decimal degrees
+            spcy:
+                Fixed latitude spacing of successive points in decimal degrees
+            vs30:
+                Fixed Vs30 assigned to all points
+            vs30measured:
+                Indicates whether the Vs30 refers to a "measured" site condition
+                (True) or "inferred" (False)
+            z1pt0:
+                Depth to Vs 1.0 km layer (in m)
+            z2pt5:
+                Depth to Vs 2.5 km layer (in km)
+
+        Returns:
+            Instance of :class`shakygroundv2.site_model.SiteModel`
+        """
+
+        # Generate the mesh of points
+        llon, llat, ulon, ulat = bbox
+        lons, lats = np.meshgrid(
+            np.arange(llon, ulon + spcx, spcx), np.arange(llat, ulat + spcy, spcy)
+        )
+        nrow, ncol = lons.shape
+        nsites = nrow * ncol
+        self = object.__new__(cls)
+        # To re-shape or export of the shakemap back into a regular grid you need the original
+        # configuration properties used to build the site model
+        bbox_properties = {
+            "bbox": bbox,
+            "spcx": spcx,
+            "spcy": spcy,
+            "ncol": ncol,
+            "nrow": nrow,
+        }
+        # Site array has a set of defined datatypes
+        site_dtype = np.dtype(list(SITE_PROPERTIES.items()))
+        site_array = np.zeros(nsites, site_dtype)
+        site_array["sids"] = np.arange(nsites, dtype=np.int64)
+        site_array["lon"] = lons.flatten()
+        site_array["lat"] = lats.flatten()
+        site_array["vs30"] = vs30 * np.ones(nsites)
+        site_array["vs30measured"] = np.ones(nsites, dtype=bool)
+        # In the cases of Z1.0 and Z2.5 these can be calculated from
+        # Vs30 if not defined explicitly
+        if z1pt0 is None:
+            site_array["z1pt0"] = cls.vs30_to_z1pt0_cy14(site_array["vs30"])
+        else:
+            site_array["z1pt0"] = z1pt0 * np.ones(nsites)
+
+        if z2pt5 is None:
+            site_array["z2pt5"] = cls.vs30_to_z2pt5_cb14(site_array["vs30"])
+        else:
+            site_array["z2pt5"] = z2pt5 * np.ones(nsites)
+
+        # Get other default options from kwargs if specified
+        site_array["xvf"] = xvf * np.ones(nsites, dtype=SITE_PROPERTIES["xvf"])
+        site_array["backarc"] = self._get_boolean_array(nsites, backarc)
+        site_array["region"] = region * np.ones(nsites, dtype=SITE_PROPERTIES["region"])
+        site_array["geology"] = self._get_string_array(
+            nsites, geology, SITE_PROPERTIES["geology"]
+        )
+        return cls(site_array, bbox_properties=bbox_properties)
+
+    @classmethod
+    def from_dataframe(cls, dframe: pd.DataFrame) -> SiteModel:
+        """
+        Builds the site collection directly from a pandas dataframe
+
+        Args:
+            dframe: Input site model as a pandas.Dataframe
+
+        Returns:
+            Instance of :class`shakygroundv2.site_model.SiteModel`
+        """
+        nsites = dframe.shape[0]
+        site_dtype = np.dtype(list(SITE_PROPERTIES.items()))
+        site_array = np.zeros(nsites, site_dtype)
+        site_array["sids"] = np.arange(nsites)
+        for col in dframe.columns:
+            if col not in SITE_PROPERTIES:
+                warnings.warn(
+                    "Input site property %s not recognised" " for use - skipping" % col
+                )
+                continue
+            site_array[col] = dframe[col].to_numpy(dtype=SITE_PROPERTIES[col])
+        # Fill in any missing information from the defaults
+        if "z1pt0" not in dframe.columns:
+            site_array["z1pt0"] = cls.vs30_to_z1pt0_cy14(site_array["vs30"])
+        if "z2pt5" not in dframe.columns:
+            site_array["z2pt5"] = cls.vs30_to_z2pt5_cb14(site_array["vs30"])
+        for key, value in SITE_DEFAULTS.items():
+            if key not in dframe.columns:
+                # Use the default
+                if SITE_PROPERTIES[key] in ((np.bytes_, 20),):
+                    site_array[key] = cls._get_string_array(
+                        nsites, value, SITE_PROPERTIES[key]
+                    )
+                elif SITE_PROPERTIES[key] == np.bool:
+                    site_array[key] = cls._get_boolean_array(nsites, value)
+                else:
+                    site_array[key] = value * np.ones(
+                        nsites, dtype=SITE_PROPERTIES[key]
+                    )
+        return cls(site_array)
+
+    @staticmethod
+    def _get_string_array(num_elem: int, key: str, dtype) -> np.ndarray:
+        """
+        Returns an array of constant string values
+
+        Args:
+            num_elem: Number of elements in the array
+            key: String to fill the array
+            dtype: Numpy datatype for the string
+
+        Returns:
+            Numpy string vector
+        """
+        arr = np.empty(num_elem, dtype=dtype)
+        arr[:] = key
+        return arr
+
+    @staticmethod
+    def _get_boolean_array(num_elem: int, value: bool) -> np.ndarray:
+        """
+        Returns an array of boolean values
+
+        Args:
+            num_elem: Number of elements in the array
+            value: Indicates if the array should be True of False
+        """
+        arr = np.empty(num_elem, dtype=bool)
+        arr[:] = value
+        return arr
+
+    @staticmethod
+    def vs30_to_z1pt0_cy14(vs30: np.ndarray, japan: bool = False) -> np.ndarray:
+        """
+        Returns the estimate depth to the 1.0 km/s velocity layer based on Vs30
+        from Chiou & Youngs (2014) California model or Japan model
+
+        Args:
+            vs30:
+                Input Vs30 values in m/s
+            japan:
+                If true then this returns the Z1.0 values according to the Japan model,
+                otherwise the California model
+        Returns:
+                Z1.0 values in m
+        """
+        if japan:
+            c1 = 412.0 ** 2.0
+            c2 = 1360.0 ** 2.0
+            return np.exp(
+                (-5.23 / 2.0) * np.log((np.power(vs30, 2.0) + c1) / (c2 + c1))
+            )
+        else:
+            c1 = 571.0 ** 4.0
+            c2 = 1360.0 ** 4.0
+            return np.exp((-7.15 / 4.0) * np.log((vs30 ** 4.0 + c1) / (c2 + c1)))
+
+    @staticmethod
+    def vs30_to_z2pt5_cb14(vs30: np.ndarray, japan: bool = False) -> np.ndarray:
+        """
+        Converts vs30 to depth to 2.5 km/s interface using model proposed by
+        Campbell & Bozorgnia (2014)
+
+        Args:
+            vs30:
+                Input Vs30 values in m/s
+            japan:
+                If true then this returns the Z1.0 values according to the Japan model,
+                otherwise the California model
+        Returns:
+                Z2.5 values in km
+        """
+        if japan:
+            return np.exp(5.359 - 1.102 * np.log(vs30))
+        else:
+            return np.exp(7.089 - 1.144 * np.log(vs30))

tests/test_site_model.py

+"""
+Test suite for the site model class
+"""
+import unittest
+import warnings
+import numpy as np
+from openquake.hazardlib.geo import Point
+from openquake.hazardlib.site import Site, SiteCollection
+from shakyground2.site_model import SiteModel, SITE_PROPERTIES
+
+
+class SiteModelTestCase(unittest.TestCase):
+    def setUp(self):
+        self.bbox = [0.0, 0.0, 2.0, 2.0]
+        self.spcx = 1.0
+        self.spcy = 1.0
+        self.vs30 = 1000.0
+        self.z1pt0 = 0.0
+        self.z2pt5 = 1.0
+
+    def test_length(self):
+        dummy_array = np.zeros([20, 10])
+        sc1 = SiteModel(dummy_array)
+        self.assertEqual(len(sc1), 20)
+
+    def test_string(self):
+        # No bounding box case
+        dummy_array = np.zeros([20, 10])
+        sc1 = SiteModel(dummy_array)
+        self.assertEqual(str(sc1), "SiteModel(20 sites)")
+        # Bounding box case
+        bbox_properties = {"bbox": [1.0, 2.0, 3.0, 4.0], "spcx": 0.5, "spcy": 0.5}
+        sc1 = SiteModel(dummy_array, bbox_properties)
+        self.assertEqual(
+            str(sc1),
+            "SiteModel(20 sites) BBOX: 1.000/2.000/3.000/4.000 DX: 5.000e-01 DY: 5.000e-01",
+        )
+
+    def test_vs30_to_z1pt0_cy14(self):
+        # California case
+        vs30 = np.array([300.0, 500.0, 800.0, 1200.0])
+        z1pt0 = SiteModel.vs30_to_z1pt0_cy14(vs30)
+        np.testing.assert_array_almost_equal(
+            z1pt0, np.array([458.77871089, 228.88509539, 31.07014865, 2.36375239])
+        )
+        # Japan case
+        z1pt0 = SiteModel.vs30_to_z1pt0_cy14(vs30, japan=True)
+        np.testing.assert_array_almost_equal(
+            z1pt0, np.array([213.3479186, 60.81637035, 10.90926822, 1.80907214])
+        )
+
+    def test_vs30_to_z2pt5_cb14(self):
+        # California case
+        vs30 = np.array([300.0, 500.0, 800.0, 1200.0])
+        z2pt5 = SiteModel.vs30_to_z2pt5_cb14(vs30)
+        np.testing.assert_array_almost_equal(
+            z2pt5, np.array([1.75746574, 0.97969727, 0.57224056, 0.35985723])
+        )
+        # Japan case
+        z2pt5 = SiteModel.vs30_to_z2pt5_cb14(vs30, japan=True)
+        np.testing.assert_array_almost_equal(
+            z2pt5, np.array([0.39591003, 0.22548579, 0.13433184, 0.08592636])
+        )
+
+    def test_get_string_array(self):
+        string_array = SiteModel._get_string_array(3, "XYZ", (np.string_, 3))
+        np.testing.assert_array_equal(
+            string_array, np.array([b"XYZ", b"XYZ", b"XYZ"], dtype=(np.string_, 3))
+        )
+
+    def test_get_boolean_array(self):
+        # True case
+        true_array = SiteModel._get_boolean_array(3, True)
+        np.testing.assert_array_almost_equal(true_array, np.ones(3, dtype=bool))
+        # False case
+        true_array = SiteModel._get_boolean_array(3, False)
+        np.testing.assert_array_almost_equal(true_array, np.zeros(3, dtype=bool))
+
+    def _compare_constant_vectors(self, test_vector, value, nelem, dtype):
+        # For vectors with constant values then test the vector against the
+        # expected target vector with a specific value, size and datatype
+        target_vector = np.array([value for i in range(nelem)], dtype)
+        if dtype is np.float64:
+            np.testing.assert_array_almost_equal(test_vector, target_vector)
+        else:
+            np.testing.assert_array_equal(test_vector, target_vector)
+
+    def test_build_from_bbox_comprehensive(self):
+        # Build the site model with a comprehensive set of inputs
+
+        sm1 = SiteModel.from_bbox(
+            self.bbox,
+            self.spcx,
+            self.spcy,
+            self.vs30,
+            z1pt0=self.z1pt0,
+            z2pt5=self.z2pt5,
+            xvf=100.0,
+            backarc=False,
+            region=1,
+            geology="ABC",
+        )
+
+        np.testing.assert_array_equal(sm1.site_array["sids"], np.arange(9))
+        np.testing.assert_array_almost_equal(
+            sm1.site_array["lon"],
+            np.array([0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.0, 1.0, 2.0]),
+        )
+        np.testing.assert_array_almost_equal(
+            sm1.site_array["lat"],
+            np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0]),
+        )
+        nsites = 9
+        key_set = [
+            "depth",
+            "vs30",
+            "vs30measured",
+            "z1pt0",
+            "z2pt5",
+            "xvf",
+            "backarc",
+            "region",
+            "geology",
+        ]
+        value_set = [0.0, 1000.0, True, 0.0, 1.0, 100.0, False, 1, "ABC"]
+        for key, value in zip(key_set, value_set):
+            self._compare_constant_vectors(
+                sm1.site_array[key], value, nsites, SITE_PROPERTIES[key]
+            )
+
+    def test_build_from_bbox_minimal(self):
+        # Build the site model from a minimal set of inputs
+        sm1 = SiteModel.from_bbox(self.bbox, self.spcx, self.spcy, self.vs30)
+        np.testing.assert_array_equal(sm1.site_array["sids"], np.arange(9))
+        np.testing.assert_array_almost_equal(
+            sm1.site_array["lon"],
+            np.array([0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 0.0, 1.0, 2.0]),
+        )
+        np.testing.assert_array_almost_equal(
+            sm1.site_array["lat"],
+            np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0]),
+        )
+        nsites = 9
+        target_z1pt0 = SiteModel.vs30_to_z1pt0_cy14(np.array([self.vs30]))[0]
+        target_z2pt5 = SiteModel.vs30_to_z2pt5_cb14(np.array([self.vs30]))[0]
+        key_set = [
+            "depth",
+            "vs30",
+            "vs30measured",
+            "z1pt0",
+            "z2pt5",
+            "xvf",
+            "backarc",
+            "region",
+            "geology",
+        ]
+        value_set = [
+            0.0,
+            1000.0,
+            True,
+            target_z1pt0,
+            target_z2pt5,
+            150.0,
+            False,
+            0,
+            "UNCLASSIFIED",
+        ]
+        for key, value in zip(key_set, value_set):
+            self._compare_constant_vectors(
+                sm1.site_array[key], value, nsites, SITE_PROPERTIES[key]
+            )
+
+    def _compare_site_arrays(self, site_array1, site_array2):
+        # Compares the site arrays within two site collections to verify equality
+        # The name order is not fixed in OpenQuake, so need to check elements one-by-one
+        for name in site_array1.dtype.names:
+            if site_array1[name].dtype not in (np.float,):
+                self.assertTrue(np.array_equal(site_array1[name], site_array2[name]))
+            else:
+                np.testing.assert_array_almost_equal(
+                    site_array1[name], site_array2[name]
+                )
+
+    def test_build_site_collection(self):
+        # Checks that the returned site collection is accurately reproduced in
+        # the same OpenQuake site collection object
+        sm1 = SiteModel.from_bbox(
+            self.bbox,
+            self.spcx,
+            self.spcy,
+            vs30=800.0,
+            z1pt0=10.0,
+            z2pt5=1.0,
+            xvf=0.0,
+            backarc=False,
+            region=0,
+            geology="UNCLASSIFIED",
+        )
+        sm_site_collection = sm1.get_site_collection()
+        # Build model directly into OpenQuake SiteCollection
+        oq_sites = []
+        for i in np.arange(0.0, 3.0, 1.0):
+            for j in np.arange(0.0, 3.0, 1.0):
+                oq_sites.append(
+                    Site(
+                        Point(j, i),
+                        vs30=800.0,
+                        vs30measured=True,
+                        z1pt0=10.0,
+                        z2pt5=1.0,
+                        xvf=0.0,
+                        backarc=False,
+                        region=0,
+                        geology="UNCLASSIFIED",
+                    )
+                )
+        oq_site_collection = SiteCollection(oq_sites)
+        self.assertIsInstance(sm_site_collection, SiteCollection)
+        self._compare_site_arrays(sm_site_collection.array, oq_site_collection.array)
+
+    def test_build_site_collecion_filter(self):
+        # Checks that the site collection is filtered appropriately
+        # Requires only a minimal model right now
+        sm1 = SiteModel.from_bbox(self.bbox, self.spcx, self.spcy, self.vs30)
+        sm_site_collection = sm1.get_site_collection(np.array([2, 5, 7]))
+        # As all other attributed are constant, only with checking that the
+        # sids, lons and lats are correct
+        np.testing.assert_array_equal(sm_site_collection["sids"], np.array([2, 5, 7]))
+        np.testing.assert_array_almost_equal(
+            sm_site_collection["lon"], np.array([2.0, 2.0, 1.0])
+        )
+        np.testing.assert_array_almost_equal(
+            sm_site_collection["lat"], np.array([0.0, 1.0, 2.0])
+        )
+
+    def test_dataframe_input_output_full(self):
+        # Test correct export to and import from pandas dataframe
+        sm1 = SiteModel.from_bbox(self.bbox, self.spcx, self.spcy, self.vs30)
+        dframe1 = sm1.dataframe
+        # Create new site model from dataframe
+        sm2 = SiteModel.from_dataframe(dframe1)
+        # Verify that the site information is the same
+        self._compare_site_arrays(sm1.site_array, sm2.site_array)
+        # Now add on a new unrecognised field and verify that the correct warning is raised
+        dframe1["xyz"] = -99.0 * np.ones(len(sm1))
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter("always")
+            SiteModel.from_dataframe(dframe1)
+        self.assertEqual(
+            str(w[-1].message),
+            "Input site property xyz not recognised for use - skipping",
+        )
+
+    def test_dataframe_input_minimal(self):
+        # Test that the site model is instantated correctly from the dataframe when
+        # default values are required
+        sm1 = SiteModel.from_bbox(self.bbox, self.spcx, self.spcy, self.vs30)
+        dframe1 = sm1.dataframe
+        # Delete all non-essential columns from the dataframe
+        dframe1.drop(
+            ["vs30measured", "z1pt0", "z2pt5", "xvf", "backarc", "region", "geology"],
+            axis=1,
+            inplace=True,
+        )
+        # Instantiate new site model from reduced dataframe and verify that the
+        # resulting site array is the same as that generated from the bounding box
+        sm2 = SiteModel.from_dataframe(dframe1)
+        self._compare_site_arrays(sm1.site_array, sm2.site_array)