Commit 8487414b authored by Graeme Weatherill's avatar Graeme Weatherill
Browse files

(Potentially) Fixes bug of shakemaps not working across antimeridian

parent cebcb8e2
Pipeline #23891 passed with stage
in 12 minutes and 17 seconds
......@@ -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()
......
......@@ -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:
......
......@@ -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()
......
......@@ -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 = {}
......
<?xml version="1.0" encoding="UTF-8"?><q:quakeml xmlns="http://quakeml.org/xmlns/bed/1.2" xmlns:q="http://quakeml.org/xmlns/quakeml/1.2"><eventParameters publicID="smi:org.gfz-potsdam.de/geofon/EventParameters"><event publicID="smi:org.gfz-potsdam.de/geofon/gfz2021jqav"><description><text>Rat Islands, Aleutian Islands</text><type>region name</type></description><creationInfo><agencyID>GFZ</agencyID><creationTime>2021-05-17T20:33:32.797114Z</creationTime></creationInfo><focalMechanism publicID="smi:org.gfz-potsdam.de/geofon/FocalMechanism/20210518060505.575219.858841"><momentTensor publicID="smi:org.gfz-potsdam.de/geofon/MomentTensor/20210518060505.575518.858844"><scalarMoment><value>1.519124002e+17</value></scalarMoment><tensor><Mrr><value>1.150862387e+17</value></Mrr><Mtt><value>-1.237443836e+17</value></Mtt><Mpp><value>8.658144886e+15</value></Mpp><Mrt><value>8.143546196e+16</value></Mrt><Mrp><value>2.588668599e+16</value></Mrp><Mtp><value>-3.820056546e+16</value></Mtp></tensor><varianceReduction>0.8169227863</varianceReduction><doubleCouple>0.765860644</doubleCouple><clvd>0.234139356</clvd><creationInfo><agencyID>GFZ</agencyID><creationTime>2021-05-18T06:05:05.574667Z</creationTime></creationInfo><derivedOriginID>smi:org.gfz-potsdam.de/geofon/Origin/20210518060505.575388.858842</derivedOriginID><momentMagnitudeID>smi:org.gfz-potsdam.de/geofon/Magnitude/20210518060505.57545.858843</momentMagnitudeID><filterID>smi:org.gfz-potsdam.de/geofon/BP_60s-150s</filterID><methodID>smi:org.gfz-potsdam.de/geofon/GEOFON_standard</methodID></momentTensor><nodalPlanes><nodalPlane1><strike><value>73.30667431</value></strike><dip><value>62.23658991</value></dip><rake><value>88.23008096</value></rake></nodalPlane1><nodalPlane2><strike><value>257.1018868</value></strike><dip><value>27.81529133</value></dip><rake><value>93.35776884</value></rake></nodalPlane2></nodalPlanes><principalAxes><tAxis><azimuth><value>339.0936181</value></azimuth><plunge><value>72.70578296</value></plunge><length><value>1.416483959e+17</value></length></tAxis><pAxis><azimuth><value>164.6168246</value></azimuth><plunge><value>17.21868289</value></plunge><length><value>-1.60429869e+17</value></length></pAxis><nAxis><azimuth><value>74.1313463</value></azimuth><plunge><value>1.566109486</value></plunge><length><value>1.878147311e+16</value></length></nAxis></principalAxes><azimuthalGap>83.76235442</azimuthalGap><misfit>0.1830772137</misfit><creationInfo><agencyID>GFZ</agencyID><creationTime>2021-05-18T06:05:05.574667Z</creationTime></creationInfo><triggeringOriginID>smi:org.gfz-potsdam.de/geofon/Origin/20210518055932.117269.1236929</triggeringOriginID></focalMechanism><origin publicID="smi:org.gfz-potsdam.de/geofon/Origin/20210518055932.117269.1236929"><time><value>2021-05-17T20:25:07.717779Z</value><uncertainty>0.2312954217</uncertainty></time><longitude><value>179.423111</value><uncertainty>1.916115284</uncertainty></longitude><latitude><value>51.32030106</value><uncertainty>3.849540949</uncertainty></latitude><quality><associatedPhaseCount>231</associatedPhaseCount><usedPhaseCount>179</usedPhaseCount><associatedStationCount>224</associatedStationCount><usedStationCount>172</usedStationCount><depthPhaseCount>7</depthPhaseCount><standardError>0.8228094851</standardError><azimuthalGap>160.0136414</azimuthalGap><maximumDistance>90.35874176</maximumDistance><minimumDistance>2.491889954</minimumDistance><medianDistance>70.25522614</medianDistance></quality><evaluationMode>manual</evaluationMode><creationInfo><agencyID>GFZ</agencyID><creationTime>2021-05-18T05:59:32.122727Z</creationTime></creationInfo><depth><value>20429.8954</value><uncertainty>1374.809742</uncertainty></depth><originUncertainty><minHorizontalUncertainty>4068.54248</minHorizontalUncertainty><maxHorizontalUncertainty>8249.555588</maxHorizontalUncertainty><azimuthMaxHorizontalUncertainty>176.0296631</azimuthMaxHorizontalUncertainty><preferredDescription>horizontal uncertainty</preferredDescription></originUncertainty><methodID>smi:org.gfz-potsdam.de/geofon/LOCSAT</methodID><earthModelID>smi:org.gfz-potsdam.de/geofon/iasp91</earthModelID><evaluationStatus>confirmed</evaluationStatus></origin><magnitude publicID="smi:org.gfz-potsdam.de/geofon/Magnitude/20210518060505.57545.858843"><stationCount>43</stationCount><azimuthalGap>83.76235442</azimuthalGap><creationInfo><agencyID>GFZ</agencyID><creationTime>2021-05-18T06:05:05.574667Z</creationTime></creationInfo><mag><value>5.387728817</value></mag><type>Mw</type><methodID>smi:org.gfz-potsdam.de/geofon/MT</methodID></magnitude><preferredOriginID>smi:org.gfz-potsdam.de/geofon/Origin/20210518055932.117269.1236929</preferredOriginID><preferredMagnitudeID>smi:org.gfz-potsdam.de/geofon/Magnitude/20210518060505.57545.858843</preferredMagnitudeID><preferredFocalMechanismID>smi:org.gfz-potsdam.de/geofon/FocalMechanism/20210518060505.575219.858841</preferredFocalMechanismID></event></eventParameters></q:quakeml>
......@@ -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)
......
......@@ -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")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment