diff --git a/openquake/calculators/views.py b/openquake/calculators/views.py index 0a94040b1d1..1fb4e53dae5 100644 --- a/openquake/calculators/views.py +++ b/openquake/calculators/views.py @@ -45,7 +45,8 @@ from openquake.risklib.scientific import ( losses_by_period, return_periods, LOSSID, LOSSTYPE) from openquake.baselib.writers import build_header, scientificformat -from openquake.calculators.getters import get_ebrupture, MapGetter, get_pmaps_gb +from openquake.calculators.getters import ( + get_ebrupture, MapGetter, get_pmaps_gb) from openquake.calculators.extract import extract TWO24 = 2**24 diff --git a/openquake/commands/mosaic.py b/openquake/commands/mosaic.py index 645a04d23eb..3ac17112a06 100644 --- a/openquake/commands/mosaic.py +++ b/openquake/commands/mosaic.py @@ -314,7 +314,9 @@ def aristotle(exposure_hdf5=None, *, trt=trt, truncation_level=truncation_level, number_of_ground_motion_fields=number_of_ground_motion_fields, asset_hazard_distance=asset_hazard_distance, - ses_seed=ses_seed, exposure_hdf5=exposure_hdf5) + ses_seed=ses_seed, exposure_hdf5=exposure_hdf5, + station_data_file=stations, + maximum_distance_stations=maximum_distance_stations) else: # assume .xml main_cmd('WithRuptureFile', rupfname, None, callback, maximum_distance=maximum_distance, diff --git a/openquake/engine/aristotle.py b/openquake/engine/aristotle.py index 806681e0556..9b50518108c 100644 --- a/openquake/engine/aristotle.py +++ b/openquake/engine/aristotle.py @@ -25,7 +25,8 @@ import numpy from openquake.baselib import config, hdf5, sap from openquake.hazardlib import geo, nrml, sourceconverter -from openquake.hazardlib.shakemap.parsers import download_rupture_dict +from openquake.hazardlib.shakemap.parsers import ( + download_rupture_dict, download_station_data_file) from openquake.commonlib import readinput from openquake.engine import engine @@ -115,6 +116,9 @@ def get_aristotle_allparams(rupture_dict, time_event, inputs = {'exposure': [exposure_hdf5], 'job_ini': ''} rupdic = get_rupture_dict(rupture_dict, ignore_shakemap) + if station_data_file is None: + # NOTE: giving precedence to the station_data_file uploaded via form + station_data_file = download_station_data_file(rupture_dict['usgs_id']) rupture_file = rupdic.pop('rupture_file') if rupture_file: inputs['rupture_model'] = rupture_file diff --git a/openquake/hazardlib/calc/conditioned_gmfs.py b/openquake/hazardlib/calc/conditioned_gmfs.py index ecb7fe38b87..6cbcd86df8b 100644 --- a/openquake/hazardlib/calc/conditioned_gmfs.py +++ b/openquake/hazardlib/calc/conditioned_gmfs.py @@ -278,8 +278,8 @@ def _create_result(target_imt, observed_imts, station_data_filtered): if num_null_values: raise ValueError( f"The station data contains {num_null_values}" - f"null values for {target_imt.string}." - "Please fill or discard these rows.") + f" null values for {target_imt.string}." + " Please fill or discard these rows.") t = TempResult(conditioning_imts=conditioning_imts, bracketed_imts=bracketed_imts, native_data_available=native_data_available) diff --git a/openquake/hazardlib/shakemap/parsers.py b/openquake/hazardlib/shakemap/parsers.py index 461c0409e32..a5c75932e81 100644 --- a/openquake/hazardlib/shakemap/parsers.py +++ b/openquake/hazardlib/shakemap/parsers.py @@ -25,6 +25,7 @@ from urllib.error import URLError from collections import defaultdict +import tempfile import io import os import pathlib @@ -32,6 +33,7 @@ import json import zipfile import pytz +import pandas as pd from datetime import datetime from shapely.geometry import Polygon import numpy @@ -234,6 +236,172 @@ def local_time_to_time_event(local_time): return 'transit' +def read_usgs_stations_json(stations_json_str): + sj = json.loads(stations_json_str) + if 'features' not in sj or not sj['features']: + raise LookupError('Station data is not available yet.') + stations = pd.json_normalize(sj, 'features') + stations['eventid'] = sj['metadata']['eventid'] + # Rename columns + stations.columns = [ + col.replace('properties.', '') for col in stations.columns] + # Extract lon and lat + stations[['lon', 'lat']] = pd.DataFrame( + stations['geometry.coordinates'].to_list()) + # Get values for available IMTs (PGA and SA) + # ========================================== + # The "channels/amplitudes" dictionary contains the values recorded at + # the seismic stations. The values could report the 3 components, in such + # cases, take the componet with maximum PGA (and in absence of PGA, the + # first IM reported). + channels = pd.DataFrame(stations.channels.to_list()) + vals = pd.Series([], dtype='object') + for row, rec_station in channels.iterrows(): + rec_station.dropna(inplace=True) + # Iterate over different columns. Each colum can be a component + data = [] + pgas = [] + for _, chan in rec_station.items(): + if chan["name"].endswith("Z") or chan["name"].endswith("U"): + continue + # logging.info(chan["name"]) + df = pd.DataFrame(chan["amplitudes"]) + if 'pga' in df.name.unique(): + pga = df.loc[df['name'] == 'pga', 'value'].values[0] + else: + pga = df['value'][0] + if pga is None or pga == "null": + continue + if isinstance(pga, str): + pga = float(pga) + pgas.append(pga) + data.append(chan["amplitudes"]) + # get values for maximum component + if pgas: + max_componet = pgas.index(max(pgas)) + vals[row] = data[max_componet] + else: + vals[row] = None + # The "pgm_from_mmi" dictionary contains the values estimated from MMI. + # Combine both dictionaries to extract the values. + # They are generally mutually exclusive (if mixed, the priority is given + # to the station recorded data). + try: + # Some events might not have macroseismic data, then skip them + vals = vals.combine_first(stations['pgm_from_mmi']).apply(pd.Series) + except Exception: + vals = vals.apply(pd.Series) + # Arrange columns since the data can include mixed positions for the IMTs + values = pd.DataFrame() + for col in vals.columns: + df = vals[col].apply(pd.Series) + df.set_index(['name'], append=True, inplace=True) + df.drop(columns=['flag', 'units'], inplace=True) + if 0 in df.columns: + df.drop(columns=[0], inplace=True) + df = df.unstack('name') + df.dropna(axis=1, how='all', inplace=True) + df.columns = [col[1]+'_'+col[0] for col in df.columns.values] + for col in df.columns: + if col in values: + # Colum already exist. Combine values in unique column + values[col] = values[col].combine_first(df[col]) + else: + values = pd.concat([values, df[col]], axis=1) + values.sort_index(axis=1, inplace=True) + # Add recording to main DataFrame + stations = pd.concat([stations, values], axis=1) + return stations + + +def usgs_to_ecd_format(stations, exclude_imts=()): + ''' + Adjust USGS format to match the ECD (Earthquake Consequence Database) + format + ''' + # Adjust column names to match format + stations.columns = stations.columns.str.upper() + stations.rename(columns={ + 'CODE': 'STATION_ID', + 'NAME': 'STATION_NAME', + 'LON': 'LONGITUDE', + 'LAT': 'LATITUDE', + 'INTENSITY': 'MMI_VALUE', + 'INTENSITY_STDDEV': 'MMI_STDDEV', + }, inplace=True) + # Identify columns for IMTs: + imts = [] + for col in stations.columns: + if 'DISTANCE_STDDEV' == col: + continue + if '_VALUE' in col or '_LN_SIGMA' in col or '_STDDEV' in col: + for imt in exclude_imts: + if imt in col: + break + else: + imts.append(col) + # Identify relevant columns + cols = ['STATION_ID', 'STATION_NAME', 'LONGITUDE', 'LATITUDE', + 'STATION_TYPE', 'VS30'] + imts + df = stations[cols].copy() + # Add missing columns + df.loc[:, 'VS30_TYPE'] = 'inferred' + df.loc[:, 'REFERENCES'] = 'Stations_USGS' + # Adjust PGA and SA untis to [g]. USGS uses [% g] + adj_cols = [imt for imt in imts + if '_VALUE' in imt and + 'PGV' not in imt and + 'MMI' not in imt] + df.loc[:, adj_cols] = round(df.loc[:, adj_cols]. + apply(pd.to_numeric, errors='coerce') / 100, 6) + df_seismic = df[df['STATION_TYPE'] == 'seismic'] + return df_seismic + + +def download_station_data_file(usgs_id): + """ + Download station data from the USGS site given a ShakeMap ID. + + :param usgs_id: ShakeMap ID + :returns: the path of a csv file with station data converted to a format + compatible with OQ + """ + # NOTE: downloading twice from USGS, but with a clearer workflow + url = SHAKEMAP_URL.format(usgs_id) + logging.info('Downloading %s' % url) + js = json.loads(urlopen(url).read()) + products = js['properties']['products'] + station_data_file = None + try: + shakemap = products['shakemap'] + except KeyError: + logging.warning('No shakemap was found') + return None + for shakemap in reversed(shakemap): + contents = shakemap['contents'] + if 'download/stationlist.json' in contents: + stationlist_url = contents.get('download/stationlist.json')['url'] + logging.info('Downloading stationlist.json') + stations_json_str = urlopen(stationlist_url).read() + try: + stations = read_usgs_stations_json(stations_json_str) + except LookupError as exc: + logging.warning(str(exc)) + else: + df = usgs_to_ecd_format(stations, exclude_imts=('SA(3.0)',)) + if len(df) < 1: + logging.warning('No seismic stations found') + else: + with tempfile.NamedTemporaryFile( + delete=False, mode='w+', newline='', + suffix='.csv') as temp_file: + station_data_file = temp_file.name + df.to_csv(station_data_file, encoding='utf8', + index=False) + logging.info(f'Wrote stations to {station_data_file}') + return station_data_file + + def download_rupture_dict(id, ignore_shakemap=False): """ Download a rupture from the USGS site given a ShakeMap ID. @@ -243,7 +411,7 @@ def download_rupture_dict(id, ignore_shakemap=False): :returns: a dictionary with keys lon, lat, dep, mag, rake """ url = SHAKEMAP_URL.format(id) - print('Downloading %s' % url) + logging.info('Downloading %s' % url) js = json.loads(urlopen(url).read()) mag = js['properties']['mag'] products = js['properties']['products'] @@ -255,7 +423,8 @@ def download_rupture_dict(id, ignore_shakemap=False): try: products['finite-fault'] except KeyError: - raise MissingLink('There is no finite-fault info for %s' % id) + raise MissingLink( + 'There is no shakemap nor finite-fault info for %s' % id) else: shakemap = [] for shakemap in reversed(shakemap): @@ -267,7 +436,7 @@ def download_rupture_dict(id, ignore_shakemap=False): ff = products['finite-fault'] except KeyError: raise MissingLink('There is no finite-fault info for %s' % id) - print('Getting finite-fault properties') + logging.info('Getting finite-fault properties') if isinstance(ff, list): if len(ff) > 1: logging.warning(f'The finite-fault list contains {len(ff)}' @@ -285,7 +454,7 @@ def download_rupture_dict(id, ignore_shakemap=False): 'is_point_rup': False, 'usgs_id': id, 'rupture_file': None} return rupdic url = contents.get('download/rupture.json')['url'] - print('Downloading rupture.json') + logging.info('Downloading rupture.json') rup_data = json.loads(urlopen(url).read()) feats = rup_data['features'] is_point_rup = len(feats) == 1 and feats[0]['geometry']['type'] == 'Point' diff --git a/openquake/server/static/js/engine.js b/openquake/server/static/js/engine.js index 955fc0b6f11..f68917b8b93 100644 --- a/openquake/server/static/js/engine.js +++ b/openquake/server/static/js/engine.js @@ -556,6 +556,12 @@ $('#local_timestamp').val(data.local_timestamp); $('#time_event').val(data.time_event); $('#is_point_rup').val(data.is_point_rup); + // NOTE: due to security restrictions in web browsers, it is not possible to programmatically + // set a specific file in an HTML file input element using JavaScript or jQuery, + // therefore we can not pre-populate the station_data_file_input with the station_data_file + // obtained converting the USGS stationlist.json, and we use a separate field referencing it + $('#station_data_file_from_usgs').val(data.station_data_file_from_usgs); + $('#station_data_file_from_usgs_loaded').val(data.station_data_file_from_usgs ? 'Loaded' : 'N.A.'); if ($('#rupture_file_input')[0].files.length == 1) { $('#dip').prop('disabled', true); $('#strike').prop('disabled', true); @@ -633,6 +639,7 @@ $('#number_of_ground_motion_fields').val()); formData.append('asset_hazard_distance', $('#asset_hazard_distance').val()); formData.append('ses_seed', $('#ses_seed').val()); + formData.append('station_data_file_from_usgs', $('#station_data_file_from_usgs').val()); formData.append('local_timestamp', $("#local_timestamp").val()); formData.append('station_data_file', $('#station_data_file_input')[0].files[0]); formData.append('maximum_distance_stations', $("#maximum_distance_stations").val()); diff --git a/openquake/server/templates/engine/index.html b/openquake/server/templates/engine/index.html index f75fb1dcabd..2a4f57e7385 100644 --- a/openquake/server/templates/engine/index.html +++ b/openquake/server/templates/engine/index.html @@ -159,6 +159,11 @@

