diff --git a/.github/workflows/engine_nightly_webui.yml b/.github/workflows/engine_nightly_webui.yml index 3e9f6a9c860d..969d64a7821d 100644 --- a/.github/workflows/engine_nightly_webui.yml +++ b/.github/workflows/engine_nightly_webui.yml @@ -51,9 +51,16 @@ jobs: set -x date wget https://downloads.openquake.org/test_data/exposure.hdf5 + wget https://downloads.openquake.org/test_data/worldcities.csv + wget https://downloads.openquake.org/test_data/countries_info.csv + wget https://downloads.openquake.org/test_data/World_Adm1_updated.gpkg + wget https://downloads.openquake.org/test_data/fonts.zip + unzip fonts.zip date mv exposure.hdf5 openquake/qa_tests_data/mosaic/ source ~/openquake/bin/activate + PY_VER=$(echo "${{ matrix.python-version }}" | tr -d '.') + pip install -r requirements-impact-py${PY_VER}-linux64.txt oq engine --upgrade-db touch ~/webui-access.log pip install https://wheelhouse.openquake.org/v3/py/pytest_django-4.9.0-py3-none-any.whl diff --git a/.github/workflows/engine_pr_test.yml b/.github/workflows/engine_pr_test.yml index 7e37fc231fe3..8e995c4f7b5d 100644 --- a/.github/workflows/engine_pr_test.yml +++ b/.github/workflows/engine_pr_test.yml @@ -37,7 +37,7 @@ jobs: run: | set -e source ~/openquake/bin/activate - pip install https://wheelhouse.openquake.org/v3/py/rtgmpy-1.0.1-py3-none-any.whl + pip install -r requirements-aelo-py-linux64.txt oq engine --upgrade-db pytest --doctest-modules -x --disable-warnings --color=yes --durations=10 openquake/calculators -k "not classical" @@ -70,7 +70,8 @@ jobs: run: | set -e source ~/openquake/bin/activate - pip install pyshp pytest flake8 ruff https://wheelhouse.openquake.org/v3/py/rtgmpy-1.0.1-py3-none-any.whl + pip install -r requirements-aelo-py-linux64.txt + pip install pyshp pytest flake8 ruff oq engine --upgrade-db ruff check openquake curl https://downloads.openquake.org/test_data/exposure.hdf5 -o openquake/qa_tests_data/mosaic/exposure.hdf5 @@ -110,7 +111,14 @@ jobs: source ~/openquake/bin/activate oq engine --upgrade-db pip install https://wheelhouse.openquake.org/v3/py/pytest_django-4.9.0-py3-none-any.whl + pip install -r requirements-impact-py311-linux64.txt + pip install -r requirements-aelo-py-linux64.txt curl https://downloads.openquake.org/test_data/exposure.hdf5 -o openquake/qa_tests_data/mosaic/exposure.hdf5; ln -s openquake/qa_tests_data/mosaic/exposure.hdf5 + wget https://downloads.openquake.org/test_data/worldcities.csv + wget https://downloads.openquake.org/test_data/countries_info.csv + wget https://downloads.openquake.org/test_data/World_Adm1_updated.gpkg + wget https://downloads.openquake.org/test_data/fonts.zip + unzip fonts.zip OQ_APPLICATION_MODE=PUBLIC pytest -vx openquake/server/tests/test_public_mode.py OQ_APPLICATION_MODE=RESTRICTED pytest -vx openquake/server/tests/test_restricted_mode.py OQ_APPLICATION_MODE=READ_ONLY pytest -vx openquake/server/tests/test_read_only_mode.py diff --git a/.github/workflows/engine_weekly_test.yml b/.github/workflows/engine_weekly_test.yml index 20527f3d0577..ac10d8917f9a 100644 --- a/.github/workflows/engine_weekly_test.yml +++ b/.github/workflows/engine_weekly_test.yml @@ -417,6 +417,9 @@ jobs: else python install.py devel --version master fi + PY_VER=`echo py${{ matrix.python-version }} | tr -d .` + echo $PY_VER + pip3 install -r requirements-impact-$PY_VER-linux64.txt - name: Actualize 'impact' templates for email notifications run: | for file in ./openquake/server/templates/registration/*.impact.tmpl; do @@ -429,6 +432,11 @@ jobs: set -x date wget https://downloads.openquake.org/test_data/exposure.hdf5 + wget https://downloads.openquake.org/test_data/worldcities.csv + wget https://downloads.openquake.org/test_data/countries_info.csv + wget https://downloads.openquake.org/test_data/World_Adm1_updated.gpkg + wget https://downloads.openquake.org/test_data/fonts.zip + unzip fonts.zip date mv exposure.hdf5 openquake/qa_tests_data/mosaic/ source ~/openquake/bin/activate diff --git a/bin/extract_impact_reports.py b/bin/extract_impact_reports.py new file mode 100644 index 000000000000..a7f66db9714a --- /dev/null +++ b/bin/extract_impact_reports.py @@ -0,0 +1,35 @@ +import sys +import subprocess +from pathlib import Path +from openquake.commonlib import datastore + + +def main(calc_id: int): + dstore = datastore.read(calc_id) + if "impact" not in dstore: + raise RuntimeError("No 'impact' group found in datastore") + impact_group = dstore["impact"] + iso3_codes = list(impact_group.keys()) + if not iso3_codes: + raise RuntimeError("No ISO3 reports found under 'impact/'") + print(f"Found ISO3 reports: {iso3_codes}") + for iso3 in iso3_codes: + country_group = impact_group[iso3] + if "report_pdf" not in country_group: + print(f"No report_pdf found for {iso3}") + continue + pdf_bytes = country_group["report_pdf"][()] + fname = Path(f"impact_report_{calc_id}_{iso3}.pdf") + with open(fname, "wb") as f: + f.write(pdf_bytes) + print(f"Wrote {fname}") + subprocess.Popen(["xdg-open", str(fname)]) + dstore.close() + + +if __name__ == "__main__": + if len(sys.argv) != 2: + print("Usage: extract_impact_reports.py ") + sys.exit(1) + + main(int(sys.argv[1])) diff --git a/openquake/calculators/country_impact_report_builder.py b/openquake/calculators/country_impact_report_builder.py new file mode 100644 index 000000000000..fc983fcb2468 --- /dev/null +++ b/openquake/calculators/country_impact_report_builder.py @@ -0,0 +1,770 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# vim: tabstop=4 shiftwidth=4 softtabstop=4 +# +# Copyright (C) 2026 GEM Foundation +# +# OpenQuake is free software: you can redistribute it and/or modify it +# under the terms of the GNU Affero General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# OpenQuake is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with OpenQuake. If not, see . + +import os +import logging +from io import BytesIO +from pathlib import Path +from PIL import Image as PILImage +from openquake import baselib +from openquake.baselib import config +from openquake.calculators.country_impact_report_utils import ( + EventContext, ReportOptions, LOSS_METADATA, _read_countries_info, + _read_world_cities, build_classifiers, load_admin_boundaries, + points_to_gdf, aggregate_losses, save_most_affected_regions) +from openquake.calculators.postproc.plots import plot_variable, MapDataElements +from openquake.hazardlib.calc.filters import upper_maxdist + +cd = Path(__file__).parent + +COUNTRY_PROFILES_BASE_URL = "https://github.com/gem/risk-profiles/tree/master" + + +class CountryImpactReportBuilder: + """ + Builds and stores a single-country impact report. + """ + # Layout constants + MARGIN = 20 + DISCLAIMER_H = 40 + HEADER_H = 80 + NOTES_H = 80 + SAFETY_BUFFER = 20 + LOGO_W = 100 + + def __init__( + self, iso3, adm_level, event: EventContext, options: ReportOptions, + losses_df, summary_data, dstore, time_of_calc, oqparam): + try: + import reportlab + from reportlab import platypus + except ImportError as exc: + raise RuntimeError( + "In order to create an impact report," + " 'reportlab' should be installed" + ) from exc + try: + import fitz # PyMuPDF + except ImportError as exc: + raise RuntimeError( + "In order to save an Impact report as PNG," + " 'PyMuPDF' should be installed" + ) from exc + + self.fitz = fitz + self.reportlab = reportlab + self.SimpleDocTemplate = platypus.SimpleDocTemplate + self.Paragraph = platypus.Paragraph + self.Table = platypus.Table + self.TableStyle = platypus.TableStyle + self.Image = platypus.Image + self.ListFlowable = platypus.ListFlowable + self.ListItem = platypus.ListItem + self.Spacer = platypus.Spacer + self.getSampleStyleSheet = reportlab.lib.styles.getSampleStyleSheet + self.ParagraphStyle = reportlab.lib.styles.ParagraphStyle + self.colors = reportlab.lib.colors + self.A4 = reportlab.lib.pagesizes.A4 + + self.iso3 = iso3 + self.adm_level = adm_level + self.losses_df = losses_df + self.summary_data = summary_data + self.dstore = dstore + self.time_of_calc = time_of_calc + self.maximum_distance = upper_maxdist(oqparam.maximum_distance) + + # Unpacking EventContext + self.event_name = event.name + self.event_date = event.date + self.shakemap_version = event.shakemap_version + self.hypocenter = event.hypocenter + + # Unpacking ReportOptions + self.disclaimer_txt = options.disclaimer_txt + self.basemap_path = options.basemap_path + self.threshold_deg = options.threshold_deg + self.no_uncertainty = options.no_uncertainty + self.loss_metric = options.loss_metric + + self.styles = self.getSampleStyleSheet() + + self.x_limits = None + self.y_limits = None + self.cities = {} + + self._load_country_info() + self.notes = self._get_notes(oqparam) + self._compute_layout() + + self._register_unicode_font() + self.styles["Normal"].fontName = "NotoSans" + self.styles["Italic"].fontName = "NotoSans-Italic" + self.styles["Heading1"].fontName = "NotoSans-Bold" + + def _register_unicode_font(self): + from reportlab.pdfbase import pdfmetrics + from reportlab.pdfbase.ttfonts import TTFont + from reportlab.pdfbase.pdfmetrics import registerFontFamily + + try: + fonts_dir = config.directory.fonts_dir + except AttributeError: + # checking if the directory is present in the oq-engine directory + if not os.path.exists( + fonts_dir := cd.parent.parent / 'fonts'): + raise AttributeError( + 'config.directory.fonts_dir is missing') + fonts_dir = Path(fonts_dir) + + # family_name -> font file prefix + font_families = { + "NotoSans": "NotoSans", # Latin, Cyrillic, Greek + "NotoSans-SC": "NotoSansSC", # Simplified Chinese + "NotoSans-TC": "NotoSansTC", # Traditional Chinese + "NotoSans-JP": "NotoSansJP", # Japanese + "NotoSans-KR": "NotoSansKR", # Korean + "NotoSans-AR": "NotoSansArabic", # Arabic + "NotoSans-Deva": "NotoSansDevanagari", # Hindi, Nepali, etc. + "NotoSans-Beng": "NotoSansBengali", # Bengali + "NotoSans-Thai": "NotoSansThai", # Thai + } + for family, name in font_families.items(): + regular = fonts_dir / f"{name}-Regular.ttf" + bold = fonts_dir / f"{name}-Bold.ttf" + if not regular.exists(): + logging.warning(f"Font not found: {regular}, skipping") + continue + bold_path = str(bold) if bold.exists() else str(regular) + pdfmetrics.registerFont(TTFont(family, str(regular))) + pdfmetrics.registerFont(TTFont(f"{family}-Bold", bold_path)) + pdfmetrics.registerFont(TTFont(f"{family}-Italic", str(regular))) + registerFontFamily( + family, + normal=family, + bold=f"{family}-Bold", + italic=f"{family}-Italic", + boldItalic=f"{family}-Bold", + ) + + def _select_font(self, text): + """Pick the right font family based on Unicode block detection.""" + text = str(text) # handle non-string input gracefully + for ch in text: + cp = ord(ch) + if 0x0600 <= cp <= 0x06FF: + return "NotoSans-AR" + if 0x0900 <= cp <= 0x097F: + return "NotoSans-Deva" + if 0x0980 <= cp <= 0x09FF: + return "NotoSans-Beng" + if 0x0E00 <= cp <= 0x0E7F: + return "NotoSans-Thai" + if 0xAC00 <= cp <= 0xD7AF: + return "NotoSans-KR" + if 0x3040 <= cp <= 0x309F: + return "NotoSans-JP" # Hiragana + if 0x30A0 <= cp <= 0x30FF: + return "NotoSans-JP" # Katakana + if 0x4E00 <= cp <= 0x9FFF: + return "NotoSans-SC" + if 0xF900 <= cp <= 0xFAFF: + return "NotoSans-TC" + return "NotoSans" + + def _one_line_paragraph( + self, text, base_style, max_width, min_font_size=8, step=0.5): + """ + Try to keep paragraph on one line by reducing font size if needed. + """ + font_size = base_style.fontSize + + font_name = self._select_font(text) + if font_name != base_style.fontName: + base_style = self.ParagraphStyle( + name="tmp_font", + parent=base_style, + fontName=font_name, + ) + + while font_size >= min_font_size: + style = self.ParagraphStyle( + name="tmp", + parent=base_style, + fontSize=font_size, + leading=font_size * 1.2, + ) + p = self.Paragraph(text, style) + w, h = p.wrap(max_width, 1000) + + # one line ≈ leading height + if h <= style.leading * 1.05: + return p + + font_size -= step + + # fallback: smallest font + style.fontSize = min_font_size + style.leading = min_font_size * 1.2 + return self.Paragraph(text, style) + + def _scaled_image(self, path, max_w, max_h): + # scale image preserving aspect ratio + if not path.exists(): + return self.Paragraph(f"Missing image: {path.name}", + self.getSampleStyleSheet()["Normal"]) + img = PILImage.open(path) + w, h = img.size + scale = min(max_w / w, max_h / h) + return self.Image(str(path), width=w*scale, height=h*scale) + + def _scaled_image_from_bytes(self, image_data, max_w, max_h): + # If it is an HDF5 dataset, read it + try: + image_data = image_data[()] + except Exception: + pass + # scale image preserving aspect ratio + img = PILImage.open(BytesIO(image_data)) + w, h = img.size + scale = min(max_w / w, max_h / h) + return self.Image(BytesIO(image_data), width=w*scale, height=h*scale) + + def _compute_viewport_from_boundaries( + self, aggregated_gdf, padding_deg=0.5): + """ + Derive the map viewport from the bounding box of admin boundaries + of all regions with at least one non-zero loss + """ + loss_labels = [meta["label"] for meta in LOSS_METADATA.values()] + mask = aggregated_gdf[loss_labels].gt(0).any(axis=1) + affected = aggregated_gdf[mask] + bounds = affected.geometry.total_bounds # (minx, miny, maxx, maxy) + return ( + [bounds[0] - padding_deg, bounds[2] + padding_deg], + [bounds[1] - padding_deg, bounds[3] + padding_deg], + ) + + def _load_country_info(self): + try: + countries_info_file = config.directory.countries_info_file + except AttributeError: + # checking if the file is present in the oq-engine directory + if not os.path.exists( + countries_info_file := cd.parent.parent / + 'countries_info.csv'): + raise AttributeError( + 'config.directory.countries_info_file is missing') + + path_str = str(Path(countries_info_file).resolve()) + df = _read_countries_info(path_str) # cached + row = df.loc[df["ISO3"] == self.iso3].iloc[0] + self.country_name = row["ENGLISH_COUNTRY"] + self.country_region = row["GEM_REGION"] + + def _get_notes(self, oqparam): + notes_data = { + "user_note": oqparam.notes if oqparam.notes else None, + "profile_link": None, + "metadata": [] + } + country_profile_link = (f'{COUNTRY_PROFILES_BASE_URL}/' + f'{self.country_region}/' + f'{self.country_name}') + notes_data["profile_link"] = ( + f"Seismic Risk Profile for the Country: " + f"" + f"{country_profile_link}" + ) + rupdic = oqparam.rupture_dict + meta = notes_data["metadata"] + meta.append(f'USGS identifier: {rupdic["usgs_id"]}') + meta.append(f'Longitude: {rupdic["lon"]}') + meta.append(f'Latitude: {rupdic["lat"]}') + meta.append(f'Depth: {rupdic["dep"]}') + meta.append(f'Magnitude: {rupdic["mag"]}') + meta.append(f'Rake: {rupdic["rake"]}') + meta.append(f'Dip: {rupdic["dip"]}') + meta.append(f'Strike: {rupdic["strike"]}') + if rupdic['approach'] != 'use_shakemap_from_usgs': + meta.append(f'Mosaic model: {oqparam.mosaic_model}') + meta.append( + f'Tectonic region type: {oqparam.tectonic_region_type}') + meta.append(f'Number of ground motion fields:' + f' {oqparam.number_of_ground_motion_fields}') + meta.append(f'Truncation level: {oqparam.truncation_level}') + meta.append(f'Considered time of the event: {oqparam.time_event}') + return notes_data + + def _get_cities_in_viewport(self, num_cities=15): + """ + Finds Top num_cities cities within the map viewport belonging + to the current country + """ + try: + # NOTE: using for the report a file structured differently with + # respect to openquake/qa_tests_data/mosaic/worldcities.csv + # We may want to replace the other file with this, changing also + # the expected column names. + world_cities_file = config.directory.world_cities_file + except AttributeError: + # checking if the file is present in the oq-engine directory + if not os.path.exists( + world_cities_file := cd.parent.parent / + 'worldcities.csv'): + raise AttributeError( + 'config.directory.world_cities_file is missing') + path_str = str(Path(world_cities_file).resolve()) + df = _read_world_cities(path_str) # cached + # Pull the pre-calculated limits + min_lon, max_lon = self.x_limits + min_lat, max_lat = self.y_limits + # Spatial query + Country filter + mask = (df['iso3'] == self.iso3) & \ + (df['lng'] >= min_lon) & (df['lng'] <= max_lon) & \ + (df['lat'] >= min_lat) & (df['lat'] <= max_lat) + # Take the biggest ones + top_cities = df[mask].sort_values( + 'population', ascending=False).head(num_cities) + return {row['city_ascii']: [row['lng'], row['lat']] + for _, row in top_cities.iterrows()} + + def _compute_layout(self): + self.page_width = self.A4[0] - (2 * self.MARGIN) + self.page_height = self.A4[1] - (2 * self.MARGIN) + + self.grid_total_h = ( + self.page_height + - self.DISCLAIMER_H + - self.HEADER_H + - self.NOTES_H + - self.SAFETY_BUFFER + ) + + self.row_h = self.grid_total_h / 2 + self.col_w = self.page_width / 2 + + def _generate_country_plots(self): + import matplotlib.pyplot as plt + tags_agg_losses = list(LOSS_METADATA) + admin_boundaries = load_admin_boundaries( + self.country_name, self.iso3, self.adm_level) + points_gdf = points_to_gdf(self.losses_df, crs=admin_boundaries.crs) + aggloss_df = aggregate_losses( + points_gdf, admin_boundaries, tags_agg_losses) + aggloss_df = aggloss_df.rename(columns={k: v["label"] + for k, v in LOSS_METADATA.items()}) + save_most_affected_regions(aggloss_df, self.dstore, self.iso3) + self.x_limits, self.y_limits = self._compute_viewport_from_boundaries( + aggloss_df) + self.cities = self._get_cities_in_viewport() + + classifiers = build_classifiers(aggloss_df, breaks=[1, 10, 100, 1000]) + images = {} + for meta in LOSS_METADATA.values(): + label = meta["label"] + plot_title = f'{self.loss_metric} {meta["title"]}' + elements = MapDataElements( + plot_title=plot_title, + # legend_title=label, # already in plot title + cities=self.cities, + x_limits=self.x_limits, + y_limits=self.y_limits, + basemap_path=self.basemap_path, + epicenter=self.hypocenter + ) + fig, ax = plot_variable( + aggloss_df, admin_boundaries, label, + classifiers[label], meta["colors"], + elements=elements + ) + buf = BytesIO() + fig.savefig(buf, format="png", dpi=300, bbox_inches="tight") + plt.close(fig) + buf.seek(0) + images[label] = buf.getvalue() + return images + + def _build_disclaimer(self): + tbl = self.Table( + [[self.Paragraph(f"DISCLAIMER: {self.disclaimer_txt}", + self.styles["Normal"])]], + colWidths=[self.page_width], + rowHeights=[self.DISCLAIMER_H], + ) + + tbl.setStyle(self.TableStyle([ + ("BACKGROUND", (0, 0), (-1, -1), self.colors.lightcoral), + ("BOX", (0, 0), (-1, -1), 1, self.colors.red), + ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), + ("LEFTPADDING", (0, 0), (-1, -1), 10), + ])) + + return tbl + + def _build_header(self): + event_style = self.ParagraphStyle( + "EventTitle", + parent=self.styles["Normal"], + fontName="NotoSans-Bold", + fontSize=12, + leading=14, + ) + + meta_style = self.ParagraphStyle( + "HeaderMeta", + parent=self.styles["Normal"], + fontSize=9, + leading=11 + ) + + # subtracting also the padding + title_width = self.page_width - self.LOGO_W - 12 + + # Line 1: Bold event name + event_text = f"{self.event_name}" + event_paragraph = self._one_line_paragraph( + event_text, + event_style, + max_width=title_width, + ) + + oq_basedir = Path(baselib.__path__[0].rsplit("/", 2)[0]) + logo_path = ( + oq_basedir + / "doc" + / "_static" + / "OQ-Logo-Standard-RGB-72DPI-01.png" # FIXME: is this logo ok? + ) + + logo_img = self._scaled_image( + logo_path, + self.LOGO_W, + self.HEADER_H - 10, + ) + + # Build individual paragraph blocks for lines 2, 3, and 4 + sm_version_txt = None + if self.shakemap_version is not None: + sm_version_txt = f'ShakeMap version: {self.shakemap_version}' + date_txt = f"Time of the event: {self.event_date}" + calc_txt = f"Time of the calculation: {self.time_of_calc}" + header_text = [event_paragraph] + if sm_version_txt: + header_text.append( + self._one_line_paragraph( + sm_version_txt, meta_style, max_width=title_width)) + header_text.extend([ + self._one_line_paragraph( + date_txt, meta_style, max_width=title_width), + self._one_line_paragraph( + calc_txt, meta_style, max_width=title_width), + ]) + + tbl = self.Table( + [[header_text, logo_img]], + colWidths=[self.page_width - self.LOGO_W, self.LOGO_W], + rowHeights=[self.HEADER_H], + ) + + tbl.setStyle(self.TableStyle([ + ("VALIGN", (0, 0), (-1, -1), "TOP"), + ("ALIGN", (1, 0), (1, 0), "RIGHT"), # logo to the right + ("TOPPADDING", (0, 0), (-1, -1), 10), + ("LEFTPADDING", (0, 0), (-1, -1), 0), + ("RIGHTPADDING", (0, 0), (-1, -1), 0), + ])) + + return tbl + + def _grid_styles(self): + body_left_style = self.ParagraphStyle( + "GridBodyTextLeft", + parent=self.styles["Normal"], + fontSize=9, + leading=11, + alignment=0, # Left-aligned for text labels + ) + body_right_style = self.ParagraphStyle( + "GridBodyTextRight", + parent=self.styles["Normal"], + fontSize=9, + leading=11, + alignment=2, # Right-aligned for numeric metrics + ) + title_style = self.ParagraphStyle( + "GridSectionTitle", + parent=self.styles["Normal"], + fontName="NotoSans-Bold", + fontSize=11, + leading=14, + ) + return body_left_style, body_right_style, title_style + + def _build_summary_table(self, body_left_style, body_right_style): + col_header = ("Estimated losses" if self.no_uncertainty + else "Range of losses (5% - 95%)") + table_data = [[ + self.Paragraph("Impact metric", body_left_style), + self.Paragraph("Exposed value", body_right_style), + self.Paragraph(f"{col_header}", body_right_style) + ]] + for meta in LOSS_METADATA.values(): + # NOTE: in order to make it easier to understand and communicate, + # we use 'residents' for both 'Fatalities' and 'Displaced' + if meta["label"] in ["Fatalities", "Displaced"]: + exposed_key = LOSS_METADATA["residents"]["label"] + "_exposed" + else: + exposed_key = meta["label"] + "_exposed" + table_data.append([ + self.Paragraph(meta["label"], body_left_style), + self.Paragraph(self.summary_data[exposed_key], + body_right_style), + self.Paragraph(self.summary_data[meta["label"]], + body_right_style) + ]) + summary_table = self.Table( + table_data, + colWidths=[self.col_w * 0.32, + self.col_w * 0.32, + self.col_w * 0.32], + hAlign="LEFT", + ) + style_cmds = [ + ("GRID", (0, 0), (-1, -1), 0.5, self.colors.grey), + ("BACKGROUND", (0, 0), (-1, 0), self.colors.whitesmoke), + ("PADDING", (0, 0), (-1, -1), 4), + ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), + ] + summary_table.setStyle(self.TableStyle(style_cmds)) + return summary_table + + def _build_left_bundle(self, summary_table, body_left_style, title_style): + most_affected = self.dstore[ + f"impact/{self.iso3}/most_affected_regions" + ] + left_bundle = [ + self.Paragraph( + f"Summary of impact for {self.country_name}:", + title_style), + self.Spacer(1, 4), + summary_table, + self.Spacer(1, 6), + ] + + exposed_value_txt = ( + f'The exposed value refers to the assets and population located' + f' within a {self.maximum_distance}km radius of the epicentre.') + left_bundle.append( + self.Paragraph(exposed_value_txt, body_left_style)) + if self.no_uncertainty: + left_bundle.extend([ + self.Spacer(1, 4), + self.Paragraph("No uncertainty was included", + body_left_style)]) + left_bundle.append(self.Spacer(1, 18)) + left_bundle.extend([ + self.Paragraph("Regions with highest number of fatalities:", + title_style), + self.ListFlowable( + [self.ListItem(self.Paragraph(region_name, self.ParagraphStyle( + "region", + parent=body_left_style, + fontName=self._select_font(region_name), + ))) for region_name in most_affected], + bulletType="bullet", + leftIndent=15, + ), + ]) + return left_bundle + + # NOTE: passing images explicitly to avoid implicit ordering dependency + def _build_grid(self, images): + body_left_style, body_right_style, title_style = self._grid_styles() + summary_table = self._build_summary_table(body_left_style, + body_right_style) + left_bundle = self._build_left_bundle( + summary_table, body_left_style, title_style) + + img_top_right = self._scaled_image_from_bytes( + images[LOSS_METADATA['number']['label']], + self.col_w - 10, + self.row_h - 10, + ) + img_bot_left = self._scaled_image_from_bytes( + images[LOSS_METADATA['occupants']['label']], + self.col_w - 10, + self.row_h - 10, + ) + img_bot_right = self._scaled_image_from_bytes( + images[LOSS_METADATA['residents']['label']], + self.col_w - 10, + self.row_h - 10, + ) + tbl = self.Table( + [ + [left_bundle, img_top_right], + [img_bot_left, img_bot_right], + ], + colWidths=[self.col_w, self.col_w], + rowHeights=[self.row_h, self.row_h], + ) + tbl.setStyle(self.TableStyle([ + ("VALIGN", (0, 0), (-1, -1), "TOP"), + ("ALIGN", (1, 0), (1, 1), "CENTER"), + ("ALIGN", (0, 1), (0, 1), "CENTER"), + ("TOPPADDING", (0, 0), (-1, -1), 5), + ("BOTTOMPADDING", (0, 0), (-1, -1), 5), + ])) + return tbl + + def _build_notes(self): + """ + Builds a bordered notes box with dedicated, dynamic rows for full-width + user notes and web links, anchoring a 3-column metadata grid below. + """ + story = [] + if not self.notes or not isinstance(self.notes, dict): + return story + grid_data = [] + + styles_to_apply = [ + ('BOX', (0, 0), (-1, -1), 1, self.reportlab.lib.colors.black), + ('VALIGN', (0, 0), (-1, -1), 'TOP'), + ('LEFTPADDING', (0, 0), (-1, -1), 5), + ('RIGHTPADDING', (0, 0), (-1, -1), 5), + ('TOPPADDING', (0, 0), (-1, -1), 3), + ('BOTTOMPADDING', (0, 0), (-1, -1), 3), + ] + + # Create a compact font style variant for the notes box + notes_style = self.ParagraphStyle( + 'NotesStyle', + parent=self.styles['Normal'], + fontSize=8, + leading=7.5 + ) + # Header and Optional User Custom Notes Row + user_note = self.notes.get("user_note") + header_text = ( + f"Notes: {user_note}" if user_note else "Notes:") + grid_data.append([self.Paragraph(header_text, notes_style), "", ""]) + styles_to_apply.append(('SPAN', (0, 0), (2, 0))) + # Keep spacing tight below title + styles_to_apply.append(('BOTTOMPADDING', (0, 0), (2, 0), 0)) + + # Optional Standalone Profile Link Row + profile_link = self.notes.get("profile_link") + if profile_link: + grid_data.append( + [self.Paragraph(profile_link, notes_style), "", ""]) + link_row_idx = len(grid_data) - 1 + styles_to_apply.append( + ('SPAN', (0, link_row_idx), (2, link_row_idx))) + styles_to_apply.append( + ('BOTTOMPADDING', (0, link_row_idx), (2, link_row_idx), 2)) + + # Spread the remaining system metadata across 3 columns + notes_items = self.notes.get("metadata", []) + for i in range(0, len(notes_items), 3): + row = [self.Paragraph(item, notes_style) + for item in notes_items[i:i+3]] + while len(row) < 3: + row.append(self.Paragraph("", notes_style)) + grid_data.append(row) + + # Add safety cushion at the bottom of the last metadata row + styles_to_apply.append(('BOTTOMPADDING', (0, -1), (2, -1), 6)) + + t = self.Table(grid_data, colWidths=[180, 180, 180]) + t.setStyle(self.TableStyle(styles_to_apply)) + story.append(t) + return story + + def build(self): + logging.info(f'Making impact report for {self.iso3}...') + images = self._generate_country_plots() + + buffer = BytesIO() + + doc = self.SimpleDocTemplate( + buffer, + pagesize=self.A4, + leftMargin=self.MARGIN, + rightMargin=self.MARGIN, + topMargin=self.MARGIN, + bottomMargin=self.MARGIN, + # A country report is deliberately a single-page document. Do + # not let ReportLab silently split the master table onto another + # page if its contents grow beyond the allocated layout. + allowSplitting=0, + ) + + master_layout = self.Table( + [ + [self._build_disclaimer()], + [self._build_header()], + [self._build_grid(images)], + [self._build_notes()], + ], + colWidths=[self.page_width], + rowHeights=[ + self.DISCLAIMER_H, + self.HEADER_H, + self.grid_total_h, + self.NOTES_H, + ], + ) + + master_layout.setStyle(self.TableStyle([ + ("LEFTPADDING", (0, 0), (-1, -1), 0), + ("RIGHTPADDING", (0, 0), (-1, -1), 0), + ("TOPPADDING", (0, 0), (-1, -1), 0), + ("BOTTOMPADDING", (0, 0), (-1, -1), 0), + ("VALIGN", (0, 0), (-1, -1), "TOP"), + ])) + + doc.build([master_layout]) + + buffer.seek(0) + pdf_bytes = buffer.getvalue() + pdf_doc = self.fitz.open(stream=pdf_bytes, filetype="pdf") + page_count = pdf_doc.page_count + if page_count != 1: + pdf_doc.close() + raise RuntimeError( + f"Impact report for {self.iso3} has " + f"{page_count} pages; expected exactly one") + + pdf_path = f'impact/{self.iso3}/report_pdf' + self.dstore[pdf_path] = pdf_bytes + logging.info( + f'The impact report in PDF format was saved into the datastore' + f' as {pdf_path}') + + # Generate and save an exact PNG duplicate of the layout + page = pdf_doc.load_page(0) + # Render to a crisp image at 3.0x scaling (~300 DPI equivalent) + pix = page.get_pixmap(matrix=self.fitz.Matrix(3.0, 3.0)) + png_path = f'impact/{self.iso3}/report_png' + self.dstore[png_path] = pix.tobytes("png") + pdf_doc.close() + logging.info( + f'The impact report in PNG format was saved into the datastore' + f' as {png_path}') diff --git a/openquake/calculators/country_impact_report_utils.py b/openquake/calculators/country_impact_report_utils.py new file mode 100644 index 000000000000..f2e1440d07e8 --- /dev/null +++ b/openquake/calculators/country_impact_report_utils.py @@ -0,0 +1,202 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# vim: tabstop=4 shiftwidth=4 softtabstop=4 +# +# Copyright (C) 2026 GEM Foundation +# +# OpenQuake is free software: you can redistribute it and/or modify it +# under the terms of the GNU Affero General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# OpenQuake is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with OpenQuake. If not, see . + +import os +import functools +import logging +import pathlib +from dataclasses import dataclass +from shapely.validation import make_valid, explain_validity +import pandas as pd +import geopandas as gpd +from openquake.baselib import config + + +cd = pathlib.Path(__file__).parent + + +@dataclass +class EventContext: + """Metadata related to the seismic event.""" + name: str + date: str + hypocenter: tuple[float, float] + shakemap_version: str = None + + +@dataclass +class ReportOptions: + """Visual, text, and threshold configurations for the report.""" + disclaimer_txt: str + basemap_path: str + threshold_deg: float + no_uncertainty: bool + loss_metric: str + + +LOSS_METADATA = { + "occupants": { + "label": "Fatalities", + "title": "fatalities", + "colors": [ + '#fff5f0', '#fcbba1', '#fb6a4a', '#cb181d', '#67000d'], + }, + "residents": { + "label": "Displaced", + "title": "displaced population", + "colors": [ + '#f1eef6', '#d7b5d8', '#df65b0', '#dd1c77', '#980043'], + }, + "number": { + "label": "Buildings destroyed", + "title": "buildings destroyed", + "colors": [ + '#ffffff', '#bdbdbd', '#737373', '#424242', '#000000'], + }, +} + + +# maxsize=1 is sufficient when only one admin-level boundary file is loaded +# per process (the common case). Increase to 2 if both adm1 and adm2 files +# are ever used within the same process. +@functools.lru_cache(maxsize=1) +def _read_admin_layer(fname): + gdf = gpd.read_file(fname) + invalid = ~gdf.is_valid + if invalid.any(): + for idx in gdf[invalid].index: + reason = explain_validity(gdf.at[idx, "geometry"]) + logging.warning("Invalid geometry at index %s: %s", idx, reason) + # fix invalid geometries + gdf["geometry"] = gdf["geometry"].apply(make_valid) + return gdf + + +@functools.lru_cache(maxsize=1) +def _read_countries_info(countries_info_path): + """ + Load and cache the countries CSV keyed on the resolved file path. + Subsequent calls with the same path return the in-memory DataFrame + without any disk I/O. + """ + return pd.read_csv(countries_info_path) + + +@functools.lru_cache(maxsize=1) +def _read_world_cities(world_cities_path): + """ + Load and cache the world-cities CSV keyed on the resolved file path. + """ + df = pd.read_csv(world_cities_path) + if 'lng' not in df.columns: + raise ValueError(f'Missing "lng" column in {world_cities_path}') + return df + + +def build_classifiers(df, *, breaks): + try: + import mapclassify + except ImportError as exc: + raise RuntimeError( + "In order to build map classifiers 'mapclassify' should" + " be installed." + ) from exc + return {meta["label"]: mapclassify.UserDefined(df[meta["label"]], + bins=breaks) + for meta in LOSS_METADATA.values()} + + +def load_admin_boundaries( + country_name, iso3, adm_level, crs="EPSG:4326"): + if adm_level == 1: + try: + fname = config.directory.admin1_boundaries_file + except AttributeError: + # checking if the file is present in the oq-engine directory + if not os.path.exists( + fname := cd.parent.parent / + 'World_Adm1_updated.gpkg'): + raise AttributeError( + 'config.directory.admin1_boundaries_file is missing') + elif adm_level == 2: + try: + fname = config.directory.admin2_boundaries_file + except AttributeError as exc: + raise AttributeError( + 'config.directory.admin2_boundaries_file is missing') from exc + else: + raise NotImplementedError(f'Admin level {adm_level} not supported') + if not fname: + raise AttributeError( + f'config.directory.admin{adm_level}_boundaries_file is missing') + # NOTE: be careful not mutating the cached object + # (in case we need to mutate it, we should make a copy + # right after reading) + gdf = _read_admin_layer(fname) # cached + if "shapeID" in gdf.columns: # geoBoundaries + iso3_col = "shapeGroup" + id_col = "shapeID" + name_col = "shapeName" + elif f"ID_{adm_level}" in gdf.columns: + iso3_col = "ID_0" + id_col = f"ID_{adm_level}" + name_col = f"NAME_{adm_level}" + else: + raise RuntimeError( + f"Unsupported admin schema. Columns: {list(gdf.columns)}" + ) + # NOTE: here we make a copy, so we don't alter the cached object + gdf = gdf[gdf[iso3_col] == iso3] + if gdf.empty: + raise ValueError( + f"No boundaries found for country '{country_name}'") + # normalize column names + gdf = gdf.rename(columns={ + iso3_col: "country_iso3", + id_col: "region_id", + name_col: "region_name", + }) + gdf["region_id"] = gdf["region_id"].astype(str) + gdf["region_name"] = gdf["region_name"].astype(str) + gdf["country_iso3"] = gdf["country_iso3"].astype(str) + return gdf.to_crs(crs) + + +def points_to_gdf(df, lon_col="lon", lat_col="lat", crs=None): + gdf = gpd.GeoDataFrame( + df, + geometry=gpd.points_from_xy(df[lon_col], df[lat_col]), + crs=crs) + return gdf + + +def aggregate_losses(points_gdf, admin_gdf, tags_agg): + joined = gpd.sjoin(points_gdf, admin_gdf, how="inner", predicate="within") + group_col = 'region_id' + merge_args = dict(on=group_col) + aggregated = joined.groupby(group_col).agg( + {col: "sum" for col in tags_agg}) + return admin_gdf.merge(aggregated, **merge_args) + + +def save_most_affected_regions(df, dstore, iso3, *, num_regions=5): + fatalities_label = LOSS_METADATA["occupants"]["label"] + regions = df.nlargest( + num_regions, fatalities_label)['region_name'].dropna().tolist() + dstore[f"impact/{iso3}/most_affected_regions"] = regions diff --git a/openquake/calculators/event_based.py b/openquake/calculators/event_based.py index fdcb7a81f691..9cad395e3e58 100644 --- a/openquake/calculators/event_based.py +++ b/openquake/calculators/event_based.py @@ -48,6 +48,7 @@ RuptureProxy, EBRupture, get_ruptures_aw) from openquake.hazardlib.shakemap.parsers import adjust_hypocenter from openquake.commonlib import util, logs, readinput, datastore +from openquake.commonlib.readinput import get_close_mosaic_models from openquake.commonlib.calc import ( gmvs_to_poes, make_hmaps, slice_dt, build_slice_by_event, RuptureImporter, SLICE_BY_EVENT_NSITES, get_proxies, get_model_lts) @@ -672,7 +673,7 @@ def read_gsim_lt(oq): elif oq.rupture_xml: hypo = readinput.get_rupture(oq).hypocenter lon, lat = [hypo.x, hypo.y] - mosaic_models = readinput.get_close_mosaic_models(lon, lat, 5) + mosaic_models = get_close_mosaic_models(lon, lat, 5) # NOTE: using the first mosaic model oq.mosaic_model = mosaic_models[0] if len(mosaic_models) > 1: diff --git a/openquake/calculators/export/__init__.py b/openquake/calculators/export/__init__.py index dff8c0077c0a..495ed0e3f1cc 100644 --- a/openquake/calculators/export/__init__.py +++ b/openquake/calculators/export/__init__.py @@ -99,13 +99,12 @@ 'economic': 'Economic loss (USD)', 'occupants': 'Fatalities', 'area': 'Floor area lost (m²)', - 'number': 'Buildings beyond repair', + 'number': 'Buildings destroyed', 'residents': 'Rendered homeless', 'injured': 'Number of injured people', - 'affectedpop': ('Number of people living in buildings ' - 'with moderate or higher damage'), + 'affectedpop': 'Affected population', 'embodied_carbon': 'Embodied carbon loss (tCO₂e)', - 'value': 'Value', + 'value': 'Exposed value', 'loss_type': 'Loss type', 'lossmea': 'Mean', 'q50': 'Median', @@ -113,6 +112,23 @@ 'q95': '95th percentile', } + +AGGRISK_FIELD_EXPLANATION = { + 'affectedpop': ('Population living in buildings that sustained at least' + ' moderate damage'), + 'area': ('Total floor area (in m²) of buildings that sustained complete' + ' damage'), + 'injured': ('Seriously or critically injured people according to the' + ' Abbreviated Injury Severity (www.aisinjuryscale.org)'), + 'number': ('Number of buildings that sustained a level of damage beyond' + ' repair (extensive damage or higher)'), + 'occupants': 'Number of fatalities', + 'residents': ('Population living in buildings that sustained at least' + ' extensive damage'), + 'economic': ('Sum of structural and contents economic loss (in USD)'), +} + + EXPOSURE_FIELD_DESCRIPTION = { 'number': 'Buildings', 'contents': 'Contents value (USD)', diff --git a/openquake/calculators/export/risk.py b/openquake/calculators/export/risk.py index c69c92ac0009..dbbbcaa5ee4a 100644 --- a/openquake/calculators/export/risk.py +++ b/openquake/calculators/export/risk.py @@ -26,11 +26,11 @@ from openquake.baselib import hdf5, writers, general, node from openquake.baselib.general import decode from openquake.hazardlib import nrml -from openquake.hazardlib.stats import compute_stats2, mean_curve +from openquake.hazardlib.stats import mean_curve from openquake.risklib import scientific from openquake.commonlib import readinput from openquake.calculators.extract import ( - extract, sanitize, avglosses, aggexp_tags) + extract, sanitize, _get_data, aggexp_tags) from openquake.calculators import base, post_risk from openquake.calculators.export import export, loss_curves from openquake.calculators.export.hazard import savez @@ -225,32 +225,6 @@ def export_exposure_by_lse(ekey, dstore): return [dest] -def _get_data(dstore, dskey, loss_types, stats): - name, kind = dskey.split('-') # i.e. ('avg_losses', 'stats') - if kind == 'stats': - try: - weights = dstore['weights'][()] - except KeyError: - # there is single realization, like in classical_risk/case_2 - weights = [1.] - if dskey in set(dstore): # precomputed - rlzs_or_stats = list(stats) - statfuncs = [stats[ros] for ros in stats] - value = avglosses(dstore, loss_types, 'stats') # shape (A, S, L) - elif dstore['oqparam'].collect_rlzs: - rlzs_or_stats = list(stats) - value = avglosses(dstore, loss_types, 'rlzs') - else: # compute on the fly - rlzs_or_stats, statfuncs = zip(*stats.items()) - value = compute_stats2( - avglosses(dstore, loss_types, 'rlzs'), statfuncs, weights) - else: # rlzs - value = avglosses(dstore, loss_types, kind) # shape (A, R, L) - R = value.shape[1] - rlzs_or_stats = ['rlz-%03d' % r for r in range(R)] - return name, value, rlzs_or_stats - - # this is used by event_based_risk, classical_risk and scenario_risk @export.add(('avg_losses-rlzs', 'csv'), ('avg_losses-stats', 'csv')) def export_avg_losses(ekey, dstore): @@ -856,7 +830,8 @@ def export_fragility_xml(dstore): crm = dstore.read_df('crm') ddic = {peril: {} for peril in crm.peril.unique()} for (peril, loss_type), df in crm.groupby(['peril', 'loss_type']): - nodeobj = convert_df_to_fragility(peril, loss_type, oq.limit_states, df) + nodeobj = convert_df_to_fragility( + peril, loss_type, oq.limit_states, df) dest = dstore.export_path('%s_%s_fragility.xml' % (peril, loss_type)) with open(dest, 'wb') as out: nrml.write([nodeobj], out) diff --git a/openquake/calculators/extract.py b/openquake/calculators/extract.py index a4b264e8d51c..68652afd7fb6 100644 --- a/openquake/calculators/extract.py +++ b/openquake/calculators/extract.py @@ -38,7 +38,7 @@ from openquake.hazardlib.contexts import ( ContextMaker, read_cmakers, read_ctx_by_grp) from openquake.hazardlib.calc import disagg, stochastic, filters -from openquake.hazardlib.stats import calc_stats +from openquake.hazardlib.stats import calc_stats, compute_stats2 from openquake.hazardlib.source import rupture from openquake.risklib.scientific import LOSSTYPE, LOSSID from openquake.risklib.asset import tagset @@ -867,15 +867,15 @@ def extract_mmi_tags(dstore, what): def ensure_npy_serializable(df): """ - Cast object-dtype columns of a DataFrame to fixed-width numpy byte - strings (dtype ``|S``) so that a structured array built from it - is serializable with ``allow_pickle=False``. + Cast object-dtype and categorical columns of a DataFrame to fixed-width + numpy byte strings (dtype ``|S``) so that a structured array built + from it is serializable with ``allow_pickle=False``. :param df: a :class:`pandas.DataFrame` - :returns: the same DataFrame with object columns replaced in-place + :returns: the same DataFrame with string columns replaced in-place """ for col in df.columns: - if df[col].dtype != object: + if df[col].dtype != object and df[col].dtype.name != 'category': continue # Explicitly encode strings to UTF-8 bytes to handle non-ASCII # characters (e.g., 'é', 'ñ') before casting to the fixed-width byte @@ -1942,6 +1942,52 @@ def extract_high_sites(dstore, what): return (max_hazard > .2).all(axis=1) # shape N +def _get_data(dstore, dskey, loss_types, stats): + name, kind = dskey.split('-') # i.e. ('avg_losses', 'stats') + if kind == 'stats': + try: + weights = dstore['weights'][()] + except KeyError: + # there is single realization, like in classical_risk/case_2 + weights = [1.] + if dskey in set(dstore): # precomputed + rlzs_or_stats = list(stats) + statfuncs = [stats[ros] for ros in stats] + value = avglosses(dstore, loss_types, 'stats') # shape (A, S, L) + elif dstore['oqparam'].collect_rlzs: + rlzs_or_stats = list(stats) + value = avglosses(dstore, loss_types, 'rlzs') + else: # compute on the fly + rlzs_or_stats, statfuncs = zip(*stats.items()) + value = compute_stats2( + avglosses(dstore, loss_types, 'rlzs'), statfuncs, weights) + else: # rlzs + value = avglosses(dstore, loss_types, kind) # shape (A, R, L) + R = value.shape[1] + rlzs_or_stats = ['rlz-%03d' % r for r in range(R)] + return name, value, rlzs_or_stats + + +@extract.add('avg_losses') +def extract_avg_losses(dstore, what): + """ + Example: + http://127.0.0.1:8800/v1/calc/30/extract/avg_losses?kind=stats + """ + oq = dstore['oqparam'] + qdict = parse(what) + dskey = f"avg_losses-{qdict['kind'][0]}" + name, value, rlzs_or_stats = _get_data( + dstore, dskey, oq.ext_loss_types, oq.hazard_stats()) + assets = util.get_assets(dstore) + dt = [(lt, F32) for lt in oq.ext_loss_types] + for ros, values in zip(rlzs_or_stats, value.transpose(1, 0, 2)): + array = numpy.zeros(len(values), dt) + for lt, ln in enumerate(oq.ext_loss_types): + array[ln] = values[:, lt] + yield ros, util.compose_arrays(assets, array) + + # ##################### extraction from the WebAPI ###################### # class WebAPIError(RuntimeError): diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 132d84cd8b1f..944c85be80a0 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -19,7 +19,11 @@ import io import base64 import numpy +import functools +from pathlib import Path from shapely.geometry import Polygon, box +from dataclasses import dataclass +from typing import Any from openquake.commonlib import readinput, datastore from openquake.hmtk.plotting.patch import PolygonPatch @@ -115,7 +119,8 @@ def redraw(event_ax=None): clipped = geom.intersection(viewport) if clipped.is_empty: continue - parts = clipped.geoms if hasattr(clipped, 'geoms') else [clipped] + parts = (clipped.geoms if hasattr(clipped, 'geoms') + else [clipped]) for part in parts: if isinstance(part, Polygon): path = _polygon_to_path(part) @@ -293,11 +298,12 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), coll = ax.scatter(shakemap_array['lon'], shakemap_array['lat'], c=gmf, cmap='jet', s=marker_size) - fig.colorbar(coll, ax=ax, shrink=0.8) + fig.colorbar(coll, ax=ax, shrink=0.6) if rupture is not None: - add_rupture(ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.9, - surf_facecolor='none', surf_linestyle='-', zorder=4) + add_rupture(ax, rupture, hypo_alpha=0.8, hypo_markersize=8, + surf_alpha=0.9, surf_facecolor='none', + surf_linestyle='-', zorder=4) ax.set_xlabel('Longitude') ax.set_ylabel('Latitude') @@ -472,3 +478,237 @@ def get_assetcol(calc_id): except AttributeError: assetcol = dstore['assetcol'].array return assetcol + + +def _resolve_limits(ax, x_limits, y_limits, epicenter=None, buffer_ratio=0.05): + """ + Resolve axis limits ensuring epicenter visibility with a margin. + :param buffer_ratio: fraction of axis span to use as padding + """ + # Start from user limits or current limits + xmin, xmax = x_limits if x_limits else ax.get_xlim() + ymin, ymax = y_limits if y_limits else ax.get_ylim() + + if epicenter is not None: + lon, lat = epicenter + + xmin = min(xmin, lon) + xmax = max(xmax, lon) + ymin = min(ymin, lat) + ymax = max(ymax, lat) + + # Calculate current directional spans + xspan = max(xmax - xmin, 1e-6) + yspan = max(ymax - ymin, 1e-6) + + # Force the spans to be identical to guarantee a square data viewport + if xspan > yspan: + center_y = (ymin + ymax) / 2 + ymin = center_y - (xspan / 2) + ymax = center_y + (xspan / 2) + yspan = xspan + elif yspan > xspan: + center_x = (xmin + xmax) / 2 + xmin = center_x - (yspan / 2) + xmax = center_x + (yspan / 2) + xspan = yspan + + # Apply buffer + xpad = xspan * buffer_ratio + ypad = yspan * buffer_ratio + + xmin -= xpad + xmax += xpad + ymin -= ypad + ymax += ypad + + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + + +def _prepare_classified_data(df, admin_boundaries, column, classifier, colors): + """ + Validates CRS and color counts, then applies classification mapping. + """ + if len(colors) > classifier.k: + colors = colors[:classifier.k] + elif len(colors) < classifier.k: + raise ValueError( + f"Not enough colors: got {len(colors)}, need {classifier.k}. " + f"Please supply at least as many colors as classifier bins." + ) + if df.crs != admin_boundaries.crs: + raise ValueError("df and admin_boundaries CRS do not match") + df = df.copy() + df["class"] = classifier(df[column]) + return df, colors + + +def _build_legend_labels(classifier, legend_digits): + """ + Generates the range strings for the map classification legend. + """ + bins = numpy.round(classifier.bins, legend_digits) + labels = [f"≤ {bins[0]:.{legend_digits}f}"] + labels += [ + f"{bins[i-1]:.{legend_digits}f} – {bins[i]:.{legend_digits}f}" + for i in range(1, len(bins)) + ] + labels[-1] = f"> {bins[-2]:.{legend_digits}f}" + if len(labels) != classifier.k: + raise RuntimeError("Generated labels do not match number of classes") + return labels + + +@functools.lru_cache(maxsize=1) +def _read_basemap(basemap_path): + import rasterio + with rasterio.open(basemap_path) as src: + return src.read(), src.transform, src.crs + + +def _overlay_basemap(ax, basemap_path, target_crs): + """ + Safely handles rasterio imports and displays the basemap overlay. + Raster data is cached after the first load via _read_basemap(). + """ + if basemap_path is None: + return + try: + import rasterio # noqa + except ImportError as exc: + raise RuntimeError( + "In order to plot raster basemaps, 'rasterio' should be installed" + ) from exc + from rasterio.plot import show + # Normalise to str for the cache key; resolve() collapses any symlinks + # or relative segments so the same file always hits the same cache slot. + path_str = str(Path(basemap_path).resolve()) + data, transform, crs = _read_basemap(path_str) + if crs != target_crs: + raise ValueError("Raster CRS does not match vector CRS") + show(data, transform=transform, ax=ax, alpha=0.8) + + +def _overlay_cities(ax, cities, city_font_size): + """ + Plots city markers and dynamically adjusts text placements to prevent + collisions. + """ + if not cities: + return + import matplotlib.patheffects as path_effects + city_scatters = [] + texts = [] + for city, (x, y) in cities.items(): + sc = ax.scatter(x, y, color='black', marker='o', s=8, zorder=6) + city_scatters.append(sc) + t = ax.text(x, y, city, fontsize=city_font_size, color="black", + zorder=7) + t.set_path_effects([ + path_effects.Stroke(linewidth=1.5, foreground="white"), + path_effects.Normal() + ]) + texts.append(t) + try: + from adjustText import adjust_text + if texts: + legend = ax.get_legend() + adjust_text( + texts, ax=ax, add_objects=city_scatters + [legend], + arrowprops=None, force_text=(0.1, 0.2), + expand_points=(1.2, 1.2), save_steps=False + ) + except ImportError: + pass + + +@dataclass +class MapStyleConfig: + """ + Encapsulates sizing, text fonts, and styling properties for rendering. + """ + font_size: int = 18 + city_font_size: int = 10 + legend_font_size: int = 10 + title_font_size: int = 20 + figsize: tuple[float, float] = (10, 10) + region_alpha: float = 0.7 + legend_digits: int = 0 + + +@dataclass +class MapDataElements: + """ + Groups optional geographical annotations, basemaps, and layout limits. + """ + plot_title: str | None = None + legend_title: str | None = None + cities: dict[str, tuple[float, float]] | None = None + x_limits: tuple[float, float] | None = None + y_limits: tuple[float, float] | None = None + basemap_path: str | Path | Any | None = None + epicenter: tuple[float, float] | None = None + + +def plot_variable(df, admin_boundaries, column, classifier, colors, *, + elements: MapDataElements = None, + style: MapStyleConfig = None): + """ + Plot a classified geospatial variable with optional basemap + and annotations. + """ + import matplotlib.pyplot as plt + import matplotlib.patches as mpatches + + # Fallback to defaults if no custom styles/elements are passed + style = style or MapStyleConfig() + elements = elements or MapDataElements() + + # Pre-process using the helper functions we built previously + df, colors = _prepare_classified_data( + df, admin_boundaries, column, classifier, colors) + labels = _build_legend_labels(classifier, style.legend_digits) + + fig, ax = plt.subplots(figsize=style.figsize) + _overlay_basemap(ax, elements.basemap_path, df.crs) + + # Plot each class with its corresponding color + for i, color in enumerate(colors): + subset = df[df["class"] == i] + if not subset.empty: + subset.plot(ax=ax, color=color, edgecolor="none", + alpha=style.region_alpha) + + epicenter_handle = None + if elements.epicenter is not None: + lon, lat = elements.epicenter + epicenter_handle = ax.scatter( + lon, lat, marker='*', s=150, color='yellow', + edgecolor='black', linewidth=1, zorder=10, label='Epicenter' + ) + + # Create legend handles + handles = [mpatches.Patch(color=col, label=lab) + for col, lab in zip(colors, labels)] + if epicenter_handle is not None: + handles.append(epicenter_handle) + + ax.legend(handles=handles, title=elements.legend_title, framealpha=0.7, + title_fontsize=style.font_size, + fontsize=style.legend_font_size, loc="best") + + admin_boundaries.plot(ax=ax, alpha=0.4, edgecolor="black", + facecolor="none", linewidth=0.4) + _resolve_limits(ax, elements.x_limits, elements.y_limits, + elements.epicenter) + _overlay_cities(ax, elements.cities, style.city_font_size) + + if elements.plot_title: + ax.set_title(elements.plot_title, fontsize=style.title_font_size) + + ax.set_xlabel("Longitude", fontsize=style.font_size) + ax.set_ylabel("Latitude", fontsize=style.font_size) + ax.set_aspect('equal', adjustable='box') + fig.tight_layout() + return fig, ax diff --git a/openquake/calculators/postrisk/make_impact_reports.py b/openquake/calculators/postrisk/make_impact_reports.py new file mode 100644 index 000000000000..5685e04d5ac7 --- /dev/null +++ b/openquake/calculators/postrisk/make_impact_reports.py @@ -0,0 +1,264 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# vim: tabstop=4 shiftwidth=4 softtabstop=4 +# +# Copyright (C) 2026 GEM Foundation +# +# OpenQuake is free software: you can redistribute it and/or modify it +# under the terms of the GNU Affero General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# OpenQuake is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with OpenQuake. If not, see . + +import pathlib +import tempfile +import logging +from datetime import datetime, timezone +import pandas as pd +from openquake.baselib import config, sap +from openquake.calculators.country_impact_report_builder import ( + CountryImpactReportBuilder) +from openquake.calculators.extract import extract +from openquake.calculators.country_impact_report_utils import ( + EventContext, ReportOptions, LOSS_METADATA) +from openquake.commonlib import logs +from openquake.commonlib.readinput import get_close_countries + +cd = pathlib.Path(__file__).parent + +LOSS_LABELS = [v["label"] for v in LOSS_METADATA.values()] + +DISCLAIMER_TXT = ''' + This is an automatically generated draft. Content has not been verified for + accuracy by a human reviewer. The metrics presented were estimated based on + ground shaking information from ShakeMap only. Impact assessments are + subject to changes as more information becomes available.''' + + +def _get_impact_summary_data(dstore, iso3, no_uncertainty): + aggrisk_tags = extract(dstore, 'aggrisk_tags') + mapping = { + meta["label"]: loss_type + for loss_type, meta in LOSS_METADATA.items() + } + rows = aggrisk_tags.loc[ + (aggrisk_tags['ID_0'] == iso3) & + (aggrisk_tags['loss_type'].isin(mapping.values())) + ] + if rows.empty: + logging.info( + f"No losses estimated for country {iso3}. Skipping report") + return None + loss_threshold = 1 + if all(r.lossmea < loss_threshold for _, r in rows.iterrows()): + logging.info(f"Estimated losses for country {iso3} are negligible" + f" (all lossmea < {loss_threshold}). Skipping report.") + return None + summary_data = {} + for label, lt in mapping.items(): + matching_rows = rows.loc[rows['loss_type'] == lt] + if not matching_rows.empty: + r = matching_rows.iloc[0] + q50 = int(round(r.get('q50', 0))) + q05 = int(round(r.get('q05', 0))) + q95 = int(round(r.get('q95', 0))) + exposed_val = int(round(r.get('value', 0))) + else: + q50 = q05 = q95 = exposed_val = 0 + # Format with thousands separators + if no_uncertainty: + # Display only the central value + summary_data[label] = f"{q50:,}" + else: + # Display the range + summary_data[label] = f"{q05:,} - {q95:,}" + summary_data[f"{label}_exposed"] = f"{exposed_val:,}" + return summary_data + + +def make_report_for_country( + iso3, adm_level, event, options, losses_df, summary_data, + dstore, time_of_calc, oqparam): + builder = CountryImpactReportBuilder( + iso3, adm_level, event, options, losses_df, summary_data, + dstore, time_of_calc, oqparam) + builder.build() + + +def to_utc_string(ts: str) -> str: + """ + Convert a timestamp with timezone offset (e.g. '+08:00') + to the format: 'YYYY-MM-DD HH:MM:SS UTC' + """ + if not ts: + return "unknown" + dt = datetime.fromisoformat(ts) + if dt.tzinfo is None: + logging.warning("Timestamp has no timezone information") + return dt.strftime('%Y-%m-%d %H:%M:%S') + dt_utc = dt.astimezone(timezone.utc) + return dt_utc.strftime('%Y-%m-%d %H:%M:%S') + ' UTC' + + +def get_dynamic_threshold(mag): + """ + Returns a search threshold in degrees based on earthquake magnitude. + """ + if mag < 5.0: + return 1.0 # ~111 km + elif mag < 6.5: + return 2.0 # ~222 km + elif mag < 7.5: + return 3.0 # ~333 km + else: + return 5.0 # ~555 km + + +def _open_dstore(dstore): + """ + Resolve the dstore argument (path/int from CLI, or an open Datastore), + returning (dstore, calc_id). + """ + if isinstance(dstore, (str, int)): + # NOTE: called from the command line + from openquake.commonlib import datastore + calc_id = int(dstore) + dstore = datastore.read(calc_id, mode='r+') + else: + calc_id = dstore.calc_id + return dstore, calc_id + + +def _get_basemap_path(): + try: + return config.directory.basemap_file + except AttributeError: + logging.error('config.directory.basemap_file is missing!') + return None + + +def _is_no_uncertainty(oqparam): + # If the ground motion is fully deterministic, we suppress uncertainty + # ranges in the report and show only the central (point) estimate. + return (oqparam.number_of_ground_motion_fields == 1 + and abs(oqparam.truncation_level) < 1e-8) + + +def _get_losses_df(avg_losses): + """ + Use the median (quantile-0.5) as the representative point estimate for + the spatial loss maps, consistent with how _get_impact_summary_data + displays the central value. Fall back to the mean only if the median is + unavailable (e.g. a calculation run without quantile outputs). + """ + if (hasattr(avg_losses, 'quantile-0.5') + and avg_losses['quantile-0.5'] is not None): + return pd.DataFrame(avg_losses['quantile-0.5']), 'Median' + elif hasattr(avg_losses, 'mean') and avg_losses.mean is not None: + logging.warning( + "Median losses not available; falling back to mean for " + "loss maps.") + return pd.DataFrame(avg_losses.mean), 'Mean' + else: + raise RuntimeError( + "avg_losses has neither 'quantile' nor 'mean' attribute; " + "cannot build losses DataFrame.") + + +def _get_event_name(rupdic): + try: + return rupdic['description'] + except KeyError: + return rupdic['title'] + + +def _get_shakemap_version(rupdic): + try: + return rupdic['shakemap_desc'] + except KeyError: + return None + + +def _get_threshold_deg(threshold_deg, mag): + if threshold_deg is None: + threshold_deg = get_dynamic_threshold(mag) + logging.info(f"Magnitude {mag} detected. Using dynamic" + f" threshold: {threshold_deg} degrees.") + return threshold_deg + return float(threshold_deg) + + +def _build_report_contexts(dstore, oqparam, calc_id, threshold_deg): + """ + Gather everything needed to build per-country reports and return + (event_ctx, report_opts, losses_df, iso3_codes, time_of_calc). + """ + mag = oqparam.rupture_dict['mag'] + lon = oqparam.rupture_dict['lon'] + lat = oqparam.rupture_dict['lat'] + no_uncertainty = _is_no_uncertainty(oqparam) + + avg_losses = extract(dstore, 'avg_losses?kind=stats') + losses_df, loss_metric = _get_losses_df(avg_losses) + + rupdic = oqparam.rupture_dict + event_name = _get_event_name(rupdic) + # FIXME: do we prefer to show UTC or perhaps it is more intuitive + # to show the local time? + event_date = to_utc_string(oqparam.local_timestamp) + shakemap_version = _get_shakemap_version(rupdic) + + job = logs.dbcmd('get_job', calc_id) + time_of_calc = job.start_time.strftime('%Y-%m-%d %H:%M:%S') + ' UTC' + + threshold_deg = _get_threshold_deg(threshold_deg, mag) + # close countries are ordered by ascending distance + iso3_codes = get_close_countries(lon, lat, buffer_radius=threshold_deg) + if not iso3_codes: + raise RuntimeError( + "No country within {threshold_deg} from the hypocenter") + + event_ctx = EventContext( + name=event_name, date=event_date, hypocenter=(lon, lat), + shakemap_version=shakemap_version) + report_opts = ReportOptions( + disclaimer_txt=DISCLAIMER_TXT, + basemap_path=_get_basemap_path(), threshold_deg=threshold_deg, + no_uncertainty=no_uncertainty, loss_metric=loss_metric) + + return event_ctx, report_opts, losses_df, iso3_codes, time_of_calc + + +def main(dstore, adm_level=1, threshold_deg=None): + """ + Create an impact report in PDF and PNG formats + """ + dstore, calc_id = _open_dstore(dstore) + adm_level = int(adm_level) + + dstore.close() + dstore.open('r+') + dstore.export_dir = config.directory.custom_tmp or tempfile.gettempdir() + oqparam = dstore['oqparam'] + + event_ctx, report_opts, losses_df, iso3_codes, time_of_calc = ( + _build_report_contexts(dstore, oqparam, calc_id, threshold_deg)) + + for iso3 in iso3_codes: + summary_data = _get_impact_summary_data( + dstore, iso3, report_opts.no_uncertainty) + if summary_data is not None: + make_report_for_country( + iso3, adm_level, event_ctx, report_opts, + losses_df, summary_data, dstore, time_of_calc, oqparam) + + +if __name__ == '__main__': + sap.run(main) diff --git a/openquake/calculators/views.py b/openquake/calculators/views.py index 734a3a12e865..580ce36f430a 100644 --- a/openquake/calculators/views.py +++ b/openquake/calculators/views.py @@ -248,13 +248,14 @@ def text_table(data, header=None, fmt=None, ext='rst'): lines.append(sepline) return '\n'.join(lines) + @view.add('high_hazard') def view_high_hazard(token, dstore): """ Returns the sites with hazard curve below max(poes) """ oq = dstore['oqparam'] - max_poe= max(oq.poes) + max_poe = max(oq.poes) max_hazard = dstore.sel('hcurves-stats', stat='mean', lvl=0)[:, 0, :, 0] # NSML1 -> NM high = (max_hazard > max_poe).all(axis=1) @@ -715,6 +716,7 @@ def view_required_params_per_trt(token, dstore): return text_table(tbl, header='trt gsims req_params'.split(), fmt=scientificformat) + def discard_small(values): """ Discard values 10x smaller than the mean @@ -807,7 +809,6 @@ def view_task_eb(token, dstore): return msg - def view_task(token, dstore, taskname): """ Display info about a given task. Here are a few examples of usage:: @@ -818,7 +819,8 @@ def view_task(token, dstore, taskname): _, index = token.split(':') if 'source_data' not in dstore: return 'Missing source_data' - data = get_array(dstore['task_info'][()], taskname=taskname.encode('ascii')) + data = get_array(dstore['task_info'][()], + taskname=taskname.encode('ascii')) if len(data) == 0: raise RuntimeError('No task_info for classical') data.sort(order='duration') @@ -1526,6 +1528,7 @@ def view_branchsets(token, dstore): return text_table(enumerate(map(repr, clt.branchsets)), header=['bsno', 'bset'], ext='org') + @view.add('sm_rlzs') def view_sm_rlzs(token, dstore): """ @@ -1540,10 +1543,12 @@ def view_sm_rlzs(token, dstore): else: sm_rlzs = dstore['full_lt'].sm_rlzs header = ['ordinal', 'lt_path', 'value', 'samples', 'weight'] + def row(rlz): value = ast.literal_eval(rlz.value.decode('utf8')) return (rlz.ordinal, '_'.join(rlz.lt_path), value, rlz.samples, rlz.weight) + return text_table(map(row, sm_rlzs), header, ext='org') @@ -1911,6 +1916,7 @@ def view_long_ruptures(token, dstore): arr.sort(order='maxlen') return arr + @view.add('msr') def view_msr(token, dstore): dic = dstore['_csm'].get_msr_by_grp() diff --git a/openquake/commonlib/calc.py b/openquake/commonlib/calc.py index 674a3a6eb58c..f8d5908c87f4 100644 --- a/openquake/commonlib/calc.py +++ b/openquake/commonlib/calc.py @@ -293,7 +293,8 @@ def _save_events(self, rup_array, idx_start_stop): trt_smr = 0 rlzs = numpy.concatenate( list(rlzs_by_gsim[trt_smr].values()), dtype=U32) - records = get_events(rup_array[start:stop], rlzs, self.scenario) + records = get_events( + rup_array[start:stop], rlzs, self.scenario) nr = len(records) events[i:i + nr] = records # (id, rup_id, rlz_id) i += nr diff --git a/openquake/commonlib/oqvalidation.py b/openquake/commonlib/oqvalidation.py index da8435e49fef..69c0b6d54ff1 100644 --- a/openquake/commonlib/oqvalidation.py +++ b/openquake/commonlib/oqvalidation.py @@ -506,6 +506,12 @@ Example: *master_seed = 1234*. Default: 123456789 +make_impact_reports: + Produce a one-page impact reports for each affected country, + as postrisk_func + Example: *make_impact_reports = true*. + Default: False + max: Compute the maximum across realizations. Akin to mean and quantiles. Example: *max = true*. @@ -600,6 +606,11 @@ Example: *mosaic_model = ZAF* Default: empty string +notes: + Additional information about the job + Example: 'Lorem ipsum' + Default: None + num_epsilon_bins: Number of epsilon bins in disaggregation calculations. Example: *num_epsilon_bins = 3*. @@ -1197,6 +1208,7 @@ class OqParam(valid.ParamSet): maximum_distance = valid.Param(valid.IntegrationDistance.new) # km maximum_distance_stations = valid.Param(valid.positivefloat, None) # km asset_hazard_distance = valid.Param(valid.floatdict, {'default': 15}) # km + make_impact_reports = valid.Param(valid.boolean, False) max = valid.Param(valid.boolean, False) max_data_transfer = valid.Param(valid.positivefloat, 2E11) max_nodes_network = valid.Param(valid.positiveint, 1000) @@ -1211,6 +1223,7 @@ class OqParam(valid.ParamSet): minimum_intensity = valid.Param(valid.floatdict, {}) # IMT -> minIML minimum_magnitude = valid.Param(valid.floatdict, {'default': 0}) # by TRT modal_damage_state = valid.Param(valid.boolean, False) + notes = valid.Param(valid.utf8, None) number_of_ground_motion_fields = valid.Param(valid.positiveint) number_of_logic_tree_samples = valid.Param(valid.positiveint, 0) num_epsilon_bins = valid.Param(valid.positiveint, 1) @@ -1664,7 +1677,7 @@ def set_loss_types(self): with datastore.read(self.hazard_calculation_id) as ds: self._parent = ds['oqparam'] if not self.total_losses: - self.total_losses = self._parent.total_losses + self.total_losses = self._parent.total_losses else: self._parent = None # set all_cost_types diff --git a/openquake/commonlib/readinput.py b/openquake/commonlib/readinput.py index 042829e8d4c7..78747867905c 100644 --- a/openquake/commonlib/readinput.py +++ b/openquake/commonlib/readinput.py @@ -127,9 +127,10 @@ def get_close_mosaic_models(lon, lat, buffer_radius): This distance is in the same units as the point's coordinates (i.e. degrees), and it defines how far from the point the buffer should extend in all directions, - creating a circular buffer region around the point - :returns: list of mosaic models intersecting the circle - centered on the given coordinates having the specified radius + creating a circular buffer region around the point. + :returns: list of mosaic models intersecting the circle, + centered on the given coordinates having the specified radius. + Models are ordered by ascending distance. """ mosaic_df = read_mosaic_df() close_mosaic_models = geo.utils.geolocate_within_buffer( @@ -144,7 +145,7 @@ def get_close_mosaic_models(lon, lat, buffer_radius): return close_mosaic_models -def get_closest_country(lon, lat, buffer_radius): +def get_close_countries(lon, lat, buffer_radius): """ :param lon: longitude :param lat: latitude @@ -152,12 +153,28 @@ def get_closest_country(lon, lat, buffer_radius): This distance is in the same units as the point's coordinates (i.e. degrees), and it defines how far from the point the buffer should extend in all directions, - creating a circular buffer region around the point - :returns: the iso3 code of the closest country or '???' + creating a circular buffer region around the point. + :returns: the iso3 codes of the close countries, ordered by + ascending distance. """ countries_df = read_countries_df() close_countries = geo.utils.geolocate_within_buffer( lon, lat, buffer_radius, countries_df) + return close_countries + + +def get_closest_country(lon, lat, buffer_radius): + """ + :param lon: longitude + :param lat: latitude + :param buffer_radius: radius of the buffer around the point. + This distance is in the same units as the point's + coordinates (i.e. degrees), and it defines how far from + the point the buffer should extend in all directions, + creating a circular buffer region around the point + :returns: the iso3 code of the closest country or '???' + """ + close_countries = get_close_countries(lon, lat, buffer_radius) if not close_countries: return '???' # close_countries are ordered by ascending distance @@ -933,7 +950,7 @@ def get_rupture(oqparam): # converting rupture_model from json to an oq-compatible xml rupture_model = convert_to_oq_xml(rupture_model, rupture_model) # NB: this is tested in aristotle_run - if rupture_model and rupture_model.endswith('.xml'): + elif rupture_model and rupture_model.endswith('.xml'): [rup_node] = nrml.read(rupture_model) conv = sourceconverter.RuptureConverter(oqparam.rupture_mesh_spacing) rup = conv.convert_node(rup_node) @@ -1901,10 +1918,15 @@ def read_countries_df(): """ :returns: a DataFrame of geometries for the world countries """ - logging.info('Reading geoBoundariesCGAZ_ADM0.gpkg') # slow - fname = os.path.join(os.path.dirname(global_risk.__file__), - 'geoBoundariesCGAZ_ADM0.gpkg') - return read_geometries(fname, 'shapeGroup') + country_boundaries_file = None + if hasattr(config.directory, 'admin0_boundaries_file'): + country_boundaries_file = config.directory.admin0_boundaries_file + if not country_boundaries_file: + country_boundaries_file = os.path.join( + os.path.dirname(global_risk.__file__), + 'geoBoundariesCGAZ_ADM0.gpkg') + logging.info(f'Reading {country_boundaries_file}') + return read_geometries(country_boundaries_file, 'shapeGroup') def read_cities_df(lon_field='longitude', lat_field='latitude', diff --git a/openquake/commonlib/tests/readinput_test.py b/openquake/commonlib/tests/readinput_test.py index ea4df83c3a10..0e22a2c293db 100644 --- a/openquake/commonlib/tests/readinput_test.py +++ b/openquake/commonlib/tests/readinput_test.py @@ -28,6 +28,8 @@ from openquake.hazardlib.calc.filters import MINMAG, MAXMAG from openquake.risklib import asset from openquake.commonlib import readinput, datastore +from openquake.commonlib.readinput import ( + get_close_mosaic_models, get_close_countries) from openquake.qa_tests_data.logictree import ( case_02, case_15, case_21, case_25) from openquake.qa_tests_data.classical import case_34, case_65 @@ -530,7 +532,8 @@ def test_read_station_data(self): oq = readinput.get_oqparam(os.path.join(DATADIR, 'job.ini')) sitecol = readinput.get_site_collection(oq) with self.assertRaises(InvalidFile) as ctx: - readinput.get_station_data(oq, sitecol, duplicates_strategy='error') + readinput.get_station_data( + oq, sitecol, duplicates_strategy='error') self.assertIn( "Stations_NIED.csv: has duplicate sites ['GIF001', 'GIF013']", str(ctx.exception)) @@ -550,7 +553,8 @@ def test_read_station_data(self): oq.inputs['station_data'], 'LONGITUDE', 'LATITUDE', 'STATION_ID', duplicates_strategy='avg') self.assertTrue('GIF001|GIF013' in df['STATION_ID'].values) - pga_avg = df[df['STATION_ID'] == 'GIF001|GIF013']['PGA_VALUE'].values[0] + pga_avg = df[df['STATION_ID'] == 'GIF001|GIF013'][ + 'PGA_VALUE'].values[0] # using the same mean operator used in read_df and expecting the same # approximation data = {'values': [pga_first, pga_last]} @@ -563,8 +567,10 @@ class ReadSourceModelsTestCase(unittest.TestCase): def test(self): base = os.path.dirname(case_65.__file__) hdf5path = general.gettemp(suffix='.hdf5') - fnames = [os.path.join(base, 'ssm.xml'), os.path.join(base, 'sections.xml')] - smodels = readinput.read_source_models(fnames, hdf5path, investigation_time=1.) + fnames = [os.path.join(base, 'ssm.xml'), + os.path.join(base, 'sections.xml')] + smodels = readinput.read_source_models( + fnames, hdf5path, investigation_time=1.) nrups = 0 for smodel in smodels: for sg in smodel.src_groups: @@ -572,3 +578,17 @@ def test(self): for rup in src.iter_ruptures(): nrups += 1 self.assertEqual(nrups, 3) + + +class GetCloseRegionsTestCase(unittest.TestCase): + def test_get_close_mosaic_models(self): + lon, lat = 124.0, 8.5 + mosaic_models = get_close_mosaic_models(lon, lat, buffer_radius=5) + self.assertEqual(mosaic_models, ['PHL', 'IDN', 'OPA']) + mosaic_models = get_close_mosaic_models(lon, lat, buffer_radius=0.5) + self.assertEqual(mosaic_models, ['PHL']) + + def test_get_close_countries(self): + lon, lat = 124.0, 8.5 + countries = get_close_countries(lon, lat, buffer_radius=0.5) + self.assertEqual(countries, ['PHL']) diff --git a/openquake/engine/impact.py b/openquake/engine/impact.py index 43b3f09bc450..45db3aec1cc6 100644 --- a/openquake/engine/impact.py +++ b/openquake/engine/impact.py @@ -70,6 +70,7 @@ def main_cmd(usgs_id, rupture_file=None, maximum_distance_stations='', msr='WC1994', approach='use_shakemap_from_usgs', loglevel='warn', + make_impact_reports=False, userlevel=1): # with userlevel=1 use shakemap, else rupture """ This script is meant to be called from the command-line @@ -140,6 +141,7 @@ def main_cmd(usgs_id, rupture_file=None, main_cmd.msr = 'Magnitude scaling relationship' main_cmd.approach = 'For instance use_shakemap_from_usgs' main_cmd.loglevel = 'Log level' +main_cmd.make_impact_reports = 'Make one-page report for each involved country' main_cmd.userlevel = 'User level' main = main_cmd diff --git a/openquake/hazardlib/shakemap/parsers.py b/openquake/hazardlib/shakemap/parsers.py index 69e1bc76ba46..b9af62ebfdfc 100644 --- a/openquake/hazardlib/shakemap/parsers.py +++ b/openquake/hazardlib/shakemap/parsers.py @@ -395,9 +395,12 @@ def convert_to_oq_xml(input_json_file, output_xml_file): def utc_to_local_time(utc_timestamp, lon, lat): """ - Convert a timestamp '%Y-%m-%dT%H:%M:%S.%fZ' into a datetime object + Convert a timestamp string or a datetime into a local datetime object """ - utc_time = datetime.strptime(utc_timestamp, '%Y-%m-%dT%H:%M:%S.%fZ') + if isinstance(utc_timestamp, str): + utc_time = datetime.strptime(utc_timestamp, '%Y-%m-%dT%H:%M:%S.%fZ') + else: + utc_time = utc_timestamp try: from timezonefinder import TimezoneFinder except ImportError: @@ -637,12 +640,14 @@ def load_rupdic_from_finite_fault(usgs_id, mag, products): ff = _get_usgs_preferred_item(products['finite-fault']) p = ff['properties'] - # TODO: we probably need to get the rupture coordinates from shakemap_polygon.txt + # TODO: we probably need to get the rupture coordinates from + # shakemap_polygon.txt # if 'shakemap_polygon.txt' in ff['contents']: # # with open(f'/tmp/{usgs_id}-shakemap_polygon.txt', 'wb') as f: # # f.write(urlopen(url).read()) # if user.testdir: # in parsers_test - # fname = os.path.join(user.testdir, f'{usgs_id}-shakemap_polygon.txt') + # fname = os.path.join(user.testdir, + # f'{usgs_id}-shakemap_polygon.txt') # text = open(fname).read() # else: # url = ff['contents']['shakemap_polygon.txt']['url'] @@ -758,12 +763,20 @@ def download_shakemap_rupture_data(usgs_id, shakemap_contents, user): def extract_event_details(ffm): # Extract event details from the geojson metadata # and return a data object with the event details - epicenter = ffm["metadata"]["epicenter"] + try: + hypocenter = ffm["metadata"]["hypocenter"] + except KeyError: + # NOTE: the field was originally named 'epicenter'. However, since + # it included also the depth, we assume it was always meant to + # represent the hypocenter. For compatibility with old events + # providing data in the old format, we look for the old name if the new + # one is missing. + hypocenter = ffm["metadata"]["epicenter"] return { - "mag": epicenter.get("mag"), - "lon": epicenter.get("lon"), - "lat": epicenter.get("lat"), - "dep": epicenter.get("depth"), + "mag": hypocenter.get("mag"), + "lon": hypocenter.get("lon"), + "lat": hypocenter.get("lat"), + "dep": hypocenter.get("depth"), } @@ -923,14 +936,18 @@ def download_mmi(usgs_id, shakemap_contents, user): return mmi_file -def convert_rup_data(rup_data, usgs_id, rup_path, shakemap_array=None): +def convert_rup_data(rup_data, usgs_id, rup_path, utc_time_ms, + shakemap_array=None): """ Convert JSON data coming from the USGS into a rupdic """ md = rup_data['metadata'] lon = md['lon'] lat = md['lat'] - local_time = utc_to_local_time(md['time'], lon, lat) + # NOTE: retrieving the local timestamp from rup_data['metadata'] is not + # reliable + utc_time = datetime.fromtimestamp(utc_time_ms / 1000.0, tz=timezone.utc) + local_time = utc_to_local_time(utc_time, lon, lat) time_event = local_time_to_time_event(local_time) return { 'lon': lon, 'lat': lat, 'dep': md['depth'], @@ -989,7 +1006,13 @@ def _contents_properties_shakemap(usgs_id, user, get_grid, monitor, return None, None, None, None, err properties = usgs_event_data['properties'] - shakemaps = properties['products']['shakemap'] + + try: + shakemaps = properties['products']['shakemap'] + except KeyError: + err = {'status': 'failed', + 'error_msg': f'No ShakeMap available for {usgs_id}'} + return None, None, None, None, err if shakemap_version == 'usgs_preferred': shakemap = _get_usgs_preferred_item(shakemaps) else: @@ -1155,8 +1178,8 @@ def _get_rup_dic_from_csv(usgs_id, user, rupture_file): def get_stations_from_usgs(usgs_id, user=User(), - monitor=performance.Monitor(), - shakemap_version='usgs_preferred'): + monitor=performance.Monitor(), + shakemap_version='usgs_preferred'): n_stations = 0 try: usgs_id = valid.simple_id(usgs_id) @@ -1175,8 +1198,8 @@ def get_stations_from_usgs(usgs_id, user=User(), def ms_to_utc_date_time(ms): - # convert from milliseconds to utc date time - dt = datetime.fromtimestamp(ms / 1000, timezone.utc) + # convert from milliseconds to timezone-aware UTC date time + dt = datetime.fromtimestamp(ms / 1000, timezone.utc) # convert to seconds return dt.strftime("%Y-%m-%d %H:%M:%S") @@ -1204,7 +1227,12 @@ def get_shakemap_versions(usgs_id, user=User(), js = json.loads(text) properties = js['properties'] - shakemaps = properties['products']['shakemap'] + try: + shakemaps = properties['products']['shakemap'] + except KeyError: + err = {'status': 'failed', + 'error_msg': f'No ShakeMap available for {usgs_id}'} + return None, None, err usgs_preferred_shakemap = _get_usgs_preferred_item(shakemaps) usgs_preferred_version = usgs_preferred_shakemap['id'] sorted_shakemaps = sorted( @@ -1260,57 +1288,105 @@ def make_rup_from_dic(inputdic, rupture_file): return rup, rupdic, rupture_issue -def _fetch_usgs_rupture(inputdic, contents, properties, rup_data, - user, monitor): - """ - Download rupture or finite fault model files from USGS if needed. - """ - usgs_id = inputdic['usgs_id'] - approach = inputdic['approach'] +def _load_rupdic_for_approach( + usgs_id, approach, inputdic, properties, contents): + # Build the initial rupdic for approaches that derive it directly from + # USGS origin/finite-fault data, before any rupture-file download step. rupdic = {} - rup = None - rupture_issue = None - + err = {} if approach in ['use_pnt_rup_from_usgs', 'build_rup_from_usgs']: - if inputdic.get('lon') is None: + if inputdic.get('lon') is None: # don't override user-inserted values rupdic, err = load_rupdic_from_origin( usgs_id, properties['products']) - if err: - return None, None, None, None, err else: rupdic = inputdic.copy() elif 'download/rupture.json' not in contents: + # happens for us6000f65h in parsers_test rupdic, err = load_rupdic_from_finite_fault( usgs_id, properties['mag'], properties['products']) - if err: - return None, None, None, None, err - - if not rup_data and approach in ['use_shakemap_from_usgs', - 'use_shakemap_fault_rup_from_usgs', - 'use_finite_fault_model_from_usgs']: - if approach == 'use_finite_fault_model_from_usgs': - with monitor('Download finite fault rupture'): - rupture_file, err = download_finite_fault_rupture( - usgs_id, user, monitor) - if err: - return None, None, None, None, err - else: - with monitor('Downloading rupture json'): - rup_data, rupture_file = download_shakemap_rupture_data( - usgs_id, contents, user) - if rupture_file: - rup, rupdic, updated_rup_data, rupture_issue = ( - _convert_rupture_file(inputdic, rupture_file, usgs_id, user) - ) - if updated_rup_data: - rup_data = updated_rup_data - elif approach in ['use_shakemap_fault_rup_from_usgs', - 'use_finite_fault_model_from_usgs']: - err = {"status": "failed", - "error_msg": 'Unable to retrieve rupture geometries'} - return None, None, None, None, err + return rupdic, err + - return rup, rupdic, rup_data, rupture_issue, {} +def _download_rup_data(usgs_id, approach, inputdic, contents, user, monitor): + # Download rupture geometry for approaches that need it from USGS + # (shakemap, shakemap+fault rupture, or finite fault model) and convert + # it via _convert_rupture_file when a rupture file is obtained. + # Only called when approach is one of: + # 'use_shakemap_from_usgs', 'use_shakemap_fault_rup_from_usgs', + # 'use_finite_fault_model_from_usgs' + rup = None + rupdic = {} + rup_data = {} + rupture_issue = None + err = {} + + if approach == 'use_finite_fault_model_from_usgs': + with monitor('Download finite fault rupture'): + rupture_file, err = download_finite_fault_rupture( + usgs_id, user, monitor) + if err: + return None, {}, {}, None, err + else: # use_shakemap_from_usgs or use_shakemap_fault_rup_from_usgs + with monitor('Downloading rupture json'): + rup_data, rupture_file = download_shakemap_rupture_data( + usgs_id, contents, user) + + if rupture_file: + rup, rupdic, updated_rup_data, rupture_issue = _convert_rupture_file( + inputdic, rupture_file, usgs_id, user) + if updated_rup_data: + rup_data = updated_rup_data + elif approach in ['use_shakemap_fault_rup_from_usgs', + 'use_finite_fault_model_from_usgs']: + err = {"status": "failed", + "error_msg": 'Unable to retrieve rupture geometries'} + + return rup, rupdic, rup_data, rupture_issue, err + + +def _finalize_rupdic(rupdic, rup_data, usgs_id, rupture_file, shakemap, + shakemap_desc, contents, inputdic, approach, properties, + user): + # Merge downloaded/converted rupture data into rupdic, backfill any + # fields still missing from inputdic, and attach metadata (mmi file, + # shakemap array, title, shakemap description). + # Returns a new dict + new_rupdic = dict(rupdic) + utc_time_ms = properties['time'] + if rup_data: + converted_rup_data = convert_rup_data( + rup_data, usgs_id, rupture_file, utc_time_ms, shakemap) + if 'rupture_file' in new_rupdic: # already converted: do not overwrite + converted_rup_data.pop('rupture_file') + new_rupdic.update(converted_rup_data) + if 'local_timestamp' not in new_rupdic: + utc_time = datetime.fromtimestamp(utc_time_ms / 1000.0, + tz=timezone.utc) + local_time = utc_to_local_time( + utc_time, new_rupdic['lon'], new_rupdic['lat']) + new_rupdic['local_timestamp'] = str(local_time) + new_rupdic['time_event'] = local_time_to_time_event(local_time) + for key in inputdic: + if inputdic[key] is not None and key not in new_rupdic: + new_rupdic[key] = inputdic[key] + + if 'mmi_file' not in new_rupdic: + new_rupdic['mmi_file'] = download_mmi(usgs_id, contents, user) + if approach == 'use_shakemap_from_usgs': + new_rupdic['shakemap_array'] = shakemap + new_rupdic['title'] = properties['title'] + new_rupdic['shakemap_desc'] = shakemap_desc + + return new_rupdic + + +def _build_planar_rupture(rupdic): + # Try to build a planar rupture from rupdic. + # Returns (rup, error_msg) where error_msg is None on success. + try: + return build_planar_rupture_from_dict(rupdic), None + except ValueError as exc: + return None, str(exc) def get_rup_dic(inputdic, user=User(), use_shakemap=False, @@ -1322,20 +1398,37 @@ def get_rup_dic(inputdic, user=User(), use_shakemap=False, NOTE: this function is called twice by impact_validate: first when retrieving rupture data, then when running the job. + + :param inputdic: + dictionary with ShakeMap ID and other parameters + :param user: + User instance + :param use_shakemap: + download the ShakeMap only if True + :param shakemap_version: + id of the ShakeMap to be used (if the ShakeMap is used) + :param rupture_file: + None + :returns: + (rupture object or None, rupture dictionary, error dictionary or {}) """ rupdic = {} rup_data = {} + err = {} usgs_id = inputdic['usgs_id'] approach = inputdic['approach'] rup = None rupture_issue = None + if approach == 'provide_rup_params': return make_rup_from_dic(inputdic, rupture_file) + if rupture_file: rup, rupdic, rup_data, rupture_issue = _convert_rupture_file( inputdic, rupture_file, usgs_id, user) if rupture_issue or usgs_id == 'FromFile': return rup, rupdic, rupture_issue + assert usgs_id get_grid = user.level == 1 or use_shakemap contents, properties, shakemap, shakemap_desc, err = ( @@ -1345,54 +1438,50 @@ def get_rup_dic(inputdic, user=User(), use_shakemap=False, if err: return None, None, err - fetched_rup, fetched_rupdic, rup_data, fetched_issue, err = ( - _fetch_usgs_rupture( - inputdic, contents, properties, rup_data, user, monitor) - ) + rupdic, err = _load_rupdic_for_approach( + usgs_id, approach, inputdic, properties, contents) if err: return None, None, err - if fetched_rup: - rup = fetched_rup - if fetched_issue: - rupture_issue = fetched_issue - if fetched_rupdic: - rupdic.update(fetched_rupdic) + if not rup_data and approach not in ['use_pnt_rup_from_usgs', + 'build_rup_from_usgs']: + if approach in ['use_shakemap_from_usgs', + 'use_shakemap_fault_rup_from_usgs', + 'use_finite_fault_model_from_usgs']: + (rup, downloaded_rupdic, rup_data, rupture_issue, + err) = _download_rup_data( + usgs_id, approach, inputdic, contents, user, monitor) + if err: + return None, None, err + if downloaded_rupdic: + rupdic = downloaded_rupdic + + rupdic = _finalize_rupdic( + rupdic, rup_data, usgs_id, rupture_file, shakemap, shakemap_desc, + contents, inputdic, approach, properties, user) - if 'lon' not in rupdic: - rupdic = convert_rup_data(rup_data, usgs_id, rupture_file, shakemap) - for key in inputdic: - if inputdic[key] is not None and key not in rupdic: - rupdic[key] = inputdic[key] - if 'mmi_file' not in rupdic: - rupdic['mmi_file'] = download_mmi(usgs_id, contents, user) - if approach == 'use_shakemap_from_usgs': - rupdic['shakemap_array'] = shakemap - rupdic['title'] = properties['title'] - rupdic['shakemap_desc'] = shakemap_desc if not rup and not rup_data: # in parsers_test if approach == 'use_pnt_rup_from_usgs': rupdic['msr'] = 'PointMSR' - try: - rup = build_planar_rupture_from_dict(rupdic) - except ValueError as exc: - err = {"status": "failed", "error_msg": str(exc)} + rup, err_msg = _build_planar_rupture(rupdic) + if err_msg: + err = {"status": "failed", "error_msg": err_msg} return rup, rupdic, err elif (not rup and len(rup_data['features']) == 1 and rup_data['features'][0]['geometry']['type'] == 'Point'): # TODO: we can remove this when OQ can handle xml with Point ruptures rupdic['msr'] = 'PointMSR' - try: - rup = build_planar_rupture_from_dict(rupdic) - except ValueError as exc: - rupture_issue = {"status": "failed", "error_msg": str(exc)} + rup, err_msg = _build_planar_rupture(rupdic) + if err_msg: + rupture_issue = {"status": "failed", "error_msg": err_msg} + if rupture_issue and user.level > 1: # in parsers_test for us6000jllz # NOTE: hiding rupture-related issues to level 1 users rupdic['rupture_issue'] = rupture_issue['error_msg'] return rup, rupdic, err -# tested in the nightly tests aristotle_run +# tested in the nightly tests impact_run # the default argument is needed to avoid an # error in is_valid_shakemap def get_array_usgs_id(kind, id, contents={}): diff --git a/openquake/hazardlib/shakemap/validate.py b/openquake/hazardlib/shakemap/validate.py index e42a65c48c62..2db2ae466a5b 100644 --- a/openquake/hazardlib/shakemap/validate.py +++ b/openquake/hazardlib/shakemap/validate.py @@ -24,6 +24,7 @@ from openquake.baselib import config, general, hdf5, performance from openquake.hazardlib import valid from openquake.commonlib import readinput +from openquake.commonlib.readinput import get_close_mosaic_models from openquake.hazardlib.shakemap.parsers import get_rup_dic from openquake.qa_tests_data import mosaic from openquake.hazardlib.geo.utils import SiteAssociationError @@ -50,6 +51,8 @@ class ImpactParam: mosaic_model: str = None trt: str = None description: str = None + notes: str = None + make_impact_reports: bool = False def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): """ @@ -80,6 +83,7 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): params = dict( base_path='', # no .ini file description=self.description, + notes=self.notes, calculation_mode='scenario_risk', rupture_dict=str(rupdic), time_event=self.time_event, @@ -91,6 +95,7 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): self.number_of_ground_motion_fields), asset_hazard_distance=str(self.asset_hazard_distance), ses_seed=str(self.ses_seed), + make_impact_reports=bool(self.make_impact_reports), inputs=inputs) if use_shakemap: fname = general.gettemp(suffix='.npy') @@ -129,6 +134,7 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): 'rupture_file': 'Rupture model XML', 'use_shakemap': 'Use the ShakeMap', 'shakemap_version': 'ShakeMap version', + 'shakemap_desc': 'ShakeMap description', 'lon': 'Longitude (degrees)', 'lat': 'Latitude (degrees)', 'dep': 'Depth (km)', @@ -152,7 +158,9 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): 'nodal_plane': 'Nodal plane', 'msr': 'Magnitude scaling relationship', 'description': 'Description', + 'notes': 'Notes', 'no_uncertainty': 'No uncertainty', + 'make_impact_reports': 'Make one-page country reports', } IMPACT_FORM_PLACEHOLDERS = { @@ -182,6 +190,7 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): 'nodal_plane': '', 'msr': '', 'description': 'Leave blank to set automatically', + 'notes': '', } IMPACT_FORM_DEFAULTS = { @@ -253,6 +262,8 @@ def get_oqparams(self, usgs_id, mosaic_models, trts, use_shakemap): 'ses_seed': valid.positiveint, 'maximum_distance_stations': valid.positivefloat, 'description': valid.utf8, # if empty, it will be set automatically + 'notes': valid.utf8, + 'make_impact_reports': valid.boolean, } @@ -271,7 +282,7 @@ def _validate(POST): value = validation_func(POST.get(field)) except Exception as exc: blankable = ['dip', 'strike', 'maximum_distance_stations', - 'local_timestamp'] + 'local_timestamp', 'notes'] if field in blankable and POST.get(field) == '': if field in inputdic: inputdic[field] = None @@ -364,7 +375,7 @@ def impact_validate(POST, user, rupture_file=None, station_data_file=None, os.path.join(MOSAIC_DIR, 'exposure.hdf5')) with monitor('get_close_mosaic_models'): try: - mosaic_models = readinput.get_close_mosaic_models( + mosaic_models = get_close_mosaic_models( rupdic['lon'], rupdic['lat'], 5) except ValueError as exc: # e.g.: @@ -379,6 +390,8 @@ def impact_validate(POST, user, rupture_file=None, station_data_file=None, rupdic['rupture_was_loaded'] = rup is not None if 'description' in inputdic and inputdic['description']: params['description'] = inputdic['description'] + if 'notes' in inputdic: + params['notes'] = inputdic['notes'] if len(params) > 1: # called by impact_run params['rupture_dict'] = rupdic params['station_data_file'] = station_data_file diff --git a/openquake/server/db/actions.py b/openquake/server/db/actions.py index 1a21fc2e5616..878bf8c19c6c 100644 --- a/openquake/server/db/actions.py +++ b/openquake/server/db/actions.py @@ -536,7 +536,7 @@ def get_calcs(db, request_get_dict, allowed_users, user_acl_on=False, id=None): :returns: list of tuples (id, user_name, status, calculation_mode, is_running, description, pid, hazard_calculation_id, size_mb, - host, start_time, relevant) + host, start_time, relevant, tags) """ # helper to get job+calculation data from the oq-engine database query_params = [] diff --git a/openquake/server/settings.py b/openquake/server/settings.py index 69ebab350fcd..aa1e16bda3cf 100644 --- a/openquake/server/settings.py +++ b/openquake/server/settings.py @@ -287,20 +287,21 @@ CONTEXT_PROCESSORS.append('openquakeplatform.utils.oq_context_processor') TEMPLATES[0]['OPTIONS']['context_processors'] = CONTEXT_PROCESSORS -try: - # Try to load a local_settings.py from the current folder; this is useful - # when packages are used. A custom local_settings.py can be placed in - # /usr/share/openquake/engine, avoiding changes inside the python package - from local_settings import * # noqa -except ImportError: - # If no local_settings.py is availble in the current folder let's try to - # load it from openquake/server/local_settings.py +if not TEST: try: - from openquake.server.local_settings import * # noqa + # Try to load a local_settings.py from the current folder; this is useful + # when packages are used. A custom local_settings.py can be placed in + # /usr/share/openquake/engine, avoiding changes inside the python package + from local_settings import * # noqa except ImportError: - # If a local_setting.py does not exist - # settings in this file only will be used - pass + # If no local_settings.py is availble in the current folder let's try to + # load it from openquake/server/local_settings.py + try: + from openquake.server.local_settings import * # noqa + except ImportError: + # If a local_setting.py does not exist + # settings in this file only will be used + pass if SUPPRESS_PERMISSION_DENIED_WARNINGS: class SuppressPermissionDeniedWarnings(logging.Filter): diff --git a/openquake/server/static/css/impact.css b/openquake/server/static/css/impact.css index b6b769176adb..17e393bad09c 100644 --- a/openquake/server/static/css/impact.css +++ b/openquake/server/static/css/impact.css @@ -39,6 +39,7 @@ display: block; overflow-x: scroll; position: relative; + margin-top: 10px; } table#impact-losses { @@ -58,11 +59,11 @@ table#impact-losses th td { .shakemap-image-container { float: left; - /* If we want to adjust the size of images instead of making them of the suitable size + /* If we want to adjust the size of images instead of making them of the suitable size */ /* width: 50%; */ } -/* If we want to adjust the size of images instead of making them of the suitable size +/* If we want to adjust the size of images instead of making them of the suitable size */ /* .shakemap-image-container img{ */ /* width: 600px; */ /* } */ @@ -116,15 +117,35 @@ table#impact-losses th td { background-color: #f1f1f1; } -/* Hide all form elements except the submit button */ -.hidden_except_submit input:not([type="submit"]), -.hidden_except_submit select, -.hidden_except_submit textarea, -.hidden_except_submit button:not([type="submit"]), -.hidden_except_submit label { - display: none; -} - .hidden { display: none; } + +.grid-table { + border-collapse: collapse; + width: 55%; +} + +/* Horizontal separators for rows */ +.grid-table tr { + border-bottom: 1px solid #e0e0e0; +} + +.grid-table thead tr { + border-bottom: 2px solid #222; /* Stronger line under header */ +} + +/* Vertical separators for columns */ +.grid-table th:not(:last-child), +.grid-table td:not(:last-child) { + border-right: 1px solid #e0e0e0; +} + +.grid-table th, +.grid-table td { + padding: 2px 14px; +} + +.impact-inputs-container { + margin-top: 10px; +} diff --git a/openquake/server/static/js/impact.js b/openquake/server/static/js/impact.js index 4ba5f59ad900..52c72c185cd9 100644 --- a/openquake/server/static/js/impact.js +++ b/openquake/server/static/js/impact.js @@ -125,8 +125,11 @@ window.initImpactForm = function() { } function set_shakemap_version_selector() { + let shakemap_selector = $("#shakemap_version"); + shakemap_selector.empty(); const usgs_id = $.trim($("#usgs_id").val()); if (usgs_id == '') return; + $('input#usgs_id').prop('disabled', true); $('#submit_impact_get_rupture').prop('disabled', true); $('#getStationDataFromUsgs').prop('disabled', true); $('#submit_impact_calc').prop('disabled', true); @@ -166,6 +169,7 @@ window.initImpactForm = function() { var err_msg = resp.error_msg; diaerror.show(false, "Error", err_msg); }).always(function (data) { + $('input#usgs_id').prop('disabled', false); $('input[name="impact_approach"]').prop('disabled', false); $('#getStationDataFromUsgs').prop('disabled', false); toggleRunCalcBtnState(); @@ -398,8 +402,8 @@ window.initImpactForm = function() { encode: true, }).done(function (data) { // console.log(data); - $('.impact_time_grp').css('display', 'inline-block'); - $('div.impact_time_grp').css('display', 'block'); + $('.after_get_rupture_btn').css('display', 'inline-block'); + $('div.after_get_rupture_btn').css('display', 'block'); $('#lon').val(data.lon); toggleRunCalcBtnState(); $('#lat').val(data.lat); @@ -479,7 +483,7 @@ window.initImpactForm = function() { $('#rupture-map').html('

