From 8f1ab5518a1940f2d0f6c82ec2d1602cefd669f4 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Wed, 7 May 2025 12:35:14 +0200 Subject: [PATCH 1/9] ISEPOS-2373 Add saving geographical coordinates to SVG files as metadata --- src/seismic_hazard_forecasting.py | 41 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index b1efe81..c711f4f 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -468,6 +468,16 @@ verbose: {verbose}") else: iml_grid_prep = iml_grid_raw + # TODO Remove. Saving coordinates to + # north, south = lat.max(), lat.min() + # east, west = lon.max(), lon.min() + # bbox_json = {"south": float(south), "west": float(west), + # "north": float(north), "east": float(east)} + # with open("overlay_bounds.json", "w", encoding="utf‑8") as fh: + # json.dump(bbox_json, fh, indent=2) + # + # logger.info(f"Saved bbox to overlay_bounds.json → {bbox_json}") + for j in range(0, len(products)): vmin = min(x for x in iml_grid_prep[j] if x is not np.nan) vmax = max(x for x in iml_grid_prep[j] if x is not np.nan) @@ -498,7 +508,38 @@ verbose: {verbose}") plt.savefig("overlay_" + str(j) + ".svg", bbox_inches="tight", pad_inches=0, transparent=True) plt.close(fig) + + import xml.etree.ElementTree as ET # <-- TODO place this with other imports + + # ----------------------------------------- + # Inject bounding box (BBOX) metadata into the SVG + # ----------------------------------------- + north, south = lat.max(), lat.min() + east, west = lon.max(), lon.min() + bbox_dict = { + "south": float(south), + "west": float(west), + "north": float(north), + "east": float(east) + } + tree = ET.parse(svg_path) + + root = tree.getroot() + + # Remove any existing tags + for old_meta in root.findall("{http://www.w3.org/2000/svg}metadata"): + root.remove(old_meta) + + # Add new element with the bounding box JSON + meta_elem = ET.SubElement(root, "metadata") + meta_elem.text = json.dumps(bbox_dict) + + # ----------------------------------------- + # END Inject bounding box (BBOX) metadata into the SVG + # ----------------------------------------- + # Make the color bar + cmap_name = 'viridis' width = 50 height = 500 From 22fc9f7c073c254974830f8ddbb38adfe922fc54 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Wed, 7 May 2025 12:52:06 +0200 Subject: [PATCH 2/9] ISEPOS-2373 Fix: name 'svg_path' is not defined --- src/seismic_hazard_forecasting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index c711f4f..dd8cb7e 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -505,7 +505,8 @@ verbose: {verbose}") # Save the figure fig.canvas.draw() - plt.savefig("overlay_" + str(j) + ".svg", bbox_inches="tight", pad_inches=0, transparent=True) + svg_path = f"overlay_{j}.svg" + plt.savefig(svg_path, bbox_inches="tight", pad_inches=0, transparent=True) plt.close(fig) From c20f7c06a7a63b2b98478bf07dfe3c9b34d24259 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Wed, 7 May 2025 12:54:42 +0200 Subject: [PATCH 3/9] ISEPOS-2373 Fix: json not defined --- src/seismic_hazard_forecasting.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index dd8cb7e..55a2807 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -511,6 +511,7 @@ verbose: {verbose}") import xml.etree.ElementTree as ET # <-- TODO place this with other imports + import json # ----------------------------------------- # Inject bounding box (BBOX) metadata into the SVG From cce0cd258dbc1d78bedd00a66f530c7d1ef8e2d7 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Wed, 7 May 2025 12:59:21 +0200 Subject: [PATCH 4/9] ISEPOS-2373 Also store the bounding box as a data attribute for quick access --- src/seismic_hazard_forecasting.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 55a2807..0d00306 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -536,6 +536,11 @@ verbose: {verbose}") meta_elem = ET.SubElement(root, "metadata") meta_elem.text = json.dumps(bbox_dict) + # (Optional) Also store the bounding box as a data attribute for quick access + root.set("data-bbox", json.dumps(bbox_dict)) + tree.write(svg_path, encoding="utf-8", xml_declaration=True) + logger.info(f"Embedded bbox into {svg_path} → {bbox_dict}") + # ----------------------------------------- # END Inject bounding box (BBOX) metadata into the SVG # ----------------------------------------- From 17cbcc8e79c8ff803385ae72c144f22222750fcf Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Fri, 9 May 2025 10:26:27 +0200 Subject: [PATCH 5/9] ISEPOS-2373 Simplify code --- src/seismic_hazard_forecasting.py | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 0d00306..10ace1e 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -468,16 +468,6 @@ verbose: {verbose}") else: iml_grid_prep = iml_grid_raw - # TODO Remove. Saving coordinates to - # north, south = lat.max(), lat.min() - # east, west = lon.max(), lon.min() - # bbox_json = {"south": float(south), "west": float(west), - # "north": float(north), "east": float(east)} - # with open("overlay_bounds.json", "w", encoding="utf‑8") as fh: - # json.dump(bbox_json, fh, indent=2) - # - # logger.info(f"Saved bbox to overlay_bounds.json → {bbox_json}") - for j in range(0, len(products)): vmin = min(x for x in iml_grid_prep[j] if x is not np.nan) vmax = max(x for x in iml_grid_prep[j] if x is not np.nan) @@ -528,16 +518,8 @@ verbose: {verbose}") root = tree.getroot() - # Remove any existing tags - for old_meta in root.findall("{http://www.w3.org/2000/svg}metadata"): - root.remove(old_meta) - - # Add new element with the bounding box JSON - meta_elem = ET.SubElement(root, "metadata") - meta_elem.text = json.dumps(bbox_dict) - # (Optional) Also store the bounding box as a data attribute for quick access - root.set("data-bbox", json.dumps(bbox_dict)) + root.set("data-map-bounds", json.dumps(bbox_dict)) tree.write(svg_path, encoding="utf-8", xml_declaration=True) logger.info(f"Embedded bbox into {svg_path} → {bbox_dict}") @@ -546,7 +528,6 @@ verbose: {verbose}") # ----------------------------------------- # Make the color bar - cmap_name = 'viridis' width = 50 height = 500 From 9c58664770d8b05cdbf1dc1ee83a57df2c29f8f3 Mon Sep 17 00:00:00 2001 From: Mieszko Makuch Date: Fri, 9 May 2025 10:30:00 +0200 Subject: [PATCH 6/9] ISEPOS-2373 Simplify code --- src/seismic_hazard_forecasting.py | 40 +++++++++++---------------------------- 1 file changed, 11 insertions(+), 29 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 10ace1e..5e415cd 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -70,6 +70,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max import matplotlib.pyplot as plt from matplotlib.ticker import MultipleLocator from matplotlib.contour import ContourSet + import xml.etree.ElementTree as ET + import json logger = getDefaultLogger('igfash') @@ -495,37 +497,17 @@ verbose: {verbose}") # Save the figure fig.canvas.draw() - svg_path = f"overlay_{j}.svg" - plt.savefig(svg_path, bbox_inches="tight", pad_inches=0, transparent=True) + overlay_filename = f"overlay_{j}.svg" + plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True) plt.close(fig) - - import xml.etree.ElementTree as ET # <-- TODO place this with other imports - import json - - # ----------------------------------------- - # Inject bounding box (BBOX) metadata into the SVG - # ----------------------------------------- - north, south = lat.max(), lat.min() - east, west = lon.max(), lon.min() - bbox_dict = { - "south": float(south), - "west": float(west), - "north": float(north), - "east": float(east) - } - tree = ET.parse(svg_path) - - root = tree.getroot() - - # (Optional) Also store the bounding box as a data attribute for quick access - root.set("data-map-bounds", json.dumps(bbox_dict)) - tree.write(svg_path, encoding="utf-8", xml_declaration=True) - logger.info(f"Embedded bbox into {svg_path} → {bbox_dict}") - - # ----------------------------------------- - # END Inject bounding box (BBOX) metadata into the SVG - # ----------------------------------------- + # Embed geographic bounding box into the SVG + map_bounds = dict(zip(("south", "west", "north", "east"), + map(float, (lat.min(), lon.min(), lat.max(), lon.max())))) + tree = ET.parse(overlay_filename) + tree.getroot().set("data-map-bounds", json.dumps(map_bounds)) + tree.write(overlay_filename, encoding="utf-8", xml_declaration=True) + logger.info(f"Saved geographic bounds to SVG metadata (data-map-bounds): {overlay_filename} → {map_bounds}") # Make the color bar cmap_name = 'viridis' From 7a39d5a07e828c90a502fba8f2aaa434f0aff7b9 Mon Sep 17 00:00:00 2001 From: ftong Date: Fri, 23 May 2025 16:43:15 +0200 Subject: [PATCH 7/9] Update src/seismic_hazard_forecasting.py --- src/seismic_hazard_forecasting.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 5e415cd..0c48502 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -241,6 +241,9 @@ verbose: {verbose}") grid_y_max = int(ceil(y_max / grid_dim) * grid_dim) grid_y_min = int(floor(y_min / grid_dim) * grid_dim) + grid_lat_max, grid_lon_max = utm.to_latlon(grid_x_max, grid_y_max, utm_zone_number, utm_zone_letter) + grid_lat_min, grid_lon_min = utm.to_latlon(grid_x_min, grid_y_min, utm_zone_number, utm_zone_letter) + # rectangular grid nx = int((grid_x_max - grid_x_min) / grid_dim) + 1 ny = int((grid_y_max - grid_y_min) / grid_dim) + 1 @@ -503,7 +506,7 @@ verbose: {verbose}") # Embed geographic bounding box into the SVG map_bounds = dict(zip(("south", "west", "north", "east"), - map(float, (lat.min(), lon.min(), lat.max(), lon.max())))) + map(float, (grid_lat_min, grid_lon_min, grid_lat_max, grid_lon_max)))) tree = ET.parse(overlay_filename) tree.getroot().set("data-map-bounds", json.dumps(map_bounds)) tree.write(overlay_filename, encoding="utf-8", xml_declaration=True) From f1472bf2508ebe2bb3c120c9ef136ae35bf3c1c8 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 4 Jun 2025 17:15:19 +0200 Subject: [PATCH 8/9] add non-parallel processing and set as default; fix colourbar --- src/seismic_hazard_forecasting.py | 50 ++++++++++++++++++++++++++------------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 0c48502..78d7211 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -439,22 +439,40 @@ verbose: {verbose}") PGA = np.zeros(shape=(nx * ny)) - # use dask parallel computing - start = timer() - pbar = ProgressBar() - pbar.register() - # iter = range(0,len(distances)) - iter = indices - iml_grid_raw = [] # raw ground motion grids - for imt in products: - logger.info(f"Estimating {imt}") - imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, - forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, - log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, - rx_label=i) for i in iter] - iml = dask.compute(*imls) - iml_grid_raw.append(list(iml)) + start = timer() + + use_pp = False + + if use_pp: # use dask parallel computing + pbar = ProgressBar() + pbar.register() + # iter = range(0,len(distances)) + iter = indices + iml_grid_raw = [] # raw ground motion grids + for imt in products: + logger.info(f"Estimating {imt}") + + imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, + forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, + log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, + rx_label=i, rtol=0.1, use_cython=False) for i in iter] + iml = dask.compute(*imls) + iml_grid_raw.append(list(iml)) + + else: + iml_grid_raw = [] + iter = indices + for imt in products: + iml = [] + for i in iter: + iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, + lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, + IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=False) + iml.append(iml_i) + logger.info(f"Estimated {imt} at rx {i} is {iml_i}") + iml_grid_raw.append(iml) + end = timer() logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds") @@ -530,7 +548,7 @@ verbose: {verbose}") tick_positions = np.linspace(vmin, vmax, num_ticks) ax.set_yticks(tick_positions) ax.set_yticklabels([f"{tick:.2f}" for tick in tick_positions]) # format tick labels - ax.set_title(imt, pad=15) + ax.set_title(products[j], pad=15) fig.subplots_adjust(left=0.25, right=0.75, bottom=0.05, top=0.95) # Adjust Layout fig.savefig("colorbar_" + str(j) + ".svg", bbox_inches='tight') plt.close(fig) From 846078352b5d9aae804ae6f386b05294ac551257 Mon Sep 17 00:00:00 2001 From: ftong Date: Wed, 4 Jun 2025 17:49:54 +0200 Subject: [PATCH 9/9] revert f1472bf2508ebe2bb3c120c9ef136ae35bf3c1c8 revert add non-parallel processing and set as default; fix colourbar --- src/seismic_hazard_forecasting.py | 48 ++++++++++++--------------------------- 1 file changed, 15 insertions(+), 33 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 78d7211..0c48502 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -439,40 +439,22 @@ verbose: {verbose}") PGA = np.zeros(shape=(nx * ny)) - + # use dask parallel computing start = timer() + pbar = ProgressBar() + pbar.register() + # iter = range(0,len(distances)) + iter = indices + iml_grid_raw = [] # raw ground motion grids + for imt in products: + logger.info(f"Estimating {imt}") - use_pp = False - - if use_pp: # use dask parallel computing - pbar = ProgressBar() - pbar.register() - # iter = range(0,len(distances)) - iter = indices - iml_grid_raw = [] # raw ground motion grids - for imt in products: - logger.info(f"Estimating {imt}") - - imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, - forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, - log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, - rx_label=i, rtol=0.1, use_cython=False) for i in iter] - iml = dask.compute(*imls) - iml_grid_raw.append(list(iml)) - - else: - iml_grid_raw = [] - iter = indices - for imt in products: - iml = [] - for i in iter: - iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, - lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, - IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=False) - iml.append(iml_i) - logger.info(f"Estimated {imt} at rx {i} is {iml_i}") - iml_grid_raw.append(iml) - + imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, + forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, + log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, + rx_label=i) for i in iter] + iml = dask.compute(*imls) + iml_grid_raw.append(list(iml)) end = timer() logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds") @@ -548,7 +530,7 @@ verbose: {verbose}") tick_positions = np.linspace(vmin, vmax, num_ticks) ax.set_yticks(tick_positions) ax.set_yticklabels([f"{tick:.2f}" for tick in tick_positions]) # format tick labels - ax.set_title(products[j], pad=15) + ax.set_title(imt, pad=15) fig.subplots_adjust(left=0.25, right=0.75, bottom=0.05, top=0.95) # Adjust Layout fig.savefig("colorbar_" + str(j) + ".svg", bbox_inches='tight') plt.close(fig)