Commit a5b9f7e7 authored by g-weatherill's avatar g-weatherill
Browse files

Fixs maximum distance bug, sets explicit maximum distance

parent b8dfdaed
Pipeline #40894 failed with stage
in 18 minutes and 54 seconds
......@@ -11,7 +11,7 @@ import matplotlib.pyplot as plt
from typing import Dict, Optional, Tuple, List, Union
from rasterio.warp import reproject
from openquake.hazardlib import const, imt
from openquake.hazardlib.contexts import ContextMaker
from openquake.hazardlib.contexts import ContextMaker, IntegrationDistance
from shakyground2 import valid
from shakyground2.earthquake import Earthquake
from shakyground2.site_model import SiteModel
......@@ -19,6 +19,9 @@ from shakyground2.synthetic_rupture_generator import FiniteRuptureSampler
from shakyground2.utils import plt_contours_to_gpd_dataframe
MAX_DIST = "1.0E5"
class Shakemap(object):
"""
The core class for calculating shakemaps. The general workflow contained within
......@@ -121,10 +124,15 @@ class Shakemap(object):
Construct the rupture, site and distances contexts from the earthquake, site and
ground motion models
"""
# param is a dummy object to mock an OpenQuake configuation
param = {
"imtls": {"PGA": []},
"maximum_distance": IntegrationDistance.new(MAX_DIST),
}
cmaker = ContextMaker(
self.tectonic_region,
[gmm[1] for gmm in self.ground_motion_model],
oq={"imtls": {"PGA": []}}, # param is a dummy object to mock an OQ config
oq=param,
)
# Always calculate rupture distance regardless of whether it is needed by the GMM
if "rrup" not in cmaker.REQUIRES_DISTANCES: # pylint: disable=no-member
......@@ -218,7 +226,6 @@ class Shakemap(object):
del fle[self.earthquake.id]
fle_eq = fle.create_group(self.earthquake.id)
self._cache_contexts(fle_eq)
# Pre-allocate the aggregated shakemaps with zeros
aggregated_means = np.zeros(self.site_model.shape, dtype=shakemap_dtypes)
aggregated_stddevs = np.zeros(self.site_model.shape, dtype=shakemap_dtypes)
......
......@@ -271,7 +271,7 @@ class SiteModel(object):
# 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:
elif SITE_PROPERTIES[key] == bool:
site_array[key] = cls._get_boolean_array(nsites, value)
else:
site_array[key] = value * np.ones(nsites, dtype=SITE_PROPERTIES[key])
......
......@@ -82,9 +82,7 @@ def plt_contours_to_gpd_dataframe(
contour_geometries.append(geometries[0])
contour_levels.append(level)
dframe = GeoDataFrame(
{name: contour_levels, "geometry": GeoSeries(contour_geometries)}
)
dframe = GeoDataFrame({name: contour_levels, "geometry": GeoSeries(contour_geometries)})
if crs:
dframe.crs = crs
return dframe
......@@ -255,7 +255,7 @@ class SiteModelTestCase(unittest.TestCase):
# 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,):
if site_array1[name].dtype not in (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])
......
......@@ -26,9 +26,7 @@ class PltContoursToGpdDataframeTestCase(unittest.TestCase):
ax = fig.add_subplot(111)
contours = ax.contour(x, y.T[0], h, levels)
# First important point: Don't raise any exception here.
dframe = plt_contours_to_gpd_dataframe(
contours, name="h", crs={"init": "epsg:4326"}
)
dframe = plt_contours_to_gpd_dataframe(contours, name="h", crs={"init": "epsg:4326"})
plt.close()
plt.ion()
# And also make sure the geometry is valid
......
Supports Markdown
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