diff --git a/openquake/calculators/base.py b/openquake/calculators/base.py index 86b6ee1c98a..9d5fd95f843 100644 --- a/openquake/calculators/base.py +++ b/openquake/calculators/base.py @@ -37,7 +37,7 @@ import configparser import collections -from openquake.commands.plot_assets import main as plot_assets +from openquake.commands.plot_assets_webmercator import main as plot_assets from openquake.baselib import general, hdf5, config from openquake.baselib import parallel from openquake.baselib.performance import Monitor, idx_start_stop diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 233cb600902..658467044e9 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -20,6 +20,8 @@ import os import base64 import numpy +import contextily +from pyproj import Transformer from shapely.geometry import MultiPolygon from openquake.commonlib import readinput, datastore from openquake.hmtk.plotting.patch import PolygonPatch @@ -33,6 +35,23 @@ def import_plt(): return plt +def add_basemap(ax, min_x, min_y, max_x, max_y, + source=contextily.providers.CartoDB.Positron): + # NOTE: another interesting option: + # source = contextily.providers.TopPlusOpen.Grey + img, extent = contextily.bounds2img(min_x, min_y, max_x, max_y, source=source) + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) + ax.text( + 0.01, 0.01, # Position: Bottom-left corner (normalized coordinates) + source['attribution'], + transform=ax.transAxes, # Place text relative to axes + fontsize=8, color="black", alpha=0.5, + ha="left", # Horizontal alignment + va="bottom", # Vertical alignment + ) + return ax + + def adjust_limits(x_min, x_max, y_min, y_max, padding=0.5): # Make the plot display all items with some margin, looking square x_min, x_max = x_min - padding, x_max + padding @@ -46,6 +65,7 @@ def adjust_limits(x_min, x_max, y_min, y_max, padding=0.5): ylim = y_center - max_range / 2, y_center + max_range / 2 return xlim, ylim + def add_borders(ax, read_df=readinput.read_countries_df, buffer=0, alpha=0.1): plt = import_plt() polys = read_df(buffer)['geom'] @@ -113,30 +133,35 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), matplotlib.use(backend) _fig, ax = plt.subplots(figsize=figsize) ax.set_aspect('equal') - ax.grid(True) + # ax.grid(True) ax.set_xlabel('Longitude') ax.set_ylabel('Latitude') title = 'Avg GMF for %s' % imt ax.set_title(title) gmf = shakemap_array['val'][imt] - markersize = 5 - coll = ax.scatter(shakemap_array['lon'], shakemap_array['lat'], c=gmf, - cmap='jet', s=markersize) - plt.colorbar(coll) - ax = add_borders(ax, alpha=0.2) - min_x = shakemap_array['lon'].min() - max_x = shakemap_array['lon'].max() - min_y = shakemap_array['lat'].min() - max_y = shakemap_array['lat'].max() + markersize = 0.005 + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) + x_webmercator, y_webmercator = transformer.transform( + shakemap_array['lon'], shakemap_array['lat']) + min_x, min_y, max_x, max_y = (min(x_webmercator), + min(y_webmercator), + max(x_webmercator), + max(y_webmercator)) if rupture is not None: - ax, rup_min_x, rup_min_y, rup_max_x, rup_max_y = add_rupture( - ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.9, + ax, rup_min_x, rup_min_y, rup_max_x, rup_max_y = add_rupture_webmercator( + ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=1, surf_facecolor='none', surf_linestyle='--') min_x = min(min_x, rup_min_x) max_x = max(max_x, rup_max_x) min_y = min(min_y, rup_min_y) max_y = max(max_y, rup_max_y) - xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y) + xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) + min_x, max_x = xlim + min_y, max_y = ylim + add_basemap(ax, min_x, min_y, max_x, max_y) + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, + alpha=0.4) + plt.colorbar(coll, ax=ax) ax.set_xlim(*xlim) ax.set_ylim(*ylim) if with_cities: @@ -151,10 +176,9 @@ def plot_avg_gmf(ex, imt): plt = import_plt() _fig, ax = plt.subplots(figsize=(10, 10)) ax.set_aspect('equal') - ax.grid(True) + # ax.grid(True) ax.set_xlabel('Lon') ax.set_ylabel('Lat') - title = 'Avg GMF for %s' % imt assetcol = get_assetcol(ex.calc_id) if assetcol is not None: @@ -162,22 +186,22 @@ def plot_avg_gmf(ex, imt): if country_iso_codes is not None: title += ' (Countries: %s)' % country_iso_codes ax.set_title(title) - avg_gmf = ex.get('avg_gmf?imt=%s' % imt) gmf = avg_gmf[imt] markersize = 5 - coll = ax.scatter(avg_gmf['lons'], avg_gmf['lats'], c=gmf, cmap='jet', - s=markersize) - plt.colorbar(coll) - - ax = add_borders(ax) - - minx = avg_gmf['lons'].min() - maxx = avg_gmf['lons'].max() - miny = avg_gmf['lats'].min() - maxy = avg_gmf['lats'].max() - - xlim, ylim = adjust_limits(minx, maxx, miny, maxy) + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) + x_webmercator, y_webmercator = transformer.transform( + avg_gmf['lons'], avg_gmf['lats']) + min_x, min_y, max_x, max_y = (min(x_webmercator), + min(y_webmercator), + max(x_webmercator), + max(y_webmercator)) + xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) + min_x, max_x = xlim + min_y, max_y = ylim + add_basemap(ax, min_x, min_y, max_x, max_y) + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) + plt.colorbar(coll, ax=ax) ax.set_xlim(*xlim) ax.set_ylim(*ylim) return plt @@ -195,6 +219,25 @@ def add_surface(ax, surface, label, alpha=0.5, facecolor=None, linestyle='-'): return surface.get_bounding_box() +def add_surface_webmercator( + ax, surface, label, alpha=0.5, facecolor=None, linestyle='-'): + transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) + lon, lat = surface.get_surface_boundaries() + x, y = transformer.transform(lon, lat) # Transform lat/lon to Web Mercator + fill_params = { + 'alpha': alpha, + 'edgecolor': 'grey', + 'label': label + } + if facecolor is not None: + fill_params['facecolor'] = facecolor + ax.fill(x, y, **fill_params) + lon_min, lon_max, lat_max, lat_min = surface.get_bounding_box() + x_min, y_max = transformer.transform(lon_min, lat_max) # Top-left corner + x_max, y_min = transformer.transform(lon_max, lat_min) # Bottom-right corner + return x_min, x_max, y_min, y_max + + def add_rupture(ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5, surf_facecolor=None, surf_linestyle='-'): if hasattr(rup.surface, 'surfaces'): @@ -220,6 +263,30 @@ def add_rupture(ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5, return ax, min_x, min_y, max_x, max_y +def add_rupture_webmercator( + ax, rup, hypo_alpha=0.5, hypo_markersize=8, surf_alpha=0.5, + surf_facecolor=None, surf_linestyle='-'): + transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) + min_x, max_x = float('inf'), float('-inf') + min_y, max_y = float('inf'), float('-inf') + if hasattr(rup.surface, 'surfaces'): + for surf_idx, surface in enumerate(rup.surface.surfaces): + min_x_, max_x_, min_y_, max_y_ = add_surface_webmercator( + ax, surface, 'Surface %d' % surf_idx, alpha=surf_alpha, + facecolor=surf_facecolor, linestyle=surf_linestyle) + min_x, max_x = min(min_x, min_x_), max(max_x, max_x_) + min_y, max_y = min(min_y, min_y_), max(max_y, max_y_) + else: + min_x, max_x, min_y, max_y = add_surface_webmercator( + ax, rup.surface, 'Surface', alpha=surf_alpha, + facecolor=surf_facecolor, linestyle=surf_linestyle) + hypo_x, hypo_y = transformer.transform(rup.hypocenter.longitude, + rup.hypocenter.latitude) + ax.plot(hypo_x, hypo_y, marker='*', color='orange', label='Hypocenter', + alpha=hypo_alpha, linestyle='', markersize=hypo_markersize) + return ax, min_x, min_y, max_x, max_y + + def plot_rupture(rup, backend=None, figsize=(10, 10), with_cities=False, return_base64=False): # NB: matplotlib is imported inside since it is a costly import @@ -245,6 +312,33 @@ def plot_rupture(rup, backend=None, figsize=(10, 10), return plt +def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), return_base64=False): + # NB: matplotlib is imported inside since it is a costly import + plt = import_plt() + if backend is not None: + # we may need to use a non-interactive backend + import matplotlib + matplotlib.use(backend) + _fig, ax = plt.subplots(figsize=figsize) + ax.set_aspect('equal') + # ax.grid(True) + ax, min_x, min_y, max_x, max_y = add_rupture_webmercator( + ax, rup, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.3, + surf_linestyle='--') + xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) + min_x, max_x = xlim + min_y, max_y = ylim + add_basemap(ax, min_x, min_y, max_x, max_y, + source=contextily.providers.TopPlusOpen.Color) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + ax.legend() + if return_base64: + return plt_to_base64(plt) + else: + return plt + + def add_surface_3d(ax, surface, label): lon, lat, depth = surface.get_surface_boundaries_3d() lon_grid = numpy.array([[lon[0], lon[1]], [lon[3], lon[2]]]) diff --git a/openquake/commands/plot.py b/openquake/commands/plot.py index 38d5fa9a5fd..2d8fc16085e 100644 --- a/openquake/commands/plot.py +++ b/openquake/commands/plot.py @@ -29,11 +29,12 @@ from openquake.hazardlib.geo.utils import PolygonPlotter from openquake.hazardlib.contexts import Effect, get_effect_by_mag from openquake.hazardlib.calc.filters import getdefault, IntegrationDistance -from openquake.calculators.getters import get_ebrupture +from openquake.calculators.getters import get_ebrupture, get_rupture_from_dstore from openquake.calculators.extract import ( Extractor, WebExtractor, clusterize) from openquake.calculators.postproc.plots import ( - plot_avg_gmf, import_plt, add_borders, plot_rupture, plot_rupture_3d) + plot_avg_gmf, import_plt, add_borders, plot_rupture, plot_rupture_webmercator, + plot_rupture_3d) from openquake.calculators.postproc.aelo_plots import ( plot_mean_hcurves_rtgm, plot_disagg_by_src, plot_governing_mce) @@ -1064,6 +1065,16 @@ def make_figure_rupture(extractors, what): return plot_rupture(ebr.rupture) +def make_figure_rupture_webmercator(extractors, what): + """ + $ oq plot "rupture_webmercator?" + """ + [ex] = extractors + dstore = ex.dstore + rup = get_rupture_from_dstore(dstore, rup_id=0) + return plot_rupture_webmercator(rup) + + def make_figure_rupture_3d(extractors, what): """ $ oq plot "rupture_3d?" diff --git a/openquake/commands/plot_assets_webmercator.py b/openquake/commands/plot_assets_webmercator.py new file mode 100644 index 00000000000..887553d50fe --- /dev/null +++ b/openquake/commands/plot_assets_webmercator.py @@ -0,0 +1,158 @@ +# -*- coding: utf-8 -*- +# vim: tabstop=4 shiftwidth=4 softtabstop=4 +# +# Copyright (C) 2017-2023 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 numpy +import math +import shapely +import logging +from openquake.commonlib import datastore +from openquake.hazardlib.geo.utils import cross_idl # , get_bbox +from openquake.calculators.getters import get_rupture_from_dstore +from openquake.calculators.postproc.plots import ( + get_assetcol, get_country_iso_codes, add_rupture_webmercator, adjust_limits, + add_basemap) + + +def main(calc_id: int = -1, site_model=False, + save_to=None, *, show=True, assets_only=False): + """ + Plot the sites and the assets + """ + + # NB: matplotlib is imported inside since it is a costly import + import matplotlib.pyplot as p + from pyproj import Transformer + from openquake.hmtk.plotting.patch import PolygonPatch + + dstore = datastore.read(calc_id) + oq = dstore['oqparam'] + try: + region = oq.region + except KeyError: + region = None + sitecol = dstore['sitecol'] + assetcol = get_assetcol(calc_id) + _fig, ax = p.subplots(figsize=(10, 10)) + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) + if region: + region_geom = shapely.wkt.loads(region) + region_proj_geom = shapely.ops.transform(transformer.transform, region_geom) + pp = PolygonPatch(region_proj_geom, alpha=0.1) + ax.add_patch(pp) + ax.set_aspect('equal') + # ax.grid(True) + # markersize = 0.005 + # if assets_only: + # markersize_site_model = markersize_assets = 5 + # else: + # markersize_site_model = markersize_sitecol = markersize_assets = 18 + # markersize_discarded = markersize_assets + min_x, max_x, min_y, max_y = math.inf, -math.inf, math.inf, -math.inf + if site_model and 'site_model' in dstore: + sm = dstore['site_model'] + sm_lons, sm_lats = sm['lon'], sm['lat'] + if len(sm_lons) > 1 and cross_idl(*sm_lons): + sm_lons %= 360 + x_webmercator, y_webmercator = transformer.transform(sm_lons, sm_lats) + min_x, min_y, max_x, max_y = (min(x_webmercator), + min(y_webmercator), + max(x_webmercator), + max(y_webmercator)) + p.scatter(x_webmercator, y_webmercator, marker='.', color='orange', + label='site model') # , s=markersize_site_model) + # p.scatter(sitecol.complete.lons, sitecol.complete.lats, marker='.', + # color='gray', label='grid') + x_assetcol_wm, y_assetcol_wm = transformer.transform( + assetcol['lon'], assetcol['lat']) + min_x, min_y, max_x, max_y = (min(min_x, min(x_assetcol_wm)), + min(min_y, min(y_assetcol_wm)), + max(max_x, max(x_assetcol_wm)), + max(max_y, max(y_assetcol_wm))) + p.scatter(x_assetcol_wm, y_assetcol_wm, marker='.', color='green', + label='assets') # , s=markersize_assets) + if not assets_only: + x_webmercator, y_webmercator = transformer.transform( + sitecol.lons, sitecol.lats) + p.scatter(x_webmercator, y_webmercator, marker='+', color='black', + label='sites') # , s=markersize_sitecol) + if 'discarded' in dstore: + disc = numpy.unique(dstore['discarded']['lon', 'lat']) + x_webmercator, y_webmercator = transformer.transform( + disc['lon'], disc['lat']) + p.scatter(x_webmercator, y_webmercator, marker='x', color='red', + label='discarded') # , s=markersize_discarded) + if oq.rupture_xml or oq.rupture_dict: + rec = dstore['ruptures'][0] + lon, lat, _dep = rec['hypo'] + xlon, xlat = [lon], [lat] + x_webmercator, y_webmercator = transformer.transform(xlon, xlat) + dist = sitecol.get_cdist(rec) + print('rupture(%s, %s), dist=%s' % (lon, lat, dist)) + if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE': + # assuming there is only 1 rupture, so rup_id=0 + rup = get_rupture_from_dstore(dstore, rup_id=0) + ax, _minx, _miny, _maxx, _maxy = add_rupture_webmercator( + ax, rup, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.3, + surf_linestyle='--') + else: + p.scatter(x_webmercator, y_webmercator, marker='*', color='orange', + label='hypocenter', alpha=.5) + else: + xlon, xlat = [], [] + + if region: + minx, miny, maxx, maxy = region_proj_geom.bounds + else: + minx, miny, maxx, maxy = ( + min(x_assetcol_wm), min(y_assetcol_wm), + max(x_assetcol_wm), max(y_assetcol_wm)) + min_x = min(min_x, _minx, minx) + max_x = max(max_x, _maxx, maxx) + min_y = min(min_y, _miny, miny) + max_y = max(max_y, _maxy, maxy) + xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) + min_x, max_x = xlim + min_y, max_y = ylim + add_basemap(ax, min_x, min_y, max_x, max_y) + ax.set_xlim(*xlim) + ax.set_ylim(*ylim) + + country_iso_codes = get_country_iso_codes(calc_id, assetcol) + if country_iso_codes is not None: + # NOTE: use following lines to add custom items without changing title + # ax.plot([], [], ' ', label=country_iso_codes) + # ax.legend() + title = 'Countries: %s' % country_iso_codes + ax.legend(title=title) + else: + ax.legend() + + if save_to: + p.savefig(save_to, alpha=True, dpi=300) + logging.info(f'Plot saved to {save_to}') + if show: + p.show() + return p + + +main.calc_id = 'a computation id' +main.site_model = 'plot the site model too' +main.save_to = 'save the plot to this filename' +main.show = 'show the plot' +main.assets_only = 'display assets only (without sites and discarded)' diff --git a/openquake/server/views.py b/openquake/server/views.py index d46ae2bdbee..1cc37d2d72c 100644 --- a/openquake/server/views.py +++ b/openquake/server/views.py @@ -735,7 +735,7 @@ def impact_get_rupture_data(request): status=400 if 'invalid_inputs' in err else 500) if rupdic.get('shakemap_array', None) is not None: shakemap_array = rupdic['shakemap_array'] - figsize = (6.3, 6.3) # fitting in a single row in the template without resizing + figsize = (6.2, 6.2) # fitting in a single row in the template without resizing rupdic['pga_map_png'] = plot_shakemap( shakemap_array, 'PGA', backend='Agg', figsize=figsize, with_cities=False, return_base64=True, rupture=rup) diff --git a/requirements-py310-linux64.txt b/requirements-py310-linux64.txt index 6be02787ad2..703ceba51d9 100644 --- a/requirements-py310-linux64.txt +++ b/requirements-py310-linux64.txt @@ -76,3 +76,4 @@ https://wheelhouse.openquake.org/v3/py/zipp-3.17.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pluggy-1.5.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pytest-8.3.3-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/django_appconf-1.0.6-py3-none-any.whl +contextily==1.6.2 # FIXME: replace with the url to our wheelhouse diff --git a/requirements-py311-linux64.txt b/requirements-py311-linux64.txt index df7c418675f..337269c0e28 100644 --- a/requirements-py311-linux64.txt +++ b/requirements-py311-linux64.txt @@ -76,3 +76,4 @@ https://wheelhouse.openquake.org/v3/py/zipp-3.17.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pluggy-1.5.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pytest-8.3.3-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/django_appconf-1.0.6-py3-none-any.whl +contextily==1.6.2 # FIXME: replace with the url to our wheelhouse diff --git a/requirements-py312-linux64.txt b/requirements-py312-linux64.txt index 62ba153bc2f..a2a40f273a9 100644 --- a/requirements-py312-linux64.txt +++ b/requirements-py312-linux64.txt @@ -76,3 +76,4 @@ https://wheelhouse.openquake.org/v3/py/zipp-3.17.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pluggy-1.5.0-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/pytest-8.3.3-py3-none-any.whl https://wheelhouse.openquake.org/v3/py/django_appconf-1.0.6-py3-none-any.whl +contextily==1.6.2 # FIXME: replace with the url to our wheelhouse