+
+ + + +
diff --git a/openquake/server/views.py b/openquake/server/views.py index 550705a92ce..176af63c013 100644 --- a/openquake/server/views.py +++ b/openquake/server/views.py @@ -57,6 +57,7 @@ from openquake.engine.aelo import ( get_params_from, PRELIMINARY_MODELS, PRELIMINARY_MODEL_WARNING) from openquake.engine.export.core import DataStoreExportError +from openquake.hazardlib.shakemap.parsers import download_station_data_file from openquake.engine.aristotle import ( get_trts_around, get_aristotle_allparams, get_rupture_dict) from openquake.server import utils @@ -129,6 +130,7 @@ 'number_of_ground_motion_fields': 'Number of ground motion fields', 'asset_hazard_distance': 'Asset hazard distance (km)', 'ses_seed': 'Random seed', + 'station_data_file_from_usgs': 'Station data from USGS', 'station_data_file': 'Station data CSV', 'maximum_distance_stations': 'Maximum distance of stations (km)', } @@ -151,6 +153,7 @@ 'number_of_ground_motion_fields': 'float ≥ 1', 'asset_hazard_distance': 'float ≥ 0', 'ses_seed': 'int ≥ 0', + 'station_data_file_from_usgs': '', 'station_data_file': 'Station data CSV', 'maximum_distance_stations': 'float ≥ 0', } @@ -712,7 +715,7 @@ def aristotle_get_rupture_data(request): res = aristotle_validate(request) if isinstance(res, HttpResponse): # error return res - [rupdic] = res + rupdic, station_data_file = res try: trts, mosaic_model = get_trts_around( rupdic, os.path.join(config.directory.mosaic_dir, 'exposure.hdf5')) @@ -732,6 +735,7 @@ def aristotle_get_rupture_data(request): content=json.dumps(response_data), content_type=JSON, status=400) rupdic['trts'] = trts rupdic['mosaic_model'] = mosaic_model + rupdic['station_data_file_from_usgs'] = station_data_file response_data = rupdic return HttpResponse(content=json.dumps(response_data), content_type=JSON, status=200) @@ -823,11 +827,6 @@ def aristotle_validate(request): if 'is_point_rup' in request.POST: dic['is_point_rup'] = request.POST['is_point_rup'] == 'true' - # FIXME: validate station_data_file - if 'lon' in request.POST: - # NOTE: if the request contains all form parameters - params['station_data_file'] = station_data_path - if validation_errs: err_msg = 'Invalid input value' err_msg += 's\n' if len(validation_errs) > 1 else '\n' @@ -851,6 +850,15 @@ def aristotle_validate(request): logging.error('', exc_info=True) return HttpResponse(content=json.dumps(response_data), content_type=JSON, status=500) + if station_data_path is not None: + # giving precedence to the user-uploaded station data file + params['station_data_file'] = station_data_path + elif request.POST.get('station_data_file_from_usgs'): + params['station_data_file'] = request.POST.get( + 'station_data_file_from_usgs') + else: + station_data_file = download_station_data_file(dic['usgs_id']) + params['station_data_file'] = station_data_file return rupdic, *params.values()