Merge branch 'master' of github.com:gem/oq-engine

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

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,

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': '<in-memory>'}
     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

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)

openquake/hazardlib/shakemap/parsers.py

@@ -25,13 +25,15 @@
 from urllib.error import URLError
 from collections import defaultdict
 
+import tempfile
 import io
 import os
 import pathlib
 import logging
 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'

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());

openquake/server/templates/engine/index.html

@@ -159,6 +159,11 @@ <h3>
                 <label for"ses_seed">{{ aristotle_form_labels.ses_seed }}</label>
                 <input class="aristotle-input" type="text" id="ses_seed" name="ses_seed" placeholder="{{ aristotle_form_placeholders.ses_seed }}" value="42"/>
               </div>
+              <div class="aristotle_form_row">
+                <label for"station_data_file_from_usgs">{{ aristotle_form_labels.station_data_file_from_usgs }}</label>
+                <input class="aristotle-input" type="hidden" id="station_data_file_from_usgs" name="station_data_file_from_usgs" placeholder="{{ aristotle_form_placeholders.station_data_file_from_usgs }}" disabled />
+                <input class="aristotle-input" type="text" id="station_data_file_from_usgs_loaded" name="station_data_file_from_usgs_loaded" disabled />
+              </div>
               <div class="aristotle_form_row">
                 <label for"station_data_file">{{ aristotle_form_labels.station_data_file }}</label>
                 <input class="aristotle-input" type="file" id="station_data_file_input" name="station_data_file" placeholder="{{ aristotle_form_placeholders.station_data_file }}" >

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()