From 82381c1e6d86eb2b03337ddd997e53e39f8a2bb0 Mon Sep 17 00:00:00 2001 From: unknown Date: Fri, 10 Jan 2025 13:05:49 +0100 Subject: [PATCH 1/2] add --- openquake/man/checks/plotting.py | 268 ++++++++++++++++++ openquake/man/tests/checks_test.py | 39 +++ .../checks/expected/fault_sections.geojson | 6 + .../tests/data/checks/in/ssm/source_model.xml | 47 +++ 4 files changed, 360 insertions(+) create mode 100644 openquake/man/tests/checks_test.py create mode 100644 openquake/man/tests/data/checks/expected/fault_sections.geojson create mode 100644 openquake/man/tests/data/checks/in/ssm/source_model.xml diff --git a/openquake/man/checks/plotting.py b/openquake/man/checks/plotting.py index 2e0aa8057..121a63295 100644 --- a/openquake/man/checks/plotting.py +++ b/openquake/man/checks/plotting.py @@ -1,6 +1,21 @@ +import os import numpy as np +import pandas as pd +import geopandas as gpd import matplotlib.pyplot as plt +from glob import glob +from shapely import Polygon as ShapelyPolygon, LineString +from openquake.commonlib import readinput +from openquake.hazardlib.nrml import GeometryModel +from openquake.hazardlib.geo.mesh import RectangularMesh +from openquake.hazardlib.geo.surface.simple_fault import SimpleFaultSurface +from openquake.hazardlib.source.multi_fault import MultiFaultSource +from openquake.hazardlib.source.complex_fault import ComplexFaultSource +from openquake.hazardlib.source.simple_fault import SimpleFaultSource +from openquake.hazardlib.source.kite_fault import KiteFaultSource +from openquake.hazardlib.source.characteristic import CharacteristicFaultSource +from openquake.hazardlib.source.complex_fault import ComplexFaultSurface from openquake.hazardlib.source.point import PointSource @@ -58,3 +73,256 @@ def plot_polygon(poly): def plot_end(): plt.show() + + +def get_ssm_files(model_dir): + """ + For a given source model get the xml ssm files + """ + # Get seismic source model XMLs + base_path = os.path.join(model_dir, 'in', 'ssm') + try: + assert os.path.exists(base_path) + except: + raise ValueError("Admitted model does not abide to required directory " + "hierarchy: model_dir/in/ssm/") + files = [] + + # Loop through each subfolder and get all XML files + for depth in range(5): # Search up to 5 levels deep + fipath = os.path.join(base_path, *['**']*depth) + files.extend(glob(fipath + "/*.xml", recursive=True)) + + # Remove duplicates from recursive search + files = pd.Series(files).drop_duplicates().values + + return files + + +def get_sources(model_dir): + """ + Load the sources in the given model and return them as a list. + """ + # Get source XMLs + files = get_ssm_files(model_dir) + + # Fault classess + faults = [MultiFaultSource, + SimpleFaultSource, + ComplexFaultSource, + CharacteristicFaultSource, + KiteFaultSource] + + # Read the XMLs for all srcs in the given model + ssm = readinput.read_source_models(files, 'tmp.hdf5', investigation_time=1, + rupture_mesh_spacing=5, + area_source_discretization=5, + width_of_mfd_bin=0.1) + # Get a list of the source objects + srcs = [] + geom_model = [] + for idx, ssm_obj in enumerate(ssm): + if isinstance(ssm_obj, GeometryModel): + geom_model.append(ssm_obj) + ssm = pd.Series(ssm).drop(idx).values + else: + for idx_sg, sg in enumerate(ssm_obj.src_groups): + for src in ssm_obj[idx_sg]: + if type(src) in faults: # Only retain the fault sources + srcs.append(src) + assert len(geom_model) < 2 # Should be one or none + if len(geom_model) == 1: + geom_model = geom_model[0] + else: + geom_model = None + + return srcs, geom_model + + +def get_boundary_2d(smsh): + """ + Returns a polygon for the given surface. + """ + coo = [] + + # Upper boundary + trace + idx = np.where(np.isfinite(smsh.lons[0, :]))[0] + tmp = [(smsh.lons[0, i], smsh.lats[0, i]) for i in idx] + trace = LineString(tmp) + coo.extend(tmp) + + # Right boundary + idx = np.where(np.isfinite(smsh.lons[:, -1]))[0] + tmp = [(smsh.lons[i, -1], smsh.lats[i, -1]) for i in idx] + coo.extend(tmp) + + # Lower boundary + idx = np.where(np.isfinite(smsh.lons[-1, :]))[0] + tmp = [(smsh.lons[-1, i], smsh.lats[-1, i]) for i in np.flip(idx)] + coo.extend(tmp) + + # Left boundary + idx = idx = np.where(np.isfinite(smsh.lons[:, 0]))[0] + tmp = [(smsh.lons[i, 0], smsh.lats[i, 0]) for i in np.flip(idx)] + coo.extend(tmp) + + return trace, ShapelyPolygon(coo) + + +def get_complex_mesh(src): + """ + Get the mesh of a ComplexFaultSource + """ + # Get the surface + sfc = ComplexFaultSurface.from_fault_data(src.edges, mesh_spacing=5) + + # Get mesh + mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths) + + return mesh + + +def get_characteristic_mesh(src): + """ + Get the mesh of a CharacteristicFaultSource + """ + # Get the surface + sfc = src.surface.mesh + + # Get mesh + mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths) + + return mesh + + +def get_simple_mesh(src): + """ + Get the mesh of a SimpleFaultSource + """ + # Get the surface + sfc = SimpleFaultSurface.from_fault_data( + src.fault_trace, src.upper_seismogenic_depth, + src.lower_seismogenic_depth, src.dip, mesh_spacing=5) + + # Get mesh + mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths) + + return mesh + + +def get_kite_mesh(src): + """ + Get the mesh of a KiteFaultSource + """ + # Get the sfc (idl set to false as not automated yet in engine) + sfc = src.surface.from_profiles(src.profiles, src.profiles_sampling, + src.rupture_mesh_spacing, idl=False, + align=False) + + # Get mesh + mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths) + + return mesh + + +def get_geoms(srcs, geom_model): + """ + Extract the geometry of each fault source and write to a geoJSON + """ + traces = [] + polys = [] + suids = [] + + # First non-MFS sources + for i, src in enumerate(srcs): + if isinstance(src, MultiFaultSource): # Skip MFS as use the geom model + continue # object to get these geometries + if isinstance(src, ComplexFaultSource): + surf = get_complex_mesh(src) + elif isinstance(src, CharacteristicFaultSource): + surf = get_characteristic_mesh(src) + elif isinstance(src, SimpleFaultSource): + surf = get_simple_mesh(src) + elif isinstance(src, KiteFaultSource): + surf = get_kite_mesh(src) + else: + raise ValueError(f"Unknown source typology admitted: ({type(src)})") + trace, poly = get_boundary_2d(surf) + traces.append(trace) + polys.append(poly) + suids.append(i) + + # Get geometries of the MFS sources too + if geom_model: + for i, key in enumerate(geom_model.sections): + surf = geom_model.sections[key] + trace, poly = get_boundary_2d(surf.mesh) + traces.append(trace) + polys.append(poly) + suids.append(i) + + daf = pd.DataFrame({'suid': suids, 'geometry': polys}) + gdaf_polys = gpd.GeoDataFrame(daf, geometry='geometry') + daf = pd.DataFrame({'suid': suids, 'geometry': traces}) + gdaf_traces = gpd.GeoDataFrame(daf, geometry='geometry') + + return gdaf_polys, gdaf_traces + + +def export_faults(gdaf_polys, gdaf_traces, model_dir): + """ + Export the GeoDataFrames of the fault traces and fault polygons + """ + out_traces = os.path.join(model_dir, 'fault_traces.geojson') + out_polys = os.path.join(model_dir, 'fault_sections.geojson') + gdaf_traces.to_file(out_traces) + gdaf_polys.to_file(out_polys) + print(f"Fault traces exported to {out_traces}") + print(f"Fault sections exported to {out_polys}") + + return out_polys, out_traces + + +def plot_faults(gdaf_polys, gdaf_traces, region, model_dir): + """ + Plot the faults within the geoJSONs + """ + import pygmt + out = os.path.join(model_dir, "faults_plot.png") + fig = pygmt.Figure() + fig.basemap(region=region, projection="M15c", frame=True) + fig.coast(land="grey", water="skyblue") + fig.plot(gdaf_polys, pen="0.5p,red") + fig.plot(gdaf_traces, pen="1p,green") + fig.savefig(out) + print(f"Faults plotted and saved to {out}") + fig.show() + + +def get_fault_geojsons(model_dir, plotting=False, plotting_region=None): + """ + Write the fault sections and fault traces within the given hazard model + model to geojsons. + + :param model_dir: directory containing the required hazard model + + :param plotting: Boolean which if True creates a plot using the geoJSONs + of the faults in the given hazard model + + :param plotting_region: list of [min lon, max lon, min lat, max lat] used + to define axis limits of the plotted geoJSONs + """ + # Get the sources in the given model + srcs, geom_model = get_sources(model_dir) + + # Now get the geometries + gdaf_polys, gdaf_traces = get_geoms(srcs, geom_model) + + # Export into geoJSONs + out_polys, out_traces = export_faults(gdaf_polys, gdaf_traces, model_dir) + + # Plot them if specified + if plotting: + plot_faults(gdaf_polys, gdaf_traces, plotting_region, model_dir) + + return gdaf_polys, gdaf_traces, out_polys, out_traces \ No newline at end of file diff --git a/openquake/man/tests/checks_test.py b/openquake/man/tests/checks_test.py new file mode 100644 index 000000000..08c052fe6 --- /dev/null +++ b/openquake/man/tests/checks_test.py @@ -0,0 +1,39 @@ +import os +import numpy as np +import geopandas as gpd +import unittest + +from openquake.man.checks.plotting import get_fault_geojsons + +BASE_DATA_PATH = os.path.dirname(__file__) + +exp_polys_path = os.path.join( + BASE_DATA_PATH, 'data', 'checks', 'expected', 'fault_sections.geojson') + + +class TestReadModel(unittest.TestCase): + + def test_get_fault_geoJSONs(self): + """ + Check execution of the functions which convert the fault sources + within a hazard model into geoJSONs. + + Uses a source XML with a single SimpleFaultSource. + """ + fname = os.path.join(BASE_DATA_PATH, 'data', 'checks') + results = get_fault_geojsons(fname) + + # Check the geoJSON contents + obs_polys = results[0] + exp_polys = gpd.read_file(exp_polys_path) + obs_p_coo = obs_polys.geometry[0].exterior.coords._coords + exp_p_coo = exp_polys.geometry[0].exterior.coords._coords + np.testing.assert_allclose(obs_p_coo, exp_p_coo) + + # Check the geojsons exist + assert os.path.exists(results[2]) # polys geojson + assert os.path.exists(results[3]) # traces geojson + + # Then remove them + os.remove(results[2]) + os.remove(results[3]) diff --git a/openquake/man/tests/data/checks/expected/fault_sections.geojson b/openquake/man/tests/data/checks/expected/fault_sections.geojson new file mode 100644 index 000000000..d1e1819d0 --- /dev/null +++ b/openquake/man/tests/data/checks/expected/fault_sections.geojson @@ -0,0 +1,6 @@ +{ +"type": "FeatureCollection", +"features": [ +{ "type": "Feature", "properties": { "suid": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 1.021395468520495, -0.274887161684799 ], [ 1.061615435442728, -0.254777787522597 ], [ 1.101835125551899, -0.234668373923735 ], [ 1.142054573532541, -0.21455891841192 ], [ 1.182273814068017, -0.194449418510723 ], [ 1.222492881840751, -0.174339871743573 ], [ 1.262711811532374, -0.154230275633758 ], [ 1.302930637823961, -0.134120627704419 ], [ 1.343149395396204, -0.114010925478548 ], [ 1.383368118929616, -0.093901166478988 ], [ 1.424116359046308, -0.074887187950001 ], [ 1.469082479230706, -0.074887358892875 ], [ 1.514048659583617, -0.074887468242194 ], [ 1.559014927801107, -0.074887515997653 ], [ 1.603981311579775, -0.074887502158835 ], [ 1.648947838616836, -0.074887426725208 ], [ 1.693914536610364, -0.074887289696126 ], [ 1.738881433259303, -0.074887091070829 ], [ 1.738881433259303, -0.074887091070829 ], [ 1.749579076381868, -0.112330675622099 ], [ 1.760276746916998, -0.149774256257488 ], [ 1.770974454002376, -0.187217831671684 ], [ 1.781672206775785, -0.224661400559353 ], [ 1.781672206775785, -0.224661400559353 ], [ 1.73670531012761, -0.224661599184538 ], [ 1.691738612134616, -0.224661736213544 ], [ 1.646772085097836, -0.224661811647129 ], [ 1.601805701319218, -0.22466182548594 ], [ 1.55683943310155, -0.224661777730507 ], [ 1.511873252748232, -0.224661668381249 ], [ 1.466907132563147, -0.224661497438471 ], [ 1.426158969044339, -0.243675465362055 ], [ 1.385940331651136, -0.263785213144987 ], [ 1.345721665490477, -0.283894904154221 ], [ 1.305502935881849, -0.304004540866913 ], [ 1.265284108144636, -0.324114125760065 ], [ 1.225065147597765, -0.344223661310532 ], [ 1.184846019559678, -0.364333149995023 ], [ 1.144626689347953, -0.384442594290103 ], [ 1.104407122279304, -0.404551996672199 ], [ 1.064187283669246, -0.424661359617599 ], [ 1.064187283669246, -0.424661359617599 ], [ 1.05348919727219, -0.387217829077592 ], [ 1.042791205374247, -0.349774285038472 ], [ 1.032093298836565, -0.312330728805726 ], [ 1.021395468520495, -0.274887161684799 ] ] ] } } +] +} diff --git a/openquake/man/tests/data/checks/in/ssm/source_model.xml b/openquake/man/tests/data/checks/in/ssm/source_model.xml new file mode 100644 index 000000000..106b5d144 --- /dev/null +++ b/openquake/man/tests/data/checks/in/ssm/source_model.xml @@ -0,0 +1,47 @@ + + + + + + + + + 1.0000000E+00 -2.0000000E-01 1.4000000E+00 0.0000000E+00 1.7000000E+00 0.0000000E+00 + + + + 3.0000000E+01 + + + 5.0000000E+00 + + + 1.5000000E+01 + + + + WC1994 + + + 2.0000000E+00 + + + + 9.0000000E+01 + + + + + From 687dc67a8d80e981d413068c40646fc95e4b410f Mon Sep 17 00:00:00 2001 From: unknown Date: Fri, 10 Jan 2025 14:41:49 +0100 Subject: [PATCH 2/2] upd --- openquake/man/tests/checks_test.py | 39 ------------------------------ 1 file changed, 39 deletions(-) delete mode 100644 openquake/man/tests/checks_test.py diff --git a/openquake/man/tests/checks_test.py b/openquake/man/tests/checks_test.py deleted file mode 100644 index 08c052fe6..000000000 --- a/openquake/man/tests/checks_test.py +++ /dev/null @@ -1,39 +0,0 @@ -import os -import numpy as np -import geopandas as gpd -import unittest - -from openquake.man.checks.plotting import get_fault_geojsons - -BASE_DATA_PATH = os.path.dirname(__file__) - -exp_polys_path = os.path.join( - BASE_DATA_PATH, 'data', 'checks', 'expected', 'fault_sections.geojson') - - -class TestReadModel(unittest.TestCase): - - def test_get_fault_geoJSONs(self): - """ - Check execution of the functions which convert the fault sources - within a hazard model into geoJSONs. - - Uses a source XML with a single SimpleFaultSource. - """ - fname = os.path.join(BASE_DATA_PATH, 'data', 'checks') - results = get_fault_geojsons(fname) - - # Check the geoJSON contents - obs_polys = results[0] - exp_polys = gpd.read_file(exp_polys_path) - obs_p_coo = obs_polys.geometry[0].exterior.coords._coords - exp_p_coo = exp_polys.geometry[0].exterior.coords._coords - np.testing.assert_allclose(obs_p_coo, exp_p_coo) - - # Check the geojsons exist - assert os.path.exists(results[2]) # polys geojson - assert os.path.exists(results[3]) # traces geojson - - # Then remove them - os.remove(results[2]) - os.remove(results[3])