diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py index a15a9845f19461c9df221b8fa2d8bd1294e46353..f274211655b3e05a360fe3e70c6797e84297b2a2 100644 --- a/shakyground2/shakemap.py +++ b/shakyground2/shakemap.py @@ -10,6 +10,7 @@ import geopandas as gpd import matplotlib.pyplot as plt from typing import Dict, Optional, Tuple, List, Union from shapely.geometry import LineString, MultiLineString +from rasterio.warp import reproject from openquake.hazardlib import const, imt from openquake.hazardlib.contexts import ContextMaker from shakyground2 import valid @@ -357,34 +358,49 @@ class Shakemap(object): * rasterio.transform.Affine.scale(spcx, spcy) ) # Export to file + kwargs0 = { + "height": results.shape[0], + "width": results.shape[1], + "count": 1, + "dtype": results.dtype, + "transform": transform, + "crs": "+proj=latlong", + "compress": "lzw", + } if filename: - with rasterio.open( - filename, - "w", - "GTiff", - height=results.shape[0], - width=results.shape[1], - count=1, - dtype=results.dtype, - crs="+proj=latlong", - transform=transform, - compress="lzw", - ) as dst: - dst.write(results, 1) + with rasterio.open(filename, "w", "GTiff", **kwargs0) as dst: + if self.site_model.cross_antimeridian: + reproj_results = np.zeros(results.shape) + reproject( + results, + reproj_results, + src_transform=kwargs0["transform"], + src_crs=kwargs0["crs"], + dst_transform=kwargs0["transform"], + dst_crs=kwargs0["crs"], + kwargs={"CENTER_LONG": 180}, + ) + dst.write(reproj_results, 1) + else: + dst.write(results, 1) return # Returns a bytes object with rasterio.MemoryFile(ext=".tif") as memfile: - with memfile.open( - driver="GTiff", - height=results.shape[0], - width=results.shape[1], - count=1, - dtype=results.dtype, - crs="+proj=latlong", - transform=transform, - compress="lzw", - ) as dst: - dst.write(results, 1) + with memfile.open(driver="GTiff", **kwargs0) as dst: + if self.site_model.cross_antimeridian: + reproj_results = np.zeros(results.shape) + reproject( + results, + reproj_results, + src_transform=kwargs0["transform"], + src_crs=kwargs0["crs"], + dst_transform=kwargs0["transform"], + dst_crs=kwargs0["crs"], + kwargs={"CENTER_LONG": 180}, + ) + dst.write(reproj_results, 1) + else: + dst.write(results, 1) buffer = memfile.getbuffer() io_buffer = io.BytesIO(buffer) @@ -460,6 +476,8 @@ class Shakemap(object): self.site_model.bbox_properties["ncol"], ] lons = self.site_model["lon"].reshape(new_shape) + if self.site_model.cross_antimeridian: + lons[lons < 0] += 360.0 lats = self.site_model["lat"].reshape(new_shape) # Surpress plotting plt.ioff() diff --git a/shakyground2/site_model.py b/shakyground2/site_model.py index 7ce9daa43ba88af80a6d5a706b08d79be362c11d..3249e2a62c2cc2e45568d77cd057fd723b828514 100644 --- a/shakyground2/site_model.py +++ b/shakyground2/site_model.py @@ -44,12 +44,37 @@ class SiteModel(object): 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) + bbox_properties: Properties of the bounding box used to create the site array, if it + has been constructed from a bounding box. If the site model was not + constructed from a bounding box then "bbox_properties" is not needed + except if the site model spans more than one hemisphere without + crossing the meridian, in which case this should contain + "cross_antimeridian: False". """ - def __init__(self, site_array: np.ndarray, bbox_properties: Dict = {}): + def __init__(self, site_array: np.ndarray, bbox_properties: Optional[Dict] = None): + if bbox_properties is None: + bbox_properties = {} self.site_array = site_array + if "cross_antimeridian" in bbox_properties: + self.cross_antimeridian = bbox_properties.pop("cross_antimeridian") + else: + bbox = [ + np.min(site_array["lon"]), + np.min(site_array["lat"]), + np.max(site_array["lon"]), + np.max(site_array["lat"]), + ] + if (bbox[2] - bbox[0]) > 180.0: + warnings.warn( + "Upper longitude of bounding box exceeds lower longitude " + "by more than 180.0. Assuming antimeridian crossing. If this is " + "not the case set `bbox_properties={'cross_antimeridian': False}" + " in the call to the invoking function" + ) + self.cross_antimeridian = True + else: + self.cross_antimeridian = False self.bbox_properties = bbox_properties def __len__(self): @@ -155,9 +180,19 @@ class SiteModel(object): # 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) - ) + if (llon > ulon) and (np.fabs(llon - ulon) > 180.0): + # Bounding box crosses the antimeridian + lons, lats = np.meshgrid( + np.arange(llon, ulon + 360.0 + spcx, spcx), np.arange(llat, ulat + spcy, spcy) + ) + cross_antimeridian = True + lons[lons > 180.0] -= 360.0 + else: + + lons, lats = np.meshgrid( + np.arange(llon, ulon + spcx, spcx), np.arange(llat, ulat + spcy, spcy) + ) + cross_antimeridian = False nrow, ncol = lons.shape nsites = nrow * ncol self = object.__new__(cls) @@ -169,6 +204,7 @@ class SiteModel(object): "spcy": spcy, "ncol": ncol, "nrow": nrow, + "cross_antimeridian": cross_antimeridian, } # Site array has a set of defined datatypes site_dtype = np.dtype(list(SITE_PROPERTIES.items())) @@ -200,7 +236,9 @@ class SiteModel(object): return cls(site_array, bbox_properties=bbox_properties) @classmethod - def from_dataframe(cls, dframe: pd.DataFrame) -> SiteModel: + def from_dataframe( + cls, dframe: pd.DataFrame, bbox_properties: Optional[Dict] = None + ) -> SiteModel: """ Builds the site collection directly from a pandas dataframe @@ -210,6 +248,8 @@ class SiteModel(object): Returns: Instance of :class`shakygroundv2.site_model.SiteModel` """ + if bbox_properties is None: + bbox_properties = {} nsites = dframe.shape[0] site_dtype = np.dtype(list(SITE_PROPERTIES.items())) site_array = np.zeros(nsites, site_dtype) @@ -235,7 +275,7 @@ class SiteModel(object): 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) + return cls(site_array, bbox_properties) @staticmethod def _get_string_array(num_elem: int, key: str, dtype) -> np.ndarray: diff --git a/shakyground2/synthetic_rupture_generator.py b/shakyground2/synthetic_rupture_generator.py index f91448145b19b86de0a5e96e5f8bed279ba4b1bf..8081e99a3f15b41f4a73e56bc8cd9c06bbf10647 100644 --- a/shakyground2/synthetic_rupture_generator.py +++ b/shakyground2/synthetic_rupture_generator.py @@ -117,10 +117,19 @@ class FiniteRuptureSampler(object): # the hypocentre spaced every 0.05 degrees - longer distances are less likely # to be significantly influenced by uncertainties in the rupture plane llon, llat, ulon, ulat = earthquake.get_maximum_distance_bbox(maximum_site_distance) - target_lons, target_lats = np.meshgrid( - np.arange(llon, ulon + site_spacing, site_spacing), - np.arange(llat, ulat + site_spacing, site_spacing), - ) + if ulon < llon: + # Happens when crosses the antimeridian + target_lons, target_lats = np.meshgrid( + np.arange(llon, ulon + 360.0 + site_spacing, site_spacing), + np.arange(llat, ulat + site_spacing, site_spacing), + ) + target_lons[target_lons > 180.0] -= 360.0 + else: + target_lons, target_lats = np.meshgrid( + np.arange(llon, ulon + site_spacing, site_spacing), + np.arange(llat, ulat + site_spacing, site_spacing), + ) + target_lons[target_lons > 180.0] -= 360.0 target_lons = target_lons.flatten() target_lats = target_lats.flatten() diff --git a/shakyground2/workflows.py b/shakyground2/workflows.py index 849fa0b7937e8b918b2d9d82e4a1beaf169eeee3..c1057acfbfa7ae080600a715ad530bb18dd2ad37 100644 --- a/shakyground2/workflows.py +++ b/shakyground2/workflows.py @@ -119,6 +119,7 @@ def shakemaps_from_quakeml( ) results.num_sites = len(site_model) results.bbox_properties = site_model.bbox_properties + results.bbox_properties["cross_antimeridian"] = site_model.cross_antimeridian # Get the ground motion models results.tectonic_region, results.ground_motion_models = regionalization(results.earthquake) shakemap_config = {} diff --git a/tests/data/gfz2021jqav.xml b/tests/data/gfz2021jqav.xml new file mode 100644 index 0000000000000000000000000000000000000000..68ea92a2e5beb5bb7a38590176544ab84b4ed37f --- /dev/null +++ b/tests/data/gfz2021jqav.xml @@ -0,0 +1 @@ +Rat Islands, Aleutian Islandsregion nameGFZ2021-05-17T20:33:32.797114Z1.519124002e+171.150862387e+17-1.237443836e+178.658144886e+158.143546196e+162.588668599e+16-3.820056546e+160.81692278630.7658606440.234139356GFZ2021-05-18T06:05:05.574667Zsmi:org.gfz-potsdam.de/geofon/Origin/20210518060505.575388.858842smi:org.gfz-potsdam.de/geofon/Magnitude/20210518060505.57545.858843smi:org.gfz-potsdam.de/geofon/BP_60s-150ssmi:org.gfz-potsdam.de/geofon/GEOFON_standard73.3066743162.2365899188.23008096257.101886827.8152913393.35776884339.093618172.705782961.416483959e+17164.616824617.21868289-1.60429869e+1774.13134631.5661094861.878147311e+1683.762354420.1830772137GFZ2021-05-18T06:05:05.574667Zsmi:org.gfz-potsdam.de/geofon/Origin/20210518055932.117269.1236929179.4231111.91611528451.320301063.84954094923117922417270.8228094851160.013641490.358741762.49188995470.25522614manualGFZ2021-05-18T05:59:32.122727Z20429.89541374.8097424068.542488249.555588176.0296631horizontal uncertaintysmi:org.gfz-potsdam.de/geofon/LOCSATsmi:org.gfz-potsdam.de/geofon/iasp91confirmed4383.76235442GFZ2021-05-18T06:05:05.574667Z5.387728817Mwsmi:org.gfz-potsdam.de/geofon/MTsmi:org.gfz-potsdam.de/geofon/Origin/20210518055932.117269.1236929smi:org.gfz-potsdam.de/geofon/Magnitude/20210518060505.57545.858843smi:org.gfz-potsdam.de/geofon/FocalMechanism/20210518060505.575219.858841 diff --git a/tests/site_model_test.py b/tests/site_model_test.py index 9e7869d3529009489063b922846b6ee88c4f7fe5..555d17d25b325566e91b125a078880355621369a 100644 --- a/tests/site_model_test.py +++ b/tests/site_model_test.py @@ -4,6 +4,7 @@ Test suite for the site model class import unittest import warnings import numpy as np +import pandas as pd from openquake.hazardlib.geo import Point from openquake.hazardlib.site import Site, SiteCollection from shakyground2.site_model import SiteModel, SITE_PROPERTIES @@ -20,16 +21,21 @@ class SiteModelTestCase(unittest.TestCase): def test_length(self): dummy_array = np.zeros([20, 10]) - sc1 = SiteModel(dummy_array) + sc1 = SiteModel(dummy_array, bbox_properties={"cross_antimeridian": False}) self.assertEqual(len(sc1), 20) def test_string(self): # No bounding box case dummy_array = np.zeros([20, 10]) - sc1 = SiteModel(dummy_array) + sc1 = SiteModel(dummy_array, bbox_properties={"cross_antimeridian": False}) 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} + bbox_properties = { + "bbox": [1.0, 2.0, 3.0, 4.0], + "spcx": 0.5, + "spcy": 0.5, + "cross_antimeridian": False, + } sc1 = SiteModel(dummy_array, bbox_properties) self.assertEqual( str(sc1), @@ -128,6 +134,81 @@ class SiteModelTestCase(unittest.TestCase): sm1.site_array[key], value, nsites, SITE_PROPERTIES[key] ) + def test_build_from_bbox_crossing_antimeridian(self): + # Check that a valid site model is built when the bbox crosses the antimeridian + sm1 = SiteModel.from_bbox([178.0, -30.0, -178.0, -28.0], 1.0, 1.0, 800) + self.assertEqual(len(sm1), 15) + self.assertTrue(sm1.cross_antimeridian) + lons = np.array( + [ + 178.0, + 179.0, + 180.0, + -179.0, + -178.0, + 178.0, + 179.0, + 180.0, + -179.0, + -178.0, + 178.0, + 179.0, + 180.0, + -179.0, + -178.0, + ] + ) + lats = np.array( + [ + -30.0, + -30.0, + -30.0, + -30.0, + -30.0, + -29.0, + -29.0, + -29.0, + -29.0, + -29.0, + -28.0, + -28.0, + -28.0, + -28.0, + -28.0, + ] + ) + np.testing.assert_array_almost_equal(sm1["lon"], lons) + np.testing.assert_array_almost_equal(sm1["lat"], lats) + + def test_build_from_dataframe_crossing_antimeridian(self): + # Test case for building site model from a dataframe with sites + # crossing the antimeridian + lons, lats = np.meshgrid(np.arange(178.0, 182.5, 0.5), np.arange(-30.0, -25.0, 0.5)) + lons = lons.flatten() + lons[lons > 180.0] -= 360.0 + lats = lats.flatten() + sids = np.arange(1, len(lons) + 1) + vs30 = 800.0 * np.ones(len(lons)) + input_site_data = pd.DataFrame({"sids": sids, "lon": lons, "lat": lats, "vs30": vs30}) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + sm1 = SiteModel.from_dataframe(input_site_data) + self.assertEqual( + str(w[-1].message), + "Upper longitude of bounding box exceeds lower longitude " + "by more than 180.0. Assuming antimeridian crossing. If this is " + "not the case set `bbox_properties={'cross_antimeridian': False}" + " in the call to the invoking function", + ) + self.assertTrue(sm1.cross_antimeridian) + + # Check the case when the crossing antimeridian check is over-ridden + sm2 = SiteModel.from_dataframe( + input_site_data, bbox_properties={"cross_antimeridian": False} + ) + self.assertFalse(sm2.cross_antimeridian) + self.assertEqual(len(sm1), len(sm2)) + 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) diff --git a/tests/workflows_test.py b/tests/workflows_test.py index 853cf90ad7200204456812717d3013455bf766db..7c7567d69212e382b6f9a2a4b577da656381499a 100644 --- a/tests/workflows_test.py +++ b/tests/workflows_test.py @@ -46,6 +46,18 @@ class GeofonWorkflowTestCase(unittest.TestCase): # Verify that none of the statistics are nans self.assertFalse(np.isnan(results.statistics["mean"][imt][stat_type])) + def test_complete_workflow_crossing_antimeridian(self): + event_id = os.path.join(DATA_PATH, "gfz2021jqav.xml") + results = shakemaps_from_quakeml(event_id, self.imts, config={"spcx": 0.1, "spcy": 0.1}) + self.assertEqual(results.num_sites, 14841) + self.assertListEqual(list(results.shakemaps), ["mean", "stddevs"]) + self.assertListEqual(list(results.shakemaps["mean"]), self.imts) + self.assertIsInstance(results, ShakemapWorkflowResult) + for imt in self.imts: + for maptype in results.shakemaps: + self.assertIsInstance(results.shakemaps[maptype][imt], bytes) + self.assertIsInstance(results.contours["mean"][imt], GeoDataFrame) + def test_workflow_with_exports(self): test_folder = os.path.join(DATA_PATH, "tmp_shakemaps_from_quakeml") event_id = os.path.join(DATA_PATH, "gfz2021eksc.xml")