From de205da248f0133bd6ef632a535469399342342c Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Fri, 8 Nov 2024 10:40:31 +0100 Subject: [PATCH 01/13] Use contextily and pyproj to render openstreetmap as basemap in the plots for avg_gmf and shakemap --- openquake/calculators/postproc/plots.py | 115 +++++++++++++++++------- 1 file changed, 81 insertions(+), 34 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 0c5b2b456d0a..ff603cab2227 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -94,42 +94,77 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), with_populated_places=False, return_base64=False, rupture=None): plt = import_plt() + import contextily as ctx + from pyproj import Proj, transform 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.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) - BUF_ANGLE = 1 - min_x = shakemap_array['lon'].min() - max_x = shakemap_array['lon'].max() - min_y = shakemap_array['lat'].min() - max_y = shakemap_array['lat'].max() - 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, - 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 = (min_x - BUF_ANGLE, max_x + BUF_ANGLE) - ylim = (min_y - BUF_ANGLE, max_y + BUF_ANGLE) + # markersize = 5 + markersize = 0.005 + + proj_wgs84 = Proj(init='epsg:4326') + proj_webmercator = Proj(init='epsg:3857') + + x_webmercator, y_webmercator = transform( + proj_wgs84, proj_webmercator, 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)) + + w, h = max_x - min_x, max_y - min_y + buf = 0.1 + min_x, max_x = min_x - buf * w, max_x + buf * w + min_y, max_y = min_y - buf * h, max_y + buf * h + + # BUF_ANGLE = 1 + # min_x = min_x - BUF_ANGLE + # max_x = max_x + BUF_ANGLE + # min_y = min_y - BUF_ANGLE + # max_y = max_y + BUF_ANGLE + xlim = (min_x, max_x) + ylim = (min_y, max_y) ax.set_xlim(*xlim) ax.set_ylim(*ylim) - if with_populated_places: - ax = add_populated_places(ax, xlim, ylim) + + img, extent = ctx.bounds2img( + min_x, min_y, max_x, max_y, source=ctx.providers.OpenStreetMap.Mapnik) + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) + + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, + alpha=1) + plt.colorbar(coll, ax=ax) + + # ax = add_borders(ax, alpha=0.2) + # BUF_ANGLE = 1 + # min_x = shakemap_array['lon'].min() + # max_x = shakemap_array['lon'].max() + # min_y = shakemap_array['lat'].min() + # max_y = shakemap_array['lat'].max() + # 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, + # 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 = (min_x - BUF_ANGLE, max_x + BUF_ANGLE) + # ylim = (min_y - BUF_ANGLE, max_y + BUF_ANGLE) + # ax.set_xlim(*xlim) + # ax.set_ylim(*ylim) + # if with_populated_places: + # ax = add_populated_places(ax, xlim, ylim) if return_base64: return plt_to_base64(plt) else: @@ -138,6 +173,8 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), def plot_avg_gmf(ex, imt): plt = import_plt() + import contextily as ctx + from pyproj import Proj, transform _fig, ax = plt.subplots(figsize=(10, 10)) ax.set_aspect('equal') ax.grid(True) @@ -155,19 +192,29 @@ def plot_avg_gmf(ex, imt): 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) + proj_wgs84 = Proj(init='epsg:4326') + proj_webmercator = Proj(init='epsg:3857') + + x_webmercator, y_webmercator = transform( + proj_wgs84, proj_webmercator, 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)) + w, h = max_x - min_x, max_y - min_y + min_x, max_x = min_x - 0.2 * w, max_x + 0.2 * w + min_y, max_y = min_y - 0.2 * h, max_y + 0.2 * h + img, extent = ctx.bounds2img( + min_x, min_y, max_x, max_y, source=ctx.providers.OpenStreetMap.Mapnik) + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=0.6) + + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) + plt.colorbar(coll, ax=ax) - minx = avg_gmf['lons'].min() - maxx = avg_gmf['lons'].max() - miny = avg_gmf['lats'].min() - maxy = avg_gmf['lats'].max() - w, h = maxx - minx, maxy - miny - ax.set_xlim(minx - 0.2 * w, maxx + 0.2 * w) - ax.set_ylim(miny - 0.2 * h, maxy + 0.2 * h) + ax.set_xlim(min_x, max_x) + ax.set_ylim(min_y, max_y) return plt From c868ff06d7d062135d2b30d9e21e15b9501278a6 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Mon, 18 Nov 2024 16:34:57 +0100 Subject: [PATCH 02/13] Replace deprecated method to project to webmercator also in plot_avg_gmf --- openquake/calculators/postproc/plots.py | 26 +++++++++++-------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 9dd41fff5f91..634c3b6f429e 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -161,13 +161,12 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), def plot_avg_gmf(ex, imt): plt = import_plt() import contextily as ctx - from pyproj import Proj, transform + from pyproj import Transformer _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: @@ -175,29 +174,26 @@ 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 - - proj_wgs84 = Proj(init='epsg:4326') - proj_webmercator = Proj(init='epsg:3857') - - x_webmercator, y_webmercator = transform( - proj_wgs84, proj_webmercator, avg_gmf['lons'], avg_gmf['lats']) - + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857') + # NOTE: pyproj expects latitudes first, so this is the correct order (though + # counter-intuitive) + x_webmercator, y_webmercator = transformer.transform( + avg_gmf['lats'], avg_gmf['lons']) 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 img, extent = ctx.bounds2img( min_x, min_y, max_x, max_y, source=ctx.providers.OpenStreetMap.Mapnik) - ax.imshow(img, extent=extent, interpolation='bilinear', alpha=0.6) - + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=0.3) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) plt.colorbar(coll, ax=ax) - - xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y) ax.set_xlim(*xlim) ax.set_ylim(*ylim) return plt From 0b5651ce1b106716b195b80e26b3af8bf5d1519e Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Wed, 20 Nov 2024 14:30:39 +0100 Subject: [PATCH 03/13] Use CartoDB.Positron and add attribution --- openquake/calculators/postproc/plots.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 634c3b6f429e..e3a98ea0a0ee 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -33,6 +33,17 @@ def import_plt(): return plt +def add_attribution(ax, attribution): + ax.text( + 0.01, 0.01, # Position: Bottom-left corner (normalized coordinates) + attribution, + transform=ax.transAxes, # Place text relative to axes + fontsize=8, color="black", alpha=0.5, + ha="left", # Horizontal alignment + va="bottom", # Vertical alignment + ) + + 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 @@ -144,9 +155,12 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - img, extent = ctx.bounds2img( - min_x, min_y, max_x, max_y, source=ctx.providers.OpenStreetMap.Mapnik) + source = ctx.providers.CartoDB.Positron + # NOTE: another interesting option: + # source = ctx.providers.TopPlusOpen.Grey + img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) + add_attribution(ax, source['attribution']) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, alpha=1) plt.colorbar(coll, ax=ax) @@ -189,9 +203,12 @@ def plot_avg_gmf(ex, imt): xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - img, extent = ctx.bounds2img( - min_x, min_y, max_x, max_y, source=ctx.providers.OpenStreetMap.Mapnik) + source = ctx.providers.CartoDB.Positron + # NOTE: another interesting option: + # source = ctx.providers.TopPlusOpen.Grey + img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) ax.imshow(img, extent=extent, interpolation='bilinear', alpha=0.3) + add_attribution(ax, source['attribution']) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) plt.colorbar(coll, ax=ax) ax.set_xlim(*xlim) From eefd9ae7210ac5c16fd971bc502da40b016265ff Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Wed, 20 Nov 2024 15:59:55 +0100 Subject: [PATCH 04/13] Implement add_rupture_webmercator and add_surface_webmercator --- openquake/calculators/postproc/plots.py | 75 +++++++++++++++++++------ 1 file changed, 57 insertions(+), 18 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index e3a98ea0a0ee..14edb057c320 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -20,6 +20,7 @@ import os import base64 import numpy +from pyproj import Transformer from shapely.geometry import MultiPolygon from openquake.commonlib import readinput, datastore from openquake.hmtk.plotting.patch import PolygonPatch @@ -134,24 +135,21 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), ax.set_title(title) gmf = shakemap_array['val'][imt] markersize = 0.005 - transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857') - # NOTE: pyproj expects latitudes first, so this is the correct order (though - # counter-intuitive) + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) x_webmercator, y_webmercator = transformer.transform( - shakemap_array['lat'], shakemap_array['lon']) + 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)) - # TODO: project the rupture to 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, - # 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) + if rupture is not None: + 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, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim @@ -162,7 +160,7 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) add_attribution(ax, source['attribution']) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, - alpha=1) + alpha=0.4) plt.colorbar(coll, ax=ax) ax.set_xlim(*xlim) ax.set_ylim(*ylim) @@ -191,11 +189,9 @@ def plot_avg_gmf(ex, imt): avg_gmf = ex.get('avg_gmf?imt=%s' % imt) gmf = avg_gmf[imt] markersize = 5 - transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857') - # NOTE: pyproj expects latitudes first, so this is the correct order (though - # counter-intuitive) + transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) x_webmercator, y_webmercator = transformer.transform( - avg_gmf['lats'], avg_gmf['lons']) + avg_gmf['lons'], avg_gmf['lats']) min_x, min_y, max_x, max_y = (min(x_webmercator), min(y_webmercator), max(x_webmercator), @@ -228,6 +224,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'): @@ -253,6 +268,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_populated_places=False, return_base64=False): # NB: matplotlib is imported inside since it is a costly import From 1904f4cff562540ed6a98fc4987869d7b14791da Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Wed, 20 Nov 2024 17:00:02 +0100 Subject: [PATCH 05/13] Add the possibility to plot the rupture with a webmercator basemap --- openquake/calculators/postproc/plots.py | 34 +++++++++++++++++++++++++ openquake/commands/plot.py | 13 +++++++++- 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 14edb057c320..2faa3ddbf368 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -317,6 +317,40 @@ def plot_rupture(rup, backend=None, figsize=(10, 10), return plt +def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), + with_populated_places=False, return_base64=False): + # NB: matplotlib is imported inside since it is a costly import + plt = import_plt() + import contextily as ctx + 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 + # NOTE: another interesting option: + # source = ctx.providers.CartoDB.Positron + # source = ctx.providers.TopPlusOpen.Grey + source = ctx.providers.TopPlusOpen.Color + img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) + add_attribution(ax, source['attribution']) + 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 2c6598587279..d15274a57f41 100644 --- a/openquake/commands/plot.py +++ b/openquake/commands/plot.py @@ -33,7 +33,8 @@ 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) @@ -1048,6 +1049,16 @@ def make_figure_rupture(extractors, what): return plot_rupture(rup) +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?" From d99175033e0da1d50dc1c5893c938b93655a8939 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Thu, 21 Nov 2024 15:56:09 +0100 Subject: [PATCH 06/13] Plot assets and rupture in webmercator with Positron basemap --- openquake/calculators/base.py | 4 +- openquake/commands/plot_assets_webmercator.py | 171 ++++++++++++++++++ 2 files changed, 173 insertions(+), 2 deletions(-) create mode 100644 openquake/commands/plot_assets_webmercator.py diff --git a/openquake/calculators/base.py b/openquake/calculators/base.py index 82594b265007..e78dbe5b83e1 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 @@ -658,7 +658,7 @@ def import_perils(self): # called in pre_execute logging.warning('No sites were affected by %s' % name) data[name] = peril names.append(name) - self.datastore['events'] = numpy.zeros(1, rupture.events_dt) + self.datastore['events'] = numpy.zeros(1, rupture.events_dt) create_gmf_data(self.datastore, [], names, data, N) def pre_execute(self): diff --git a/openquake/commands/plot_assets_webmercator.py b/openquake/commands/plot_assets_webmercator.py new file mode 100644 index 000000000000..f7c6896e6f18 --- /dev/null +++ b/openquake/commands/plot_assets_webmercator.py @@ -0,0 +1,171 @@ +# -*- 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) + + +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 + import contextily as ctx + 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 + source = ctx.providers.CartoDB.Positron + # NOTE: another interesting option: + # source = ctx.providers.TopPlusOpen.Grey + img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) + ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) + # FIXME: use add_attribution, moved where it does not cause a circular import + 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 + ) + 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)' From 5b2ea4764a9fcec59f1c1df06e4757ef3896fd25 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Fri, 22 Nov 2024 10:05:22 +0100 Subject: [PATCH 07/13] Extract add_basemap function adding tiles and attribution --- openquake/calculators/postproc/plots.py | 40 ++++++------------- openquake/commands/plot_assets_webmercator.py | 18 ++------- 2 files changed, 15 insertions(+), 43 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 2faa3ddbf368..9ae251404530 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -20,6 +20,7 @@ import os import base64 import numpy +import contextily as ctx from pyproj import Transformer from shapely.geometry import MultiPolygon from openquake.commonlib import readinput, datastore @@ -34,15 +35,20 @@ def import_plt(): return plt -def add_attribution(ax, attribution): +def add_basemap(ax, min_x, min_y, max_x, max_y, source=ctx.providers.CartoDB.Positron): + # NOTE: another interesting option: + # source = ctx.providers.TopPlusOpen.Grey + img, extent = ctx.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) - attribution, + 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): @@ -120,8 +126,6 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), with_populated_places=False, return_base64=False, rupture=None): plt = import_plt() - import contextily as ctx - from pyproj import Transformer if backend is not None: # we may need to use a non-interactive backend import matplotlib @@ -153,12 +157,7 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - source = ctx.providers.CartoDB.Positron - # NOTE: another interesting option: - # source = ctx.providers.TopPlusOpen.Grey - img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) - ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) - add_attribution(ax, source['attribution']) + 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) @@ -172,8 +171,6 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), def plot_avg_gmf(ex, imt): plt = import_plt() - import contextily as ctx - from pyproj import Transformer _fig, ax = plt.subplots(figsize=(10, 10)) ax.set_aspect('equal') # ax.grid(True) @@ -199,12 +196,7 @@ def plot_avg_gmf(ex, imt): xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - source = ctx.providers.CartoDB.Positron - # NOTE: another interesting option: - # source = ctx.providers.TopPlusOpen.Grey - img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) - ax.imshow(img, extent=extent, interpolation='bilinear', alpha=0.3) - add_attribution(ax, source['attribution']) + 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) @@ -317,11 +309,9 @@ def plot_rupture(rup, backend=None, figsize=(10, 10), return plt -def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), - with_populated_places=False, return_base64=False): +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() - import contextily as ctx if backend is not None: # we may need to use a non-interactive backend import matplotlib @@ -335,13 +325,7 @@ def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - # NOTE: another interesting option: - # source = ctx.providers.CartoDB.Positron - # source = ctx.providers.TopPlusOpen.Grey - source = ctx.providers.TopPlusOpen.Color - img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) - ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) - add_attribution(ax, source['attribution']) + add_basemap(ax, min_x, min_y, max_x, max_y, source=ctx.providers.TopPlusOpen.Color) ax.set_xlim(*xlim) ax.set_ylim(*ylim) ax.legend() diff --git a/openquake/commands/plot_assets_webmercator.py b/openquake/commands/plot_assets_webmercator.py index f7c6896e6f18..5408475e90e8 100644 --- a/openquake/commands/plot_assets_webmercator.py +++ b/openquake/commands/plot_assets_webmercator.py @@ -25,7 +25,8 @@ 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) + get_assetcol, get_country_iso_codes, add_rupture_webmercator, adjust_limits, + add_basemap) def main(calc_id: int = -1, site_model=False, @@ -129,20 +130,7 @@ def main(calc_id: int = -1, site_model=False, xlim, ylim = adjust_limits(min_x, max_x, min_y, max_y, padding=1E5) min_x, max_x = xlim min_y, max_y = ylim - source = ctx.providers.CartoDB.Positron - # NOTE: another interesting option: - # source = ctx.providers.TopPlusOpen.Grey - img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) - ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) - # FIXME: use add_attribution, moved where it does not cause a circular import - 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 - ) + add_basemap(ax, min_x, min_y, max_x, max_y) ax.set_xlim(*xlim) ax.set_ylim(*ylim) From 0a4d4c8790bbaec42178b6ea711c2a77b511ee10 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Fri, 22 Nov 2024 10:21:26 +0100 Subject: [PATCH 08/13] Small reduction of figure size when creating shakemap plots --- openquake/server/views.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openquake/server/views.py b/openquake/server/views.py index 0b0a3a7c6124..7fd2fd5d922f 100644 --- a/openquake/server/views.py +++ b/openquake/server/views.py @@ -771,7 +771,7 @@ def aristotle_get_rupture_data(request): del rupdic['oq_rup'] if 'shakemap_array' in rupdic: 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_populated_places=False, return_base64=True, rupture=oq_rup) From 51900b217611795fb9628bd49e715367b0213573 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Fri, 22 Nov 2024 16:48:48 +0100 Subject: [PATCH 09/13] Add contextily to linux requirements (FIXME: add url to our wheelhouse instead) --- requirements-py310-linux64.txt | 1 + requirements-py311-linux64.txt | 1 + requirements-py312-linux64.txt | 1 + 3 files changed, 3 insertions(+) diff --git a/requirements-py310-linux64.txt b/requirements-py310-linux64.txt index 6be02787ad2a..703ceba51d93 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 df7c418675fe..337269c0e28a 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 62ba153bc2f0..a2a40f273a98 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 From 914272417ecc90a1fe015e70a7ed81b11dea1055 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Fri, 22 Nov 2024 17:01:58 +0100 Subject: [PATCH 10/13] Comment out some variables that might probably be discarded --- openquake/commands/plot_assets_webmercator.py | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/openquake/commands/plot_assets_webmercator.py b/openquake/commands/plot_assets_webmercator.py index 5408475e90e8..887553d50fea 100644 --- a/openquake/commands/plot_assets_webmercator.py +++ b/openquake/commands/plot_assets_webmercator.py @@ -22,7 +22,7 @@ import shapely import logging from openquake.commonlib import datastore -from openquake.hazardlib.geo.utils import cross_idl, get_bbox +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, @@ -37,7 +37,6 @@ def main(calc_id: int = -1, site_model=False, # NB: matplotlib is imported inside since it is a costly import import matplotlib.pyplot as p - import contextily as ctx from pyproj import Transformer from openquake.hmtk.plotting.patch import PolygonPatch @@ -58,12 +57,12 @@ def main(calc_id: int = -1, site_model=False, 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 + # 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'] @@ -86,18 +85,18 @@ def main(calc_id: int = -1, site_model=False, 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) + 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) + 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) + label='discarded') # , s=markersize_discarded) if oq.rupture_xml or oq.rupture_dict: rec = dstore['ruptures'][0] lon, lat, _dep = rec['hypo'] From e97d11b925ed3fcb7ae37662113d9d65141a203e Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Mon, 2 Dec 2024 17:04:25 +0100 Subject: [PATCH 11/13] Do 'import contextily' instead of 'import contextily as ctx' --- openquake/calculators/postproc/plots.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index ba1e2a7f4da5..0e2dba5c2a4d 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -20,7 +20,7 @@ import os import base64 import numpy -import contextily as ctx +import contextily from pyproj import Transformer from shapely.geometry import MultiPolygon from openquake.commonlib import readinput, datastore @@ -35,10 +35,11 @@ def import_plt(): return plt -def add_basemap(ax, min_x, min_y, max_x, max_y, source=ctx.providers.CartoDB.Positron): +def add_basemap(ax, min_x, min_y, max_x, max_y, + source=contextily.providers.CartoDB.Positron): # NOTE: another interesting option: - # source = ctx.providers.TopPlusOpen.Grey - img, extent = ctx.bounds2img(min_x, min_y, max_x, max_y, source=source) + # 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) @@ -327,7 +328,8 @@ def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), return_base64= 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=ctx.providers.TopPlusOpen.Color) + 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() From 9dab7257be3c2f10897eb526d9a6c7e90d089d19 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Mon, 30 Jun 2025 15:33:49 +0200 Subject: [PATCH 12/13] WIP: plotting with webmercator (FIXME: fix limits of tiles) --- openquake/calculators/postproc/aelo_plots.py | 9 ++-- openquake/calculators/postproc/plots.py | 43 +++++++------------ openquake/commands/plot_assets_webmercator.py | 36 +++++----------- openquake/server/views.py | 6 +-- 4 files changed, 34 insertions(+), 60 deletions(-) diff --git a/openquake/calculators/postproc/aelo_plots.py b/openquake/calculators/postproc/aelo_plots.py index 2b3034f2b716..e28c40bcb08a 100644 --- a/openquake/calculators/postproc/aelo_plots.py +++ b/openquake/calculators/postproc/aelo_plots.py @@ -20,11 +20,10 @@ import numpy import matplotlib as mpl from scipy import interpolate -from openquake.commonlib import readinput from openquake.hazardlib.calc.mean_rates import to_rates from openquake.hazardlib.imt import from_string from openquake.calculators.extract import get_info -from openquake.calculators.postproc.plots import add_borders, adjust_limits, auto_limits +from openquake.calculators.postproc.plots import adjust_limits, auto_limits, add_basemap from PIL import Image ASCE_version = 'ASCE7-22' @@ -440,7 +439,7 @@ def plot_sites(dstore, update_dstore=False): if len(sites) == 1: markersize = 30 marker = 'x' - padding = 20 + padding = 0.005 elif len(sites) < 50: markersize = 1 marker = 'o' @@ -455,7 +454,9 @@ def plot_sites(dstore, update_dstore=False): padding = 0 plt.scatter(lons, lats, c='black', marker=marker, s=markersize) xlim, ylim = auto_limits(ax) - add_borders(ax, readinput.read_countries_df, buffer=0.) + x_min, x_max = xlim + y_min, y_max = ylim + add_basemap(ax, x_min, y_min, x_max, y_max) adjust_limits(ax, xlim, ylim, padding=padding) if update_dstore: bio = io.BytesIO() diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index e9f5841cad61..118866f0625b 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -35,11 +35,11 @@ def import_plt(): return plt -def add_basemap(ax, min_x, min_y, max_x, max_y, +def add_basemap(ax, x_min, y_min, x_max, y_max, 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) + img, extent = contextily.bounds2img(x_min, y_min, x_max, y_max, source=source) ax.imshow(img, extent=extent, interpolation='bilinear', alpha=1) ax.text( 0.01, 0.01, # Position: Bottom-left corner (normalized coordinates) @@ -155,10 +155,6 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), 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)) coll = ax.scatter(shakemap_array['lon'], shakemap_array['lat'], c=gmf, cmap='jet', s=markersize) plt.colorbar(coll) @@ -167,7 +163,9 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), ax, rupture, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=1, surf_facecolor='none', surf_linestyle='--') xlim, ylim = auto_limits(ax) - add_basemap(ax, min_x, min_y, max_x, max_y) + x_min, x_max = xlim + y_min, y_max = ylim + add_basemap(ax, x_min, y_min, x_max, y_max) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, alpha=0.4) plt.colorbar(coll, ax=ax) @@ -201,7 +199,9 @@ def plot_avg_gmf(ex, imt): x_webmercator, y_webmercator = transformer.transform( avg_gmf['lons'], avg_gmf['lats']) xlim, ylim = auto_limits(ax) - add_basemap(ax, min_x, min_y, max_x, max_y) + x_min, x_max = xlim + y_min, y_max = ylim + add_basemap(ax, x_min, y_min, x_max, y_max) coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) plt.colorbar(coll, ax=ax) adjust_limits(ax, xlim, ylim, padding=1E5) @@ -232,10 +232,6 @@ def add_surface_webmercator( 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, @@ -256,24 +252,19 @@ 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( + 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( + 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), @@ -316,17 +307,15 @@ def plot_rupture_webmercator(rup, backend=None, figsize=(10, 10), return_base64= _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( + 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, + xlim, ylim = auto_limits(ax) + x_min, x_max = xlim + y_min, y_max = ylim + add_basemap(ax, x_min, y_min, x_max, y_max, source=contextily.providers.TopPlusOpen.Color) - ax.set_xlim(*xlim) - ax.set_ylim(*ylim) - ax.legend() + adjust_limits(ax, xlim, ylim, padding=1E5) if return_base64: return plt_to_base64(plt) else: diff --git a/openquake/commands/plot_assets_webmercator.py b/openquake/commands/plot_assets_webmercator.py index 637e8c7d75b2..debd9c2e6b0b 100644 --- a/openquake/commands/plot_assets_webmercator.py +++ b/openquake/commands/plot_assets_webmercator.py @@ -18,7 +18,6 @@ import os import numpy -import math import shapely import logging from openquake.commonlib import datastore @@ -26,7 +25,7 @@ from openquake.calculators.getters import get_ebrupture from openquake.calculators.postproc.plots import ( get_assetcol, get_country_iso_codes, add_rupture_webmercator, adjust_limits, - add_basemap) + auto_limits, add_basemap) def main(calc_id: int = -1, site_model=False, @@ -63,27 +62,18 @@ def main(calc_id: int = -1, site_model=False, # 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: @@ -107,7 +97,7 @@ def main(calc_id: int = -1, site_model=False, if os.environ.get('OQ_APPLICATION_MODE') == 'ARISTOTLE': # assuming there is only 1 rupture, so rup_id=0 rup = get_ebrupture(dstore, rup_id=0) - ax, _minx, _miny, _maxx, _maxy = add_rupture_webmercator( + add_rupture_webmercator( ax, rup, hypo_alpha=0.8, hypo_markersize=8, surf_alpha=0.3, surf_linestyle='--') else: @@ -117,21 +107,15 @@ def main(calc_id: int = -1, site_model=False, xlon, xlat = [], [] if region: - minx, miny, maxx, maxy = region_proj_geom.bounds + x_min, y_min, x_max, y_max = region_proj_geom.bounds + xlim = (x_min, x_max) + ylim = (y_min, y_max) 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) + xlim, ylim = auto_limits(ax) + x_min, x_max = xlim + y_min, y_max = ylim + add_basemap(ax, x_min, y_min, x_max, y_max) + adjust_limits(ax, xlim, ylim, padding=1E5) country_iso_codes = get_country_iso_codes(calc_id, assetcol) if country_iso_codes is not None: diff --git a/openquake/server/views.py b/openquake/server/views.py index 23b7008def08..035311931135 100644 --- a/openquake/server/views.py +++ b/openquake/server/views.py @@ -58,7 +58,7 @@ from openquake.calculators.export import ( export, AGGRISK_FIELD_DESCRIPTION, EXPOSURE_FIELD_DESCRIPTION) from openquake.calculators.extract import extract as _extract -from openquake.calculators.postproc.plots import plot_shakemap, plot_rupture +from openquake.calculators.postproc.plots import plot_shakemap, plot_rupture_webmercator from openquake.engine import __version__ as oqversion from openquake.engine.export import core from openquake.engine import engine, aelo, impact @@ -793,8 +793,8 @@ def impact_get_rupture_data(request): with_cities=False, return_base64=True, rupture=rup) del rupdic['shakemap_array'] elif rup is not None: - img_base64 = plot_rupture(rup, backend='Agg', figsize=(8, 8), - return_base64=True) + img_base64 = plot_rupture_webmercator(rup, backend='Agg', figsize=(8, 8), + return_base64=True) rupdic['rupture_png'] = img_base64 return JsonResponse(rupdic, status=200) From 8387fc9ec084891a45257e3b34e306419af09c19 Mon Sep 17 00:00:00 2001 From: Paolo Tormene Date: Mon, 30 Jun 2025 17:47:17 +0200 Subject: [PATCH 13/13] Fix plot_shakemap and plot_avg_gmf --- openquake/calculators/postproc/plots.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/openquake/calculators/postproc/plots.py b/openquake/calculators/postproc/plots.py index 118866f0625b..a5bd0a5b6bf6 100644 --- a/openquake/calculators/postproc/plots.py +++ b/openquake/calculators/postproc/plots.py @@ -155,8 +155,8 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) x_webmercator, y_webmercator = transformer.transform( shakemap_array['lon'], shakemap_array['lat']) - coll = ax.scatter(shakemap_array['lon'], shakemap_array['lat'], c=gmf, - cmap='jet', s=markersize) + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, + alpha=0.4) plt.colorbar(coll) if rupture is not None: add_rupture_webmercator( @@ -166,10 +166,7 @@ def plot_shakemap(shakemap_array, imt, backend=None, figsize=(10, 10), x_min, x_max = xlim y_min, y_max = ylim add_basemap(ax, x_min, y_min, x_max, y_max) - coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize, - alpha=0.4) - plt.colorbar(coll, ax=ax) - adjust_limits(ax, xlim, ylim, padding=1E5) + adjust_limits(ax, xlim, ylim, padding=0.005) if with_cities: add_cities(ax, xlim, ylim) if return_base64: @@ -198,12 +195,12 @@ def plot_avg_gmf(ex, imt): transformer = Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True) x_webmercator, y_webmercator = transformer.transform( avg_gmf['lons'], avg_gmf['lats']) + coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) + plt.colorbar(coll, ax=ax) xlim, ylim = auto_limits(ax) x_min, x_max = xlim y_min, y_max = ylim add_basemap(ax, x_min, y_min, x_max, y_max) - coll = ax.scatter(x_webmercator, y_webmercator, c=gmf, cmap='jet', s=markersize) - plt.colorbar(coll, ax=ax) adjust_limits(ax, xlim, ylim, padding=1E5) return plt