No rupture image available

'); } } - var desc = $('#usgs_id').val() + ': '; + var desc = ''; if (data.title) { desc += data.title; } @@ -610,6 +614,8 @@ window.initImpactForm = function() { formData.append('msr', msr_selector.find(':selected').val()); } formData.append('description', $('#description').val()); + formData.append('notes', $('#notes').val()); + formData.append('make_impact_reports', $('#make_impact_reports').is(':checked')); $.ajax({ type: "POST", url: gem_oq_server_url + "/v1/calc/impact_run", diff --git a/openquake/server/templates/engine/get_outputs_impact.html b/openquake/server/templates/engine/get_outputs_impact.html index 0a0747af916e..c2042fb9fa24 100644 --- a/openquake/server/templates/engine/get_outputs_impact.html +++ b/openquake/server/templates/engine/get_outputs_impact.html @@ -16,7 +16,7 @@
-

Outputs from calculation {{ calc_id }}: {{ description }}{% if local_timestamp %} (event time: {{ local_timestamp }}){% endif %}

+

Outputs from calculation {{ calc_id }}: {% if usgs_id %}{{ usgs_id}}: {% endif %}{{ description }}{% if local_timestamp %} (event time: {{ local_timestamp }}){% endif %}

{% if time_job_after_event %}

Results computed {{ time_job_after_event }} after the event


{% endif %} {% if warnings %}
@@ -48,7 +48,11 @@

Outputs from calculation {{ calc_id }}: {{ description }}{% if local_timesta {% if aggrisk_tags %}
- Show impact table + Show impact full table +
+
+ + Show impact summary table
{% endif %} {% if mmi_tags %} @@ -70,6 +74,12 @@

Outputs from calculation {{ calc_id }}: {{ description }}{% if local_timesta

{% endif %} {% endif %} + {% for iso3 in impact_iso3_list %} +
+ + Show impact report ({{ iso3 }}) +
+ {% endfor %}
@@ -94,9 +104,25 @@

Outputs from calculation {{ calc_id }}: {{ description }}{% if local_timesta {% block templates %}