Newer
Older
Graeme Weatherill
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
"""
Defines complete end-to-end shakemap calculation workflows
"""
import os
from typing import Optional, Dict, List, Tuple, Union
from shakyground2.earthquake import Earthquake
from shakyground2.site_model import SiteModel
from shakyground2.regionalization import DEFAULT_GLOBAL_REGIONALIZATION, RegionalizationSet
from shakyground2.shakemap import Shakemap
DEFAULT_CONFIGURATION = {
# Bounding Box
"spcx": 1.0 / 120.0, # 30 arc-seconds
"spcy": 1.0 / 120.0, # 30 arc-seconds
"max_distance_bbox": None,
# Site Properties
"default_vs30": 600.0,
"vs30measured": True,
"z1pt0": None,
"z2pt5": None,
"xvf": 0.0,
"backarc": False,
"region": 0,
"geology": "UNCLASSIFIED",
# Shakemap
"num_rupture_samples": 100,
"rdim": 0.0,
"synth_dist_weights": [0.25, 0.25, 0.25, 0.25],
"synthetic_rupture_max_site_distance": 200.0,
"synthetic_rupture_site_spacing": 0.05,
}
class ShakemapWorkflowResult(object):
"""
Object to hold output information from shakemap calculations, including
configuration metadata
"""
__slots__ = [
"shakemaps",
"contours",
"num_sites",
"bbox_properties",
"ground_motion_models",
"tectonic_region",
"earthquake",
]
def __getitem__(self, key):
return getattr(self, key)
def shakemaps_from_quakeml(
event_id: str,
imts: List = ["PGA", "SA(0.3)", "SA(1.0)"],
config: Optional[Dict] = None,
export_folder: Optional[str] = None,
cache_file: Optional[str] = None,
contour_levels_mean: Union[int, List] = 10,
contour_levels_stddev: Union[int, List] = 10,
regionalization: RegionalizationSet = DEFAULT_GLOBAL_REGIONALIZATION,
) -> Tuple[Dict, Dict]:
"""
Implements the complete shakemap workflow for use with the GEOFON FDSN Web Service
Args:
event_id: GEOFON event ID or path to QuakeML file for event
imts: List of intensity measure types e.g., PGA, SA(0.1), SA(0.2), ... etc.
config: Dictionary of configuration properties (will use the GEOFON default
configuration if not supplied
export_folder: Path to export the geotiff and contour data (if storing locally)
cache_file: File to cache the shakemap data (optional)
contour_levels_mean: Number of levels for contouring the mean shakemaps, or list of
pre-defined levels
contour_levels_stddev: Number of levels for contouring the standard deviation
shakemaps, or list of pre-defined levels
Returns:
Tuple of dictionaries containing the mean and standard deviation shakemaps
as instances of BytesIO objects, organised by intensity measure type
"""
if not config:
config = DEFAULT_CONFIGURATION
if export_folder:
if os.path.exists(export_folder):
raise IOError("Designated export folder %s already exists!" % export_folder)
os.mkdir(export_folder)
# Create the event from the GEOFON event ID (or the path to the QuakeML file)
results = ShakemapWorkflowResult()
results.earthquake = Earthquake.from_quakeml(event_id)
# Split the configuration into those parts relating to the bounding box,
# the site model and the shakemap
bbox_config = {}
for key in ["spcx", "spcy", "max_distance_bbox"]:
bbox_config[key] = config.get(key, DEFAULT_CONFIGURATION[key])
site_config = {}
for key in [
"default_vs30",
"vs30measured",
"z1pt0",
"z2pt5",
"xvf",
"backarc",
"region",
"geology",
]:
site_config[key] = config.get(key, DEFAULT_CONFIGURATION[key])
# Build the site model
bbox = results.earthquake.get_maximum_distance_bbox(bbox_config["max_distance_bbox"])
vs30 = site_config.pop("default_vs30")
site_model = SiteModel.from_bbox(
bbox, bbox_config["spcx"], bbox_config["spcy"], vs30, **site_config
)
results.num_sites = len(site_model)
results.bbox_properties = site_model.bbox_properties
# Get the ground motion models
results.tectonic_region, results.ground_motion_models = regionalization(results.earthquake)
shakemap_config = {}
for key in [
"num_rupture_samples",
"rdim",
"synth_dist_weights",
"synthetic_rupture_max_site_distance",
"synthetic_rupture_site_spacing",
]:
shakemap_config[key] = config.get(key, DEFAULT_CONFIGURATION[key])
# Run the shakemap
shakemap = Shakemap(
results.earthquake,
site_model,
results.ground_motion_models,
results.tectonic_region,
cache_file=cache_file,
**shakemap_config
)
mean_shakemaps, stddev_shakemaps, _ = shakemap.get_shakemap(imts)
# Export to file (if an export directory is given) or to a dictionary of byte arrays
results.shakemaps = {"mean": {}, "stddevs": {}}
results.contours = {"mean": {}, "stddevs": {}}
for imt in imts:
results.shakemaps["mean"][imt] = shakemap.to_geotiff(mean_shakemaps, imt)
results.shakemaps["stddevs"][imt] = shakemap.to_geotiff(
stddev_shakemaps, imt, is_stddev=True
)
results.contours["mean"][imt] = shakemap.get_contours(
imt, mean_shakemaps, contour_levels_mean
)
results.contours["stddevs"][imt] = shakemap.get_contours(
imt, stddev_shakemaps, contour_levels_stddev, is_stddev=True
)
if export_folder:
filestem = "{:s}_{:s}".format(results.earthquake.id, imt)
# Export the bytes to raster files
fname_raster_mean = os.path.join(export_folder, filestem + "_mean.tif")
with open(fname_raster_mean, "wb") as f:
f.write(results.shakemaps["mean"][imt])
fname_raster_stddev = os.path.join(export_folder, filestem + "_stddev.tif")
with open(fname_raster_stddev, "wb") as f:
f.write(results.shakemaps["stddevs"][imt])
# Export the contour dataframes to geojson
fname_contour_mean = os.path.join(export_folder, filestem + "_contour_mean.geojson")
results.contours["mean"][imt].to_file(fname_contour_mean, driver="GeoJSON")
if results.contours["stddevs"][imt].shape[0]:
# If all the sites have the same standard deviation then skip this as the
# contours will yield an empty dataframe
fname_contour_stddev = os.path.join(
export_folder, filestem + "_contour_stddev.geojson"
)
results.contours["stddevs"][imt].to_file(fname_contour_stddev, driver="GeoJSON")
return results