Commit 6a3991d7 authored by Graeme Weatherill's avatar Graeme Weatherill
Browse files

Merge branch 'pevens_import_fdsnws' into 'import_fdsnws'

Cleanup and fixes to fdsnws importer

See merge request pevans/shakyground2!1
parents 5c5db22a cfc499e3
Pipeline #22389 passed with stage
in 8 minutes and 48 seconds
...@@ -15,6 +15,7 @@ setup( ...@@ -15,6 +15,7 @@ setup(
"geopandas", "geopandas",
"pyyaml", "pyyaml",
"Rtree", "Rtree",
"rasterio",
], ],
packages=find_packages(), packages=find_packages(),
python_requires=">=3.7", python_requires=">=3.7",
......
...@@ -17,10 +17,9 @@ Can run standalone: ...@@ -17,10 +17,9 @@ Can run standalone:
""" """
import datetime
import os import os
import sys import sys
from typing import List, Optional, Union from typing import List, Union
import urllib.request as ul import urllib.request as ul
from xml.etree import ElementTree as ET from xml.etree import ElementTree as ET
import yaml import yaml
...@@ -63,15 +62,16 @@ def fetch_fm(root: ET.Element, ns: dict, preferredfmID: str) -> list: ...@@ -63,15 +62,16 @@ def fetch_fm(root: ET.Element, ns: dict, preferredfmID: str) -> list:
d: List[dict] = [dict(), dict()] d: List[dict] = [dict(), dict()]
for k in range(2): for k in range(2):
plane = np.find("ns:nodalPlane" + str(k + 1), ns) plane = np.find("ns:nodalPlane" + str(k + 1), ns)
if plane == None: if plane is None:
continue continue
for child in plane: for child in plane:
found = child.find("ns:value", ns) found = child.find("ns:value", ns)
if found == None: if found is None:
continue continue
v = found.text v = found.text
tag = child.tag tag = child.tag
tag = tag[tag.index("}") + 1 :] tag_idx = tag.index("}") + 1
tag = tag[tag_idx:]
try: try:
d[k][tag] = float(str(v)) d[k][tag] = float(str(v))
except ValueError: except ValueError:
...@@ -87,10 +87,10 @@ def fetch_magnitude( ...@@ -87,10 +87,10 @@ def fetch_magnitude(
continue continue
child = m.find("ns:mag", ns) child = m.find("ns:mag", ns)
if child == None: if child is None:
return None return None
value = child.find("ns:value", ns) value = child.find("ns:value", ns)
if value == None: if value is None:
return None return None
v = value.text v = value.text
try: try:
...@@ -115,7 +115,8 @@ def fetch_origin(root, ns: dict, preferredoriginID: str) -> dict: ...@@ -115,7 +115,8 @@ def fetch_origin(root, ns: dict, preferredoriginID: str) -> dict:
d = dict() d = dict()
for child in o: for child in o:
tag = child.tag tag = child.tag
tag = tag[tag.index("}") + 1 :] tag_idx = tag.index("}") + 1
tag = tag[tag_idx:]
if tag in ("depth", "latitude", "longitude", "time"): if tag in ("depth", "latitude", "longitude", "time"):
v = child.find("ns:value", ns).text v = child.find("ns:value", ns).text
...@@ -152,14 +153,11 @@ def fetch_quakeml_ws(evid: str) -> str: ...@@ -152,14 +153,11 @@ def fetch_quakeml_ws(evid: str) -> str:
) )
) )
# print(url)
req = ul.Request(url) req = ul.Request(url)
u = ul.urlopen(req) u = ul.urlopen(req)
buf = u.read().decode("utf8") buf = u.read().decode("utf8")
print("Got", len(buf), "char(s).") print("Got", len(buf), "char(s).")
with open("tmp.xml", "w") as fid:
fid.write(str(buf))
return buf return buf
...@@ -186,7 +184,6 @@ def fetch_quakeml(path: str) -> Union[dict, None]: ...@@ -186,7 +184,6 @@ def fetch_quakeml(path: str) -> Union[dict, None]:
ns = QUAKEML_NS ns = QUAKEML_NS
event = root[0][0] event = root[0][0]
# .find("ns:event", ns)
if not event: if not event:
print("Couldn't get an event!") print("Couldn't get an event!")
return None return None
...@@ -196,33 +193,26 @@ def fetch_quakeml(path: str) -> Union[dict, None]: ...@@ -196,33 +193,26 @@ def fetch_quakeml(path: str) -> Union[dict, None]:
try: try:
# e.g. "smi:org.gfz-potsdam.de/geofon/gfz2021ekhv" # e.g. "smi:org.gfz-potsdam.de/geofon/gfz2021ekhv"
evid = event.attrib["publicID"].split("/")[-1] evid = event.attrib["publicID"].split("/")[-1]
# print("INFO Event ID from response:", evid)
except AttributeError: except AttributeError:
_perror("Oops, couldn't get event id from " + event.attrib["publicID"]) _perror("Oops, couldn't get event id from " + event.attrib["publicID"])
# descr = root[0][0].find("ns:description", ns)
# text = descr.find("ns:text", ns)
# print("INFO Found description:", text.text)
preferredoriginIDElem = root[0][0].find("ns:preferredOriginID", ns) preferredoriginIDElem = root[0][0].find("ns:preferredOriginID", ns)
if preferredoriginIDElem == None: if preferredoriginIDElem is None:
_perror("Oops, couldn't find the preferred origin ID") _perror("Oops, couldn't find the preferred origin ID")
return None return None
preferredoriginID = preferredoriginIDElem.text preferredoriginID = preferredoriginIDElem.text
# print("DEBUG Found", preferredoriginID)
preferredmagID = None preferredmagID = None
preferredMagnitudeIDElem = root[0][0].find("ns:preferredMagnitudeID", ns) preferredMagnitudeIDElem = root[0][0].find("ns:preferredMagnitudeID", ns)
if preferredMagnitudeIDElem != None: if preferredMagnitudeIDElem is not None:
try: try:
preferredmagID = preferredMagnitudeIDElem.text preferredmagID = preferredMagnitudeIDElem.text
except AttributeError: except AttributeError:
pass pass
preferredfmID = None preferredfmID = None
try: try:
preferredfmIDElem = root[0][0].find("ns:preferredFocalMechanismID", ns) preferredfmIDElem = root[0][0].find("ns:preferredFocalMechanismID", ns)
if preferredfmIDElem != None: if preferredfmIDElem is not None:
preferredfmID = preferredfmIDElem.text preferredfmID = preferredfmIDElem.text
except AttributeError: except AttributeError:
pass pass
...@@ -231,7 +221,6 @@ def fetch_quakeml(path: str) -> Union[dict, None]: ...@@ -231,7 +221,6 @@ def fetch_quakeml(path: str) -> Union[dict, None]:
print("Oops, no preferredOriginID was found") print("Oops, no preferredOriginID was found")
return None return None
origin = fetch_origin(root, ns, preferredoriginID) origin = fetch_origin(root, ns, preferredoriginID)
(d, t) = origin.pop("time").split("T", 2) (d, t) = origin.pop("time").split("T", 2)
# There's little point in forcing these into datetimes, # There's little point in forcing these into datetimes,
# since writing to YAML requires they be converted back # since writing to YAML requires they be converted back
...@@ -283,7 +272,7 @@ http://geofon.gfz-potsdam.de/eqinfo/list.php ...@@ -283,7 +272,7 @@ http://geofon.gfz-potsdam.de/eqinfo/list.php
sys.exit(1) sys.exit(1)
ev_info = fetch_quakeml(evid) ev_info = fetch_quakeml(evid)
if ev_info == None: if ev_info is None:
print("Got no dictionary from QuakeML") print("Got no dictionary from QuakeML")
outfile = evid + ".yaml" outfile = evid + ".yaml"
with open(outfile, "w") as fid: with open(outfile, "w") as fid:
......
""" """
Implements the core shakemap class Implements the core shakemap class
""" """
import io
import h5py import h5py
import numpy as np import numpy as np
import warnings import warnings
from typing import Dict, Optional, Tuple, List import rasterio
import geopandas as gpd
import matplotlib.pyplot as plt
from typing import Dict, Optional, Tuple, List, Union
from shapely.geometry import LineString, MultiLineString
from openquake.hazardlib import const, imt from openquake.hazardlib import const, imt
from openquake.hazardlib.contexts import ContextMaker from openquake.hazardlib.contexts import ContextMaker
from shakyground2 import valid from shakyground2 import valid
...@@ -286,3 +291,229 @@ class Shakemap(object): ...@@ -286,3 +291,229 @@ class Shakemap(object):
"stddev", aggregated_stddevs.shape, shakemap_dtypes "stddev", aggregated_stddevs.shape, shakemap_dtypes
) )
aggregated_stddev_dset[:] = aggregated_stddevs aggregated_stddev_dset[:] = aggregated_stddevs
def _transform_to_raster(
self, imt: str, shakemap: np.ndarray, is_stddev: bool = False
) -> np.ndarray:
"""
Transforms a shakemap output into a 2D array for raster export based
on the bounding box properties in the site model
Args:
imt: Intensity measure type
shakemap: shakemap outputs from `get_shakemap(...)`
is_stddev: shakemap values are standard deviations (True) or expected ground motion
Returns:
shakemap as 2D numpy array
"""
if not self.site_model.bbox_properties:
raise ValueError(
"No bounding bbox properties found in SiteModel - "
"cannot export to raster format"
)
if is_stddev:
# If the input is a standard deviation then no transformation of the data
results = np.copy(shakemap[imt])
else:
# For absolute ground motion take the natural exponent
results = np.exp(shakemap[imt])
assert (
results.shape[0] == self.site_model.shape
), "Shakemap dimensions do not correspond to site model"
# Reshape and export to raster
shakemap_shape = [
self.site_model.bbox_properties["nrow"],
self.site_model.bbox_properties["ncol"],
]
return results.reshape(shakemap_shape)
def to_geotiff(
self,
shakemap: np.ndarray,
imt: str,
filename: Optional[str] = None,
is_stddev: bool = False,
) -> Union[io.BytesIO, None]:
"""
Exports the shakemap for a given intensity measure type to GeoTiff format
Args:
shakemap: shakemap outputs from `get_shakemap(...)`
imt: Intensity measure type
filename: Path to geotiff file for export. If None then returns the geotiff as a
bytes object
is_stddev: shakemap values are standard deviations (True) or expected ground motion
"""
results = self._transform_to_raster(imt, shakemap, is_stddev)
spcx = self.site_model.bbox_properties["spcx"]
spcy = self.site_model.bbox_properties["spcy"]
# Build the transformation object (holds the geospatial information about the raster)
transform = (
rasterio.transform.Affine.translation(
self.site_model.bbox_properties["bbox"][0] - (spcx / 2.0),
self.site_model.bbox_properties["bbox"][1] - (spcy / 2.0),
)
* rasterio.transform.Affine.scale(spcx, spcy)
)
# Export to file
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)
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)
buffer = memfile.getbuffer()
io_buffer = io.BytesIO(buffer)
as_bytes = io_buffer.read()
return as_bytes
def to_esri_ascii(
self,
shakemap: np.ndarray,
imt: str,
filename: str,
is_stddev: bool = False,
nodata: Optional[float] = None,
fmt: str = "%.4e",
):
"""
Exports the shakemap for a given intensity measure type to ESRI ascii format
Args:
shakemap: shakemap outputs from `get_shakemap(...)`
imt: Intensity measure type
filename: Path to ESRI Ascii file for export
is_stddev: shakemap values are standard deviations (True) or expected ground motion
"""
if not np.isclose(
self.site_model.bbox_properties["spcx"], self.site_model.bbox_properties["spcy"]
):
# ESRI Grid files cannot accommodate this
raise IOError(
"ESRI ascii files cannot accommodate unequal spacing "
"in the X and Y dimension - try GeoTIFF instead"
)
results = self._transform_to_raster(imt, shakemap, is_stddev)
# Export the raster to file
with open(filename, "w") as f:
# Write header information
f.write("NCOLS %g\n" % self.site_model.bbox_properties["ncol"])
f.write("NROWS %g\n" % self.site_model.bbox_properties["nrow"])
f.write("XLLCENTER %.8f\n" % self.site_model.bbox_properties["bbox"][0])
f.write("YLLCENTER %.8f\n" % self.site_model.bbox_properties["bbox"][1])
f.write("CELLSIZE %.8e\n" % self.site_model.bbox_properties["spcx"])
if nodata is not None:
f.write("NODATA_VALUE %.8f\n" % nodata)
# Write the data
np.savetxt(f, results, fmt=fmt)
return
def get_contours(
self,
imt: str,
shakemap: np.ndarray,
levels: Union[int, np.ndarray] = 10,
is_stddev: bool = False,
) -> gpd.GeoDataFrame:
"""
For a shakemap of a given intensity measure type, retrieves the
set of contours and exports to a geopandas GeoDataframe
Args:
imt: Intensity measure type
shakemap: shakemap outputs from `get_shakemap(...)`
levels: Number of levels for the contours or pre-defined set of levels (see
documentation for matplotlib.pyplot.contour)
is_stddev: shakemap values are standard deviations (True) or expected ground motion
Returns:
Contour set as geopandas GeoDataframe
"""
# For contouring we need the shakemap in 2D
results = self._transform_to_raster(imt, shakemap, is_stddev)
new_shape = [
self.site_model.bbox_properties["nrow"],
self.site_model.bbox_properties["ncol"],
]
lons = self.site_model["lon"].reshape(new_shape)
lats = self.site_model["lat"].reshape(new_shape)
# Surpress plotting
plt.ioff()
fig = plt.figure()
ax = fig.add_subplot(111)
# Get the contours
contours = ax.contour(lons, lats, results, levels)
contour_levels = []
contour_geometries = []
# Retrieve the x, y coordinates of each contour
for level, segs in zip(contours.levels, contours.allsegs):
if not len(segs):
continue
else:
geometries = []
for seg in segs:
geometries.append(LineString([(row[0], row[1]) for row in seg]))
if len(geometries) > 1:
# Has multiple contour lines for this level
contour_geometries.append(MultiLineString(geometries))
else:
# Only a single contour for this level
contour_geometries.append(geometries[0])
contour_levels.append(level)
plt.close(fig)
# Re-allow plotting
plt.ion()
# Create geodataframe
dframe = gpd.GeoDataFrame(
{imt: contour_levels, "geometry": gpd.GeoSeries(contour_geometries)}
)
dframe.crs = {"init": "epsg:4326"}
return dframe
def contours_to_file(
self,
shakemap: np.ndarray,
imt: str,
filename: str,
levels: Union[int, np.ndarray] = 10,
is_stddev: bool = False,
driver: str = "GeoJSON",
):
"""
Exports the contours of a shapefile to a vector spatial format
Args:
imt: Intensity measure type
shakemap: shakemap outputs from `get_shakemap(...)`
levels: Number of levels for the contours or pre-defined set of levels (see
documentation for matplotlib.pyplot.contour)
is_stddev: shakemap values are standard deviations (True) or expected ground motion
driver: Any Vector spatial file driver supported by Fiona
https://github.com/Toblerity/Fiona/blob/master/fiona/drvsupport.py
"""
contours = self.get_contours(imt, shakemap, levels, is_stddev)
contours.to_file(filename, driver=driver)
return
This diff is collapsed.
This diff is collapsed.
...@@ -4,7 +4,9 @@ Set of test cases for the shakemap class ...@@ -4,7 +4,9 @@ Set of test cases for the shakemap class
import os import os
import unittest import unittest
import h5py import h5py
import rasterio
import numpy as np import numpy as np
from geopandas import GeoDataFrame
from shakyground2.earthquake import Earthquake from shakyground2.earthquake import Earthquake
from shakyground2.site_model import SiteModel from shakyground2.site_model import SiteModel
from shakyground2.shakemap import Shakemap from shakyground2.shakemap import Shakemap
...@@ -105,3 +107,191 @@ class ShakemapTestCase(unittest.TestCase): ...@@ -105,3 +107,191 @@ class ShakemapTestCase(unittest.TestCase):
def tearDown(self): def tearDown(self):
# Delete the cache file # Delete the cache file
os.remove(self.cache_test_file) os.remove(self.cache_test_file)
class ShakemapExportersTestCase(unittest.TestCase):
"""
General test case for the shakemap exporters - synthetic 2D shakemap based on the function
matplotlib contour demo (https://tinyurl.com/22sdmajk)
"""
def setUp(self):
self.site_model = SiteModel.from_bbox(
[-3.0, -2.0, 3.0, 2.0], 0.025, 0.025, vs30=800.0, z1pt0=10.0
)
self.shakemap = Shakemap(
Earthquake("XYZ", 0.0, 0.0, 10.0, 6.0),
self.site_model,
[("B14", GSIM_LIST["BindiEtAl2014Rjb"](), 1.0)],
"ASC",
)
# 2D smooth function for contouring
z1 = np.exp(-self.site_model["lon"] ** 2 - self.site_model["lat"] ** 2)
z2 = np.exp(-((self.site_model["lon"] - 1) ** 2) - (self.site_model["lat"] - 1) ** 2)
self.dummy = (z1 - z2) * 2
def test_transform_to_raster(self):
# Test transformation to raster
# Map the 2D function (double peaked) into a shakemap output format
results = np.zeros(self.site_model.shape, dtype=np.dtype([("PGA", np.float64)]))
results["PGA"] = np.log(np.fabs(self.dummy) + 0.001)
results_as_raster = self.shakemap._transform_to_raster("PGA", results)
# Verify that this gives the expected results
delta = 0.025
x = np.arange(-3.0, 3.0 + delta, delta)
y = np.arange(-2.0, 2.0 + delta, delta)
X, Y = np.meshgrid(x, y)
Z1 = np.exp(-(X ** 2.0) - Y ** 2)
Z2 = np.exp(-((X - 1) ** 2) - (Y - 1) ** 2)
Z = np.fabs((Z1 - Z2) * 2) + 0.001
np.testing.assert_array_almost_equal(results_as_raster, Z)
def test_incorrect_raster_formulations(self):
# Tests the cases where the raster cannot be built
# Case 1 - No bbox properties in Site Model
site_model = SiteModel.from_bbox(
[-3.0, -2.0, 3.0, 2.0], 0.025, 0.025, vs30=800.0, z1pt0=10.0
)
# Re-set bounding box properties
site_model.bbox_properties = {}
shakemap = Shakemap(
Earthquake("XYZ", 0.0, 0.0, 10.0, 6.0),
site_model,
[("B14", GSIM_LIST["BindiEtAl2014Rjb"](), 1.0)],
"ASC",
)
results = np.zeros(site_model.shape, dtype=np.dtype([("PGA", np.float64)]))
results["PGA"] = np.log(np.fabs(self.dummy) + 0.001)
with self.assertRaises(ValueError) as ve:
shakemap._transform_to_raster("PGA", results)
self.assertEqual(
str(ve.exception),
"No bounding bbox properties found in SiteModel - cannot export to raster format",
)
# Case 2 - Site Model has difference size to shakemap
site_model2 = SiteModel.from_bbox(
[-3.0, -2.0, 3.0, 2.0], 0.5, 0.5, vs30=800.0, z1pt0=10.0
)
# Re-set bounding box properties
shakemap2 = Shakemap(
Earthquake("XYZ", 0.0, 0.0, 10.0, 6.0),
site_model2,
[("B14", GSIM_LIST["BindiEtAl2014Rjb"](), 1.0)],
"ASC",
)
with self.assertRaises(AssertionError) as ae:
shakemap2._transform_to_raster("PGA", results)
self.assertEqual(
str(ae.exception), "Shakemap dimensions do not correspond to site model"
)
def test_export_to_raster_geotiff_file(self):
# Test raster export to geotiff as a round trip
results = np.zeros(self.site_model.shape, dtype=np.dtype([("PGA", np.float64)]))
results["PGA"] = np.copy(self.dummy)
results_as_raster = self.shakemap._transform_to_raster("PGA", results)
geotiff_fname = os.path.join(DATA_PATH, "dummy_raster.tif")
self.shakemap.to_geotiff(results, "PGA", geotiff_fname)
# Load back in the geotiff and verify that the data content is the same
data = rasterio.open(geotiff_fname)
self.assertEqual(data.height, results_as_raster.shape[0])
self.assertEqual(data.width, results_as_raster.shape[1])
np.testing.assert_array_almost_equal(data.read(1), results_as_raster)
os.remove(geotiff_fname)
def test_export_to_raster_geotiff_bytes(self):
# Test raster export to geotiff as an in memory file - compares this against the
# case that the raster is exported directly to file and verifies that the relevant
# content is the same
results = np.zeros(self.site_model.shape, dtype=np.dtype([("PGA", np.float64)]))
results["PGA"] = np.copy(self.dummy)
data_bytes = self.shakemap.to_geotiff(results, "PGA")
byte_shakemap_file = os.path.join(DATA_PATH, "dummy_shakemap_from_bytes.tif")
with open(byte_shakemap_file, "wb") as f:
f.write(data_bytes)
# Verify that this produces the same shakemap as direct export
shakemap_file = os.path.join(DATA_PATH, "dummy_shakemap.tif")
self.shakemap.to_geotiff(results, "PGA", shakemap_file)
# Compare the two files
raster_1 = rasterio.open(byte_shakemap_file)
raster_2 = rasterio.open(shakemap_file)
# Compare bounding boxes
for pos in ["left", "right", "top", "bottom"]:
self.assertAlmostEqual(getattr(raster_1.bounds, pos), getattr(raster_2.bounds, pos))
# Compare dimensions