diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 9e43458..100dc14 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -1,9 +1,264 @@ # -*- coding: utf-8 -*- +from eqdist.rate import datenum_to_datetime +import Rbeast as rb; +from scipy.stats import bootstrap +from matplotlib.dates import DateFormatter, AutoDateLocator +from matplotlib.ticker import MultipleLocator +import matplotlib.pyplot as plt +import numpy as np + +global ncp_choice, tcp_max, torder_min, torder_max +ncp_choice = 'default' +tcp_max = 5 +torder_min = 0 +torder_max = 1 + +def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries, + bin_dur, unit, multiplicator, + rate_forecast, rate_unc_high, rate_unc_low, + datenum_data, mag_data): + + end_date = bin_edges[-1] + + fig, ax = plt.subplots(figsize=(14, 6)) + ax.plot(bin_edges_dt[1:], act_rate, '-o', linewidth=2.5, markersize=6.5, label='Activity rate') + + if rate_forecast is not None: + next_date = end_date + (bin_dur / multiplicator) + ax.plot(datenum_to_datetime(next_date), rate_forecast, + 'ro', label='Forecasted Rate', markersize=6.5) + ax.plot([bin_edges_dt[-1], datenum_to_datetime(next_date)], [act_rate[-1], rate_forecast], 'r-', linewidth=2.5) + ax.vlines(datenum_to_datetime(next_date), rate_unc_low, rate_unc_high, colors='r', + linewidth=2, label='Bootstrap uncertainty') + + ax.xaxis.set_major_locator(AutoDateLocator()) + ax.xaxis.set_major_formatter(DateFormatter('%d-%b-%Y')) + plt.xticks(rotation=45) + plt.title(f'Activity rate (Time Unit: {unit}, Bin Duration: {bin_dur} {unit})',fontsize=18) + # plt.title(f'Activity rate (Bin Duration: {bin_dur} {unit})',fontsize=18) + plt.xlabel('Time (Bin Center Date)', fontsize=16) + ax.set_ylabel('Activity rate per selected time period',fontsize=16) + plt.grid(True) + + if len(rt) > 0: + for i in range(len(rt)): + ax.plot(bin_edges_dt[1:][boundaries[i]:boundaries[i+1]], + [rt[i]] * (boundaries[i+1] - boundaries[i]), + linewidth=2, label=f'Rate period {i+1}') + + # ---- Magnitude scatter on right y-axis ---- + ax2 = ax.twinx() + event_dates = [datenum_to_datetime(d) for d in datenum_data] + + #-------------extract magnitude bins from data--------------------- + mags = np.array(mag_data) + min_mag = mags.min() + max_mag = mags.max() + + low_thresh = int(np.floor(min_mag)) + high_thresh = int(np.floor(max_mag)) + + thresholds = list(range(low_thresh, high_thresh + 1)) + + base_size = 15 + size_step = 35 + + bins_def = [] + for idx, t in enumerate(thresholds): + low = t + if idx < len(thresholds) - 1: + high = thresholds[idx + 1] + label = f'{low:.1f} \u2264 M < {high:.1f}' + else: + high = np.inf + label = f'M \u2265 {low:.1f}' + size = base_size + idx * size_step + bins_def.append((low, high, size, label)) + + for low, high, size, label in bins_def: + mask = (mags >= low) & (mags < high) + if np.any(mask): + sel_dates = [d for d, m in zip(event_dates, mask) if m] + sel_mags = mags[mask] + ax2.scatter(sel_dates, sel_mags, s=size, + facecolor='purple', edgecolor='black', + alpha=0.15, linewidth=1, label=label) + + ax2.set_ylabel('Magnitude', color='purple',fontsize=16) + ax2.yaxis.set_major_locator(MultipleLocator(0.5)) + ax2.yaxis.set_minor_locator(MultipleLocator(0.1)) + ax2.spines['right'].set_color('purple') + ax2.tick_params(axis='y', colors='purple') + + h1, l1 = ax.get_legend_handles_labels() + h2, l2 = ax2.get_legend_handles_labels() + handles = h1 + h2 + labels = l1 + l2 + n_legend = len(handles) + + ncols = max(1, int(np.ceil(n_legend / 5))) # ~5 entries per column + + #-------add 20% headroom above the data to make space for legend------ + ymin, ymax = ax.get_ylim() + ax.set_ylim(ymin, ymax * 1.20) + + ax.legend(handles, labels, loc='best', + ncol=ncols, borderaxespad=0,framealpha=0.7) + + ax.set_zorder(ax2.get_zorder() + 1) # put scatter plot behind the line plot + ax.patch.set_visible(False) + + fig.tight_layout() + plt.savefig("activity_rate.png", dpi=600) + plt.show() + +def bootstrap_forecast(data): + + window_data=data + + if len(window_data) >= 5: + res95 = bootstrap((window_data,), np.mean, confidence_level=0.95, + method='BCa', n_resamples=1000) + else: + res95 = bootstrap((window_data,), np.mean, confidence_level=0.95, + method='BCa', n_resamples=int(len(window_data) ** len(window_data))) + + forecast = np.mean(res95.bootstrap_distribution) + bca_conf95 = res95.confidence_interval + return forecast, bca_conf95 + +def calc_rates(act_rate, cps): + """ + Calculates mean activity rates between changepoints. + cps : sorted array of changepoint indices into act_rate + Returns rt (list of rates) and segment boundaries + """ + boundaries = [0] + list(cps.astype(int)) + [len(act_rate)] + rt = [np.mean(act_rate[boundaries[i]:boundaries[i+1]]) + for i in range(len(boundaries)-1)] + return rt, boundaries + +def apply_beast(act_rate): + """ + Applies BEAST to the smmothed rate data using different smoothing windows. + Input + act_rate : The activity rate data array to smooth and apply BEAST. + Output + out : A list of BEAST results for each smoothed rate array. + prob : A list of probabilities and change points extracted from BEAST results. + """ + + #mirror_len = int(np.ceil(0.20 * len(act_rate))) + #left_mirror = act_rate[:mirror_len][::-1] + #right_mirror = act_rate[-mirror_len:][::-1] + #act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror]) + + mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15)) + beast_result = rb.beast(act_rate, period=0, + tcp_minmax=[0, tcp_max], + torder_minmax=[torder_min, torder_max], + tseg_minlength=2, mcmc_chains=10, + mcmc_thin=mcmc_th, mcmc_seed=10) + + # User-driven ncp selection + if ncp_choice == 'median': + ncp = beast_result.trend.ncp_median + if np.isnan(ncp) or ncp == 0: + return beast_result, np.array([]) + elif ncp_choice == 'mode': + ncp = beast_result.trend.ncp_mode + if np.isnan(ncp) or ncp == 0: + return beast_result, np.array([]) + elif ncp_choice == 'pct90': + ncp = beast_result.trend.ncp_pct90 + if np.isnan(ncp) or ncp == 0: + return beast_result, np.array([]) + else: # default: median with mode and pct90 fallback + ncp = beast_result.trend.ncp_median + if np.isnan(ncp) or ncp == 0: + ncp = beast_result.trend.mode + if np.isnan(ncp) or ncp == 0: + ncp = beast_result.trend.ncp_pct90 + if np.isnan(ncp) or ncp == 0: + return beast_result, np.array([]) + + ncp = int(ncp) + # Filter NaNs first — BEAST fills unused cp slots with nan + # valid_cps = beast_result.trend.cp[~np.isnan(beast_result.trend.cp)] + cps = beast_result.trend.cp[:ncp] + + # Discard mirrored zone changepoints and correct indices + #valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate)) + #cps = cps[valid_mask] - mirror_len + + # Discard changepoints too close to the start or end (artifacts of mirroring). + # bins_after_cp / bins_before_cp set the minimum buffer bins required at each end. + #bins_before_cp = 2 + #bins_after_cp = 2 + #if len(cps) > 0: + # cps = cps[(cps >= bins_before_cp) & (cps <= len(act_rate) - bins_after_cp)] + + return beast_result, np.sort(cps) + +def bins_and_beast(dates, unit, bin_dur, multiplicator): + start_date = dates.min() + end_date = dates.max() + + valid_units = ['hours', 'days'] + if unit not in valid_units: + unit = 'days' + bin_dur = 15 + if (end_date - start_date) < 15 and unit == 'days': + unit = 'hours' + bin_dur = 12 + + bin_edges = [end_date] + while bin_edges[-1] > start_date: + bin_edges.append(bin_edges[-1] - (bin_dur / multiplicator)) + bin_edges = bin_edges[::-1] + + #-------Drop first bin or keep it if >80% of set duration------ + first_width_days = bin_edges[1] - start_date + first_width_units = first_width_days * multiplicator + + if first_width_units >= 0.8 * bin_dur: + bin_edges[0] = start_date # edge of first bin is at data start + else: + bin_edges = bin_edges[1:] # drop bin 0 (and its events) + + #------------Error if remaining bins are fewer than 2------------ + if len(bin_edges) < 2: + raise ValueError( + f"Not enough data to form at least one full bin of duration " + f"{bin_dur} {unit}(s) after dropping the partial first bin " + f"({first_width_units:.2f} {unit}(s), below the 80% threshold). " + f"Try a shorter bin_dur or check your input data range." + ) + + bin_edges_dt = [datenum_to_datetime(d) for d in bin_edges] + bin_counts, _ = np.histogram(dates, bins=bin_edges) + + act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator / bin_dur) + for i, count in enumerate(bin_counts)] + + out, cps = apply_beast(act_rate) + + if len(cps) > 0: + rt, boundaries = calc_rates(act_rate, cps) + print(f'Changepoints detected at bins: {cps}') + else: + rt = [] + boundaries = [] + print('-----------------------------------------------------') + print('No changepoints detected by BEAST (Zhao et al., 2019)') + print('-----------------------------------------------------') + + return act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, unit def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max, m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration, - forecast_select, custom_rate, forecast_len, time_unit, model, products_string, verbose): + forecast_select, custom_rate, forecast_len, time_unit, AOI_extent, model, products_string, verbose): """ Python application that reads an earthquake catalog and performs seismic hazard forecasting. Arguments: @@ -33,6 +288,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max forecasting. forecast_len: Length of the forecast for seismic hazard assessment. time_unit: Times units for the inputs Time Window Duration, Custom Activity Rate, and Forecast Length. + AOI_extent: The forecast geographical area of interest specified as a latitude and longitude range in decimal degrees + in the form [lat_min, lat_max, lon_min, lon_max]. model: Select from the following ground motion models available. Other models in the Openquake library are available but have not yet been tested. products_string: The ground motion intensity types to output. Use a space between names to select more than @@ -49,26 +306,23 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max import logging from base_logger import getDefaultLogger from timeit import default_timer as timer - from math import ceil, floor, isnan import numpy as np import dask import kalepy as kale - import utm - from skimage.transform import resize import igfash - from igfash.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf, read_csv - from igfash.window import win_CTL, win_CNE + from igfash.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf + from igfash.window import win_CNE import igfash.kde as kde from igfash.gm import compute_IMT_exceedance - from igfash.compute import get_cdf, hellinger_dist, cols_to_rows - from igfash.rate import lambda_probs, calc_bins, bootstrap_forecast_rolling - from igfash.mc import estimate_mc + from igfash.compute import get_cdf import matplotlib.pyplot as plt from matplotlib.ticker import MultipleLocator from matplotlib.contour import ContourSet import xml.etree.ElementTree as ET import json import multiprocessing as mp + import geopandas as gpd + import shapely logger = getDefaultLogger('igfash') @@ -88,9 +342,9 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max else: logger.setLevel(logging.INFO) - exclude_low_fxy = True # skip low probability areas of the map - thresh_fxy = 1e-3 # minimum fxy value (location PDF) needed to do PGA estimation (to skip low probability areas); also should scale according to number of grid points - + AOI_lat = np.array(AOI_extent[:2]) + AOI_lon = np.array(AOI_extent[2:]) + # log user selections logger.debug(f"User input files\n Catalog: {catalog_file}\n Mc: {mc_file}\n Mag_PDF: {pdf_file}\n Mag: {m_file}") logger.debug( @@ -98,7 +352,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max xy_select: {xy_select}\n grid_dim: {grid_dim}\n xy_win_method: {xy_win_method}\n rate_select: {rate_select}\n time_win_duration: {time_win_duration}\n \ forecast_select: {forecast_select}\n custom_rate: {custom_rate}\n forecast_len: {forecast_len}\n time_unit: {time_unit}\n model: {model}\n products: {products_string}\n \ verbose: {verbose}") - + logger.debug(f"Area of interest selected by user - Latitude: {AOI_lat}, Longitude: {AOI_lon}") # print key package version numbers logger.debug(f"Python version {sys.version}") logger.debug(f"Numpy version {version('numpy')}") @@ -207,45 +461,66 @@ verbose: {verbose}") time, mag, lat, lon, depth = read_mat_cat(catalog_file) - # convert to UTM - u = utm.from_latlon(lat, lon) - x = u[0] - y = u[1] - utm_zone_number = u[2] - utm_zone_letter = u[3] - logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}") + # Create a GeoDataFrame for the catalog data (initially WGS84 EPSG:4326) + catalog_gdf = gpd.GeoDataFrame( + {'depth': depth, 'time': time}, + geometry=gpd.points_from_xy(lon, lat), + crs="EPSG:4326" + ) - # define corners of grid based on global dataset - x_min = x.min() - y_min = y.min() - x_max = x.max() - y_max = y.max() + utm_crs = catalog_gdf.estimate_utm_crs() # Find the UTM EPSG code for this location + catalog_gdf_utm = catalog_gdf.to_crs(utm_crs) # Convert the catalog entirely to UTM meters + logger.debug(f"Latitude / Longitude event coordinates converted to UTM zone {utm_crs}") + + # Extract event coordinates directly from the vector geometry + x = catalog_gdf_utm.geometry.x.values + y = catalog_gdf_utm.geometry.y.values - grid_x_max = int(ceil(x_max / grid_dim) * grid_dim) - grid_x_min = int(floor(x_min / grid_dim) * grid_dim) - grid_y_max = int(ceil(y_max / grid_dim) * grid_dim) - grid_y_min = int(floor(y_min / grid_dim) * grid_dim) + # Handle Area of Interest (AOI) + if (None not in AOI_lat) and (None not in AOI_lon): + use_AOI = True + + # Create an AOI GeoDataFrame and project it to the same UTM CRS + aoi_gdf_utm = gpd.GeoDataFrame( + geometry=gpd.points_from_xy(AOI_lon, AOI_lat), + crs="EPSG:4326" + ).to_crs(utm_crs) + + # Combine dataframes to extract the collective bounding box limits + combined_gdf = gpd.GeoDataFrame(geometry=gpd.pd.concat([catalog_gdf_utm.geometry, aoi_gdf_utm.geometry])) + x_min, y_min, x_max, y_max = combined_gdf.total_bounds + else: + use_AOI = False + x_min, y_min, x_max, y_max = catalog_gdf_utm.total_bounds - grid_lat_max, grid_lon_max = utm.to_latlon(grid_x_max, grid_y_max, utm_zone_number, utm_zone_letter) - grid_lat_min, grid_lon_min = utm.to_latlon(grid_x_min, grid_y_min, utm_zone_number, utm_zone_letter) + # round up grid dimensions + grid_x_min = (x_min // grid_dim) * grid_dim + grid_x_max = ((x_max + grid_dim - 1) // grid_dim) * grid_dim + grid_y_min = (y_min // grid_dim) * grid_dim + grid_y_max = ((y_max + grid_dim - 1) // grid_dim) * grid_dim + + # expand extent until it is square + ext_w, ext_h = grid_x_max - grid_x_min, grid_y_max - grid_y_min + delta = abs(ext_w - ext_h) / 2 + + # Shift the shorter axis outward symmetrically + grid_x_min, grid_x_max = (grid_x_min - delta, grid_x_max + delta) if ext_w < ext_h else (grid_x_min, grid_x_max) + grid_y_min, grid_y_max = (grid_y_min - delta, grid_y_max + delta) if ext_h < ext_w else (grid_y_min, grid_y_max) + + # make grid points + x_range = np.arange(grid_x_min, grid_x_max + grid_dim, grid_dim) + nx = len(x_range) + y_range = np.arange(grid_y_min, grid_y_max + grid_dim, grid_dim) + ny = len(y_range) - # rectangular grid - nx = int((grid_x_max - grid_x_min) / grid_dim) + 1 - ny = int((grid_y_max - grid_y_min) / grid_dim) + 1 - - # ensure a square grid is used - if nx > ny: # enlarge y dimension to match x - ny = nx - grid_y_max = int(grid_y_min + (ny - 1) * grid_dim) - - else: # enlarge x dimension to match y - nx = ny - grid_x_max = int(grid_x_min + (nx - 1) * grid_dim) - - # new x and y range - x_range = np.linspace(grid_x_min, grid_x_max, nx) - y_range = np.linspace(grid_y_min, grid_y_max, ny) + X, Y = np.meshgrid(x_range, y_range) + cells = shapely.box(X, Y, X + grid_dim, Y + grid_dim) + grid_gdf_utm = gpd.GeoDataFrame(geometry=cells.flatten(), crs=utm_crs) + grid_gdf_latlon = grid_gdf_utm.to_crs("EPSG:4326") + logger.debug(f"Grid extent in UTM XY {grid_gdf_utm.total_bounds}") + logger.debug(f"Grid extent in lat lon {grid_gdf_latlon.total_bounds}") + t_windowed = time r_windowed = [[x, y]] @@ -268,7 +543,8 @@ verbose: {verbose}") xy_kale = output_kale[0] xy_kde = output_kde[0] - + grid_gdf_latlon['location_PDF'] = xy_kde[0].flatten() # insert location PDF as a column of the latlon GDF + # plot location PDF xy_kale_km = type(xy_kale)(xy_kale.dataset / 1000) corner = kale.corner(xy_kale_km, quantiles=[0.025, 0.16, 0.50, 0.84, 0.975], cmap='hot') @@ -320,11 +596,14 @@ verbose: {verbose}") elif rate_select: logger.info(f"Activity rate modeling selected") + + datenum_data, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True) - time, mag_dummy, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True) - - datenum_data = time # REMEMBER THE DECIMAL DENOTES DAYS - + if trim_to_mc: + indices = np.argwhere(mag_data < mc) + mag_data = np.delete(mag_data, indices) + datenum_data = np.delete(datenum_data, indices) + if time_unit == 'hours': multiplicator = 24 elif time_unit == 'days': @@ -340,32 +619,46 @@ verbose: {verbose}") logger.error(msg) raise Exception(msg) - # Selects dates in datenum format and procceeds to forecast value - start_date = datenum_data[-1] - (2 * time_win_duration / multiplicator) - dates_calc = [date for date in datenum_data if start_date <= date <= datenum_data[-1]] - forecasts, bca_conf95, rate_mean = bootstrap_forecast_rolling(dates_calc, multiplicator) + #-----------data are sorted in case they were not----------------- + sorted_pairs = sorted(zip(datenum_data, mag_data), key=lambda x: x[0]) + datenum_data, mag_data = map(list, zip(*sorted_pairs)) - # FINAL VALUES OF RATE AND ITS UNCERTAINTY IN THE 5-95 PERCENTILE - unc_bca05 = [ci.low for ci in bca_conf95]; - unc_bca95 = [ci.high for ci in bca_conf95] - rate_unc_high = multiplicator / np.array(unc_bca05); - rate_unc_low = multiplicator / np.array(unc_bca95); - rate_forecast = multiplicator / np.median(forecasts) # [per time unit] - - # Plot of forecasted activity rate with previous binned activity rate - act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit, - time_win_duration, dates_calc, - rate_forecast, rate_unc_high, rate_unc_low, - multiplicator, quiet=True, figsize=(14,9)) - - # Assign probabilities - lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges) - lambdas = np.array(lambdas, dtype='d') - lambdas_perc = np.array(lambdas_perc, dtype='d') + #-------split the data into bins and apply BEAST for changepoint detection-------------------- + act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast( + np.array(datenum_data), time_unit, time_win_duration, multiplicator) - # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1]) - logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}") - np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage", + #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected----- + if len(cps) > 0: + rate_forecast = rt[-1] + last_cp_bin = int(cps[-1]) + else: + rate_forecast = act_rate[-1] + last_cp_bin = len(act_rate) - 1 + + last_cp_datenum = bin_edges[last_cp_bin] + dates_calc = [date for date in datenum_data if last_cp_datenum <= date <= datenum_data[-1]] + interevent_times = np.diff(dates_calc) + + #------------Use BCa for uncertainty intervals----------------- + forecast, bca_conf95 = bootstrap_forecast(interevent_times) + rate_unc_high = bin_dur / (bca_conf95.low * multiplicator) + rate_unc_low = bin_dur / (bca_conf95.high * multiplicator) + + #----------------------Plot------------------------------------ + plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries, + bin_dur, time_unit, multiplicator, + rate_forecast, rate_unc_high, rate_unc_low, + datenum_data, mag_data) + + logger.info("\n----------------- Forecast Summary -----------------") + logger.info(f"Forecasted activity rate (next {bin_dur} {time_unit}(s)): {rate_forecast:.4f}") + logger.info(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]") + logger.info("------------------------------------------------------") + + lambdas = np.array([rate_forecast/bin_dur], dtype='d') + lambdas_perc = np.array([1], dtype='d') + + np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit[:-1]})", delimiter=',', fmt='%1.4f') if forecast_select: @@ -392,7 +685,7 @@ verbose: {verbose}") logger.error(msg) raise Exception(msg) - if lambdas[0] == None: + if lambdas == None: msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..." logger.error(msg) raise Exception(msg) @@ -407,52 +700,41 @@ verbose: {verbose}") m_cdf = get_cdf(m_pdf) - fxy = xy_kde[0] - logger.debug(f"Normalization check; sum of all f(x,y) values = {np.sum(fxy)}") + centroids_utm = grid_gdf_utm.geometry.centroid.values #extract the centroid of each cell + num_points = len(grid_gdf_utm) - xx, yy = np.meshgrid(x_range, y_range, indexing='ij') # grid points + distances = np.array([shapely.distance(centroids_utm[i], centroids_utm) for i in range(num_points)]) #compute distance between every grid point + grid_gdf_latlon['distance_matrix'] = [distances[i] for i in range(num_points)] #store the distance matrix in the GDF - # set every grid point to be a receiver - x_rx = xx.flatten() - y_rx = yy.flatten() - - # compute distance matrix for each receiver - distances = np.zeros(shape=(nx * ny, nx, ny)) - rx_lat = np.zeros(nx * ny) - rx_lon = np.zeros(nx * ny) - - for i in range(nx * ny): - # Compute the squared distances directly using NumPy's vectorized operations - squared_distances = (xx - x_rx[i]) ** 2 + (yy - y_rx[i]) ** 2 - distances[i] = np.sqrt(squared_distances) - - # create context object for receiver and append to list - rx_lat[i], rx_lon[i] = utm.to_latlon(x_rx[i], y_rx[i], utm_zone_number, - utm_zone_letter) # get receiver location as lat,lon - - # convert distances from m to km because openquake ground motion models take input distances in kilometres - distances = distances/1000.0 - - # compute ground motion only at grid points that have minimum probability density of thresh_fxy - if exclude_low_fxy: - indices = list(np.where(fxy.flatten() > thresh_fxy)[0]) + # Select only cells of the grid that are inside the AOI + if use_AOI: + centroids_latlon = grid_gdf_latlon.geometry.centroid + + # Mark grid cells that are within the AOI using vectorized boundary checks + grid_gdf_latlon['AOI'] = ( + (centroids_latlon.x >= AOI_lon[0]) & (centroids_latlon.x <= AOI_lon[1]) & + (centroids_latlon.y >= AOI_lat[0]) & (centroids_latlon.y <= AOI_lat[1]) + ) else: - indices = range(0, len(distances)) + grid_gdf_latlon['AOI']=True #set entire grid to be the area of interest - fr = fxy.flatten() + distances_sel = grid_gdf_latlon.loc[grid_gdf_latlon['AOI']]['distance_matrix'].to_numpy() + centroids_sel = grid_gdf_latlon.loc[grid_gdf_latlon['AOI']].centroid + + loc_pdf = grid_gdf_latlon['location_PDF'].to_numpy() # extract the previously created location PDF from the GDF + + # convert distances from m to km because openquake ground motion models take input distances in kilometres + #distances_sel = distances_sel/1000.0 # For each receiver compute estimated ground motion values - logger.info(f"Estimating ground motion intensity at {len(indices)} grid points...") - - PGA = np.zeros(shape=(nx * ny)) - - start = timer() + logger.info(f"Estimating ground motion intensity at {len(distances_sel)} grid points...") use_pp = True - + + start = timer() if use_pp: # use dask parallel computing mp.set_start_method("fork", force=True) - iter = indices + iter = range(0,len(distances_sel)) iml_grid_raw = [] # raw ground motion grids for imt in products: logger.info(f"Estimating {imt}") @@ -462,7 +744,7 @@ verbose: {verbose}") else: IMT_max = 2.0 # search interval max for acceleration (g) - imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, + imls = [dask.delayed(compute_IMT_exceedance)(centroids_sel.iloc[i].y, centroids_sel.iloc[i].x, distances_sel[i].flatten(), loc_pdf, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i, rtol=0.1, use_cython=True) for i in iter] @@ -471,9 +753,9 @@ verbose: {verbose}") else: iml_grid_raw = [] - iter = indices + iter = range(0,len(distances_sel)) + for imt in products: - if imt == "PGV": IMT_max = 200 # search interval max for velocity (cm/s) else: @@ -481,7 +763,7 @@ verbose: {verbose}") iml = [] for i in iter: - iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, + iml_i = compute_IMT_exceedance(centroids_sel.iloc[i].y, centroids_sel.iloc[i].x, distances_sel[i].flatten(), loc_pdf, p, lambdas, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True) iml.append(iml_i) @@ -491,71 +773,51 @@ verbose: {verbose}") end = timer() logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds") + + if np.isnan(iml_grid_raw).all(): msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model." logger.error(msg) raise Exception(msg) + + for j, imt in enumerate(products): #generate image overlay for each IMT product + logger.debug(f"{products[j]} values: {iml_grid_raw[j]}") + grid_gdf_latlon.loc[grid_gdf_latlon['AOI'], imt] = iml_grid_raw[j] # insert computed imt grid into GDF - # create list of one empty list for each imt - iml_grid = [[] for _ in range(len(products))] # final ground motion grids - iml_grid_prep = iml_grid.copy() # temp ground motion grids + grid_gdf_latlon_clean = grid_gdf_latlon.dropna(subset=[imt]) # remove null values from grid - if exclude_low_fxy: - for i in range(0, len(distances)): - if i in indices: - for j in range(0, len(products)): - iml_grid_prep[j].append(iml_grid_raw[j].pop(0)) - else: - list(map(lambda lst: lst.append(np.nan), - iml_grid_prep)) # use np.nan to indicate grid point excluded - else: - iml_grid_prep = iml_grid_raw + x_plot = grid_gdf_latlon_clean.geometry.centroid.x.values + y_plot = grid_gdf_latlon_clean.geometry.centroid.y.values + z_plot = grid_gdf_latlon_clean[imt].values - for j in range(0, len(products)): - vmin = min(x for x in iml_grid_prep[j] if x is not np.nan) - vmax = max(x for x in iml_grid_prep[j] if x is not np.nan) - iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype( - dtype=np.float64) # this reduces values to 8 decimal places - iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes - - # upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1 - if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1: - up_factor = 4 - iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), - mode='reflect', anti_aliasing=False) - trim_thresh = vmin - iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan - else: - iml_grid_hd = iml_grid_tmp - - iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan - - #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x)) - vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x)) - - # generate image overlay - north, south = lat.max(), lat.min() # Latitude range - east, west = lon.max(), lon.min() # Longitude range - bounds = [[south, west], [north, east]] - - map_center = [np.mean([north, south]), np.mean([east, west])] - - # Create an image from the grid - cmap_name = 'YlOrRd' - cmap = plt.get_cmap(cmap_name) - fig, ax = plt.subplots(figsize=(6, 6)) - ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear') - ax.axis('off') - - # Save the figure - fig.canvas.draw() + vmin = np.nanmin(z_plot) + vmax = np.nanmax(z_plot) + + # Generate Image Overlay + fig, ax = plt.subplots() + contour = ax.tricontourf( + x_plot, + y_plot, + z_plot, + levels=200, #linear scale + cmap="YlOrRd", + ) + ax.set_aspect('equal') # keep geographic coordinates from stretching + ax.set_axis_off() + fig.patch.set_visible(False); ax.patch.set_visible(False) overlay_filename = f"overlay_{j}.svg" - plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True) + plt.savefig(overlay_filename, pad_inches=0, bbox_inches="tight", transparent=True) plt.close(fig) + # set image map extent in geographic coordinates + north = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[3] + south = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[1] + east = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[2] + west = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[0] + # Embed geographic bounding box into the SVG map_bounds = dict(zip(("south", "west", "north", "east"), - map(float, (grid_lat_min, grid_lon_min, grid_lat_max, grid_lon_max)))) + map(float, (south, west, north, east)))) tree = ET.parse(overlay_filename) tree.getroot().set("data-map-bounds", json.dumps(map_bounds)) tree.write(overlay_filename, encoding="utf-8", xml_declaration=True) @@ -569,17 +831,17 @@ verbose: {verbose}") gradient = np.vstack((gradient, gradient)).T gradient = np.tile(gradient, (1, width)) - colorbar_title = products[j] + colorbar_title = imt if "PGA" in colorbar_title or "SA" in colorbar_title: colorbar_title = colorbar_title + " (g)" fig, ax = plt.subplots(figsize=((width + 40) / 100.0, (height + 20) / 100.0), dpi=100) # Increase fig size for labels ax.imshow(gradient, aspect='auto', cmap=cmap.reversed(), - extent=[0, 1, vmin, vmax_hd]) # Note: extent order is different for vertical + extent=[0, 1, vmin, vmax]) # Note: extent order is different for vertical ax.set_xticks([]) # Remove x-ticks for vertical colorbar num_ticks = 11 # Show more ticks - tick_positions = np.linspace(vmin, vmax_hd, num_ticks) + tick_positions = np.linspace(vmin, vmax, num_ticks) ax.set_yticks(tick_positions) ax.set_yticklabels([f"{tick:.2f}" for tick in tick_positions]) # format tick labels ax.set_title(colorbar_title, loc='right', pad=15) diff --git a/src/shf_wrapper.py b/src/shf_wrapper.py index 7675fb3..3a033cf 100644 --- a/src/shf_wrapper.py +++ b/src/shf_wrapper.py @@ -31,6 +31,9 @@ def main(argv): return False else: raise argparse.ArgumentTypeError("Boolean value expected.") + + def float_or_none(v): + return None if v.lower() == "none" else float(v) parser = argparse.ArgumentParser() @@ -55,6 +58,7 @@ def main(argv): parser.add_argument("--time_unit", type=str) parser.add_argument("--model", type=str) parser.add_argument("--products_string", type=str) + parser.add_argument("--AOI_extent", nargs=4, type=float_or_none, default=[None] * 4, required=False) parser.add_argument("--verbose", type=str2bool) args = parser.parse_args()