diff --git a/shakyground2/site_model.py b/shakyground2/site_model.py new file mode 100644 index 0000000000000000000000000000000000000000..5fba99da74a1aa074ee7d6af53fd7038179e90e0 --- /dev/null +++ b/shakyground2/site_model.py @@ -0,0 +1,306 @@ +""" +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)) diff --git a/tests/test_site_model.py b/tests/test_site_model.py new file mode 100644 index 0000000000000000000000000000000000000000..490338481a70a09a63db00c6723b6ba1ba363c1d --- /dev/null +++ b/tests/test_site_model.py @@ -0,0 +1,268 @@ +""" +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)