diff --git a/LICENSE b/LICENSE index 839fee2aa9d3796247f073aa9f36b65800fab5c6..3e802333558816c73a3b9a1251d2cf0a251f5335 100644 --- a/LICENSE +++ b/LICENSE @@ -629,8 +629,8 @@ to attach them to the start of each source file to most effectively state the exclusion of warranty; and each file should have at least the "copyright" line and a pointer to where the full notice is found. - ramsis-sfm - Copyright (C) 2021 Indu + hermes-sfm + Copyright (C) 2021 ETH Zürich This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published diff --git a/em1/ramsis.py b/em1/hermes.py similarity index 93% rename from em1/ramsis.py rename to em1/hermes.py index 0a98712b1e07f6532a74d928628e0f773dafb8f7..238e0dc79ad5b905b998c25fe71487db2b3e8b22 100644 --- a/em1/ramsis.py +++ b/em1/hermes.py @@ -1,4 +1,4 @@ -from ramsis_model import ModelInput, validate_entrypoint +from hermes_model import ModelInput, validate_entrypoint from em1.model import EM1 diff --git a/em1/model.py b/em1/model.py index b757f0e5d009e19f43569fdc1fe6a0bafc75dd3f..e9eb40d9e465eb1d0c7c2809abe68c6296cc6d9c 100644 --- a/em1/model.py +++ b/em1/model.py @@ -3,12 +3,11 @@ from datetime import timedelta import numpy as np import scipy.optimize as opt +from hermes_model import ModelInput from hydws.parser import BoreholeHydraulics -from ramsis_model import ModelInput from scipy import stats from scipy.special import lambertw -from seismostats.seismicity.catalog import Catalog -from seismostats.seismicity.rategrid import ForecastGRRateGrid +from seismostats import Catalog, ForecastGRRateGrid from shapely import wkt from em1.likelihood import log_lhood_comp, log_lhood_comp_FF, posterior, prior @@ -60,8 +59,10 @@ class EM1(): # get dataframes from HYDSJONParser self.hydraulic_df = self.hydraulics[section_id].hydraulics - self.hydraulic_df = self.hydraulic_df[~self.hydraulic_df['topflow'].isna()] - self.hydraulic_df = self.hydraulic_df[~self.hydraulic_df['toppressure'].isna()] + self.hydraulic_df = self.hydraulic_df[ + ~self.hydraulic_df['topflow'].isna()] + self.hydraulic_df = self.hydraulic_df[ + ~self.hydraulic_df['toppressure'].isna()] self.scenario_df = self.scenario[section_id].hydraulics @@ -86,7 +87,8 @@ class EM1(): self.seismicity['associatedphasecount'].fillna(0) self.seismicity = \ - self.seismicity[self.seismicity['associatedphasecount'] >= self.parameters['n_phases']] + self.seismicity[self.seismicity['associatedphasecount'] + >= self.parameters['n_phases']] self.seismicity = self.seismicity[~self.seismicity['magnitude'].isna()] @@ -109,8 +111,10 @@ class EM1(): sc_all_identical = self.scenario_df['topflow'].nunique() == 1 if sc_all_negative and sc_all_identical: - self.scenario_df["topflow"]=(-self.scenario_df["topflow"].values[0]* - np.ones(len(self.scenario_df))*self.hydraulic_df["topflow"].values[-1]) + self.scenario_df["topflow"] = \ + (-self.scenario_df["topflow"].values[0] + * np.ones(len(self.scenario_df)) + * self.hydraulic_df["topflow"].values[-1]) def run(self): """ diff --git a/model_plots.py b/model_plots.py index fd3e081efd68ef8ea41667d79fbe8c907895fc9d..58eb24acf06eb462240cbb9306457007d4c6770a 100644 --- a/model_plots.py +++ b/model_plots.py @@ -2,18 +2,17 @@ import logging import os from datetime import datetime -import pandas as pd import numpy as np +import pandas as pd from hydws.parser import BoreholeHydraulics -from seismostats.seismicity.catalog import Catalog +from seismostats import Catalog from seismostats.utils.coordinates import CoordinateTransformer from shapely import Polygon from shapely.wkt import dumps -from em1.ramsis import ModelInput, run_em1, run_em1_validated # noqa +from em1.hermes import ModelInput, run_em1, run_em1_validated # noqa from em1.model import EM1 -from em1.utils import plotting_data,plotting_results - +from em1.utils import plotting_data, plotting_results # Set up logger date_strftime_format = "%d-%m-%y %H:%M:%S" @@ -53,8 +52,6 @@ def main(): parser[section_id].hydraulics = hydraulics_plan scenario_hydjson = parser.to_json() - - ############################################################ # SEISMICITY @@ -135,20 +132,22 @@ def main(): results = model.predict() - #results = model.parse_results(results) + # results = model.parse_results(results) catalog_df_full = pd.read_hdf(Data_Path + 'seismic_Data_full_located.h5') - catalog_df_full.rename(columns={"X": "Xs", "Y": "Ys", "Depth": "Zs", "MomMag": "magnitude", "Time(UTC)": "time"}, inplace=True) + catalog_df_full.rename(columns={"X": "Xs", "Y": "Ys", "Depth": "Zs", + "MomMag": "magnitude", "Time(UTC)": "time"}, + inplace=True) catalog_df_full = catalog_df_full[catalog_df_full["magnitude"] > -5] catalog_df_full.set_index("time", inplace=True) - catalog_df_full = catalog_df_full[catalog_df_full.index >= hydraulics.index.values[0]] + catalog_df_full = catalog_df_full[catalog_df_full.index + >= hydraulics.index.values[0]] catalog_df_full = catalog_df_full[catalog_df_full.index <= forecast_end] catalog_df_full = catalog_df_full.reset_index() - plotting_results(model,results,catalog_df_full) + plotting_results(model, results, catalog_df_full) return model if __name__ == "__main__": result = main() - diff --git a/model_plots_forge.py b/model_plots_forge.py index 55c735538ac51566861feb3b47334bc4e0200da7..c851c6a0456443927d3d09b94787f7295d2d074c 100644 --- a/model_plots_forge.py +++ b/model_plots_forge.py @@ -2,18 +2,17 @@ import logging import os from datetime import datetime -import pandas as pd import numpy as np +import pandas as pd from hydws.parser import BoreholeHydraulics -from seismostats.seismicity.catalog import Catalog +from seismostats import Catalog from seismostats.utils.coordinates import CoordinateTransformer from shapely import Polygon from shapely.wkt import dumps -from em1.ramsis import ModelInput, run_em1, run_em1_validated # noqa +from em1.hermes import ModelInput, run_em1, run_em1_validated # noqa from em1.model import EM1 -from em1.utils import plotting_data,plotting_results - +from em1.utils import plotting_data, plotting_results # Set up logger date_strftime_format = "%d-%m-%y %H:%M:%S" @@ -28,33 +27,43 @@ def main(): # HYDRAULICS dirpath = os.path.dirname(os.path.abspath(__file__)) Data_Path = os.path.join(dirpath, "data/") - #hydraulics = pd.read_hdf(Data_Path + 'TrainingData.h5') + # hydraulics = pd.read_hdf(Data_Path + 'TrainingData.h5') hydraulics = pd.read_hdf(Data_Path + 'Hydr_FORGE22.h5') - hydraulics_plan = hydraulics[hydraulics.index.values >= np.array(datetime(2022,4,21,15,0)).astype('datetime64')] - hydraulics_plan = hydraulics_plan[hydraulics_plan.index.values < np.array(datetime(2022,4,21,18)).astype('datetime64')] - + hydraulics_plan = hydraulics[hydraulics.index.values >= np.array( + datetime(2022, 4, 21, 15, 0)).astype('datetime64')] + hydraulics_plan = hydraulics_plan[hydraulics_plan.index.values < np.array( + datetime(2022, 4, 21, 18)).astype('datetime64')] # define start and end of forecast - forecast_start =datetime(2022,4,21,15, 0)#datetime.utcfromtimestamp(hydraulics_plan.index.values[0].tolist() / 1e9) - forecast_end = datetime(2022,4,21,18,0)#datetime.utcfromtimestamp(hydraulics_plan.index.values[-1].tolist() / 1e9) - - hydraulics = hydraulics[hydraulics.index.values >= np.array(datetime(2022,4,21,13,00,23)).astype('datetime64')] - hydraulics = hydraulics[hydraulics.index.values < np.array(datetime(2022,4,21,15,0)).astype('datetime64')] - - hydraulics_plan["topflow"]=-1*np.ones((len(hydraulics_plan),1)) - hydraulics_plan.index.values[-1]=np.array(datetime(2022,4,21,18)).astype('datetime64') + # datetime.utcfromtimestamp( + # hydraulics_plan.index.values[0].tolist() / 1e9) + forecast_start = datetime(2022, 4, 21, 15, 0) + # datetime.utcfromtimestamp( + # hydraulics_plan.index.values[-1].tolist() / 1e9) + forecast_end = datetime(2022, 4, 21, 18, 0) + + hydraulics = hydraulics[hydraulics.index.values >= np.array( + datetime(2022, 4, 21, 13, 00, 23)).astype('datetime64')] + hydraulics = hydraulics[hydraulics.index.values < np.array( + datetime(2022, 4, 21, 15, 0)).astype('datetime64')] + + hydraulics_plan["topflow"] = -1 * np.ones((len(hydraulics_plan), 1)) + hydraulics_plan.index.values[-1] = np.array( + datetime(2022, 4, 21, 18)).astype('datetime64') hydraulics_plan.topflow.values[-1] = hydraulics_plan.topflow.values[-2] - hydraulics_plan=hydraulics_plan[~hydraulics_plan.index.duplicated(keep='first')] - - #hydraulics_plan = pd.read_hdf(Data_Path + 'Scenario.h5') + hydraulics_plan = hydraulics_plan[~hydraulics_plan.index.duplicated( + keep='first')] + # hydraulics_plan = pd.read_hdf(Data_Path + 'Scenario.h5') - observation_start = datetime.utcfromtimestamp(hydraulics.index.values[0].tolist() / 1e9) - observation_end = datetime.utcfromtimestamp(hydraulics.index.values[-1].tolist() / 1e9) + observation_start = datetime.utcfromtimestamp( + hydraulics.index.values[0].tolist() / 1e9) + observation_end = datetime.utcfromtimestamp( + hydraulics.index.values[-1].tolist() / 1e9) # rename columns to match hydws - #hydraulics = hydraulics.drop('time', axis=1).rename( + # hydraulics = hydraulics.drop('time', axis=1).rename( # columns={'pressure': 'toppressure'}) hydraulics_plan.index.names = ['datetime'] @@ -67,7 +76,6 @@ def main(): parser[section_id].hydraulics = hydraulics_plan scenario_hydjson = parser.to_json() - ############################################################ # SEISMICITY @@ -85,10 +93,13 @@ def main(): transformer = CoordinateTransformer(proj_string, *top_section_swiss) # Read seismics, transform to local coordinate system - #catalog_df = pd.read_hdf(Data_Path + 'TrainingTarget.h5') + # catalog_df = pd.read_hdf(Data_Path + 'TrainingTarget.h5') catalog_df = pd.read_hdf(Data_Path + 'Cat_FORGE22.h5') - catalog_df = catalog_df[catalog_df.index.values < np.array(observation_end).astype('datetime64')] - catalog_df = catalog_df[catalog_df.index.values >= np.array(observation_start).astype('datetime64')] + catalog_df = catalog_df[catalog_df.index.values + < np.array(observation_end).astype('datetime64')] + catalog_df = \ + catalog_df[catalog_df.index.values + >= np.array(observation_start).astype('datetime64')] catalog_df['Xs'], catalog_df['Ys'], catalog_df['Zs'] = \ transformer.from_local_coords(catalog_df.Xs.values, catalog_df.Ys.values, @@ -152,8 +163,8 @@ def main(): results = model.predict() - #results = model.parse_results(results) - #catalog_df_full = pd.read_hdf(Data_Path + 'seismic_Data_full_located.h5') + # results = model.parse_results(results) + # catalog_df_full = pd.read_hdf(Data_Path + 'seismic_Data_full_located.h5') catalog_df_full = pd.read_hdf(Data_Path + 'Cat_FORGE22.h5') catalog_df_full.index.names = ['datetime'] catalog_df_full = catalog_df_full.reset_index().rename( @@ -164,19 +175,22 @@ def main(): 'mag': 'magnitude', 'datetime': 'time'}) - #catalog_df_full.rename(columns={"X": "Xs", "Y": "Ys", "Depth": "Zs", "MomMag": "magnitude", "Time(UTC)": "time"}, inplace=True) - #catalog_df_full.rename(columns={"mag": "magnitude"}, inplace=True) + # catalog_df_full.rename( + # columns={"X": "Xs", "Y": "Ys", "Depth": "Zs", + # "MomMag": "magnitude", "Time(UTC)": "time"}, + # inplace=True) + # catalog_df_full.rename(columns={"mag": "magnitude"}, inplace=True) catalog_df_full = catalog_df_full[catalog_df_full["magnitude"] > -5] catalog_df_full.set_index("time", inplace=True) - catalog_df_full = catalog_df_full[catalog_df_full.index >= hydraulics.index.values[0]] + catalog_df_full = catalog_df_full[catalog_df_full.index + >= hydraulics.index.values[0]] catalog_df_full = catalog_df_full[catalog_df_full.index <= forecast_end] catalog_df_full = catalog_df_full.reset_index() - plotting_results(model,results,catalog_df_full) + plotting_results(model, results, catalog_df_full) return model if __name__ == "__main__": result = main() - diff --git a/run_model.py b/run_model.py index 1ec5c6e1ad49a7c36abf7f3e67b9a580dc9a7039..a2caacd544954b67b559e24d5e5fffa312300ccb 100644 --- a/run_model.py +++ b/run_model.py @@ -4,12 +4,12 @@ from datetime import datetime import pandas as pd from hydws.parser import BoreholeHydraulics -from seismostats.seismicity.catalog import Catalog +from seismostats import Catalog from seismostats.utils.coordinates import CoordinateTransformer from shapely import Polygon from shapely.wkt import dumps -from em1.ramsis import run_em1, run_em1_validated # noqa +from em1.hermes import run_em1, run_em1_validated # noqa # Set up logger date_strftime_format = "%d-%m-%y %H:%M:%S" @@ -115,7 +115,7 @@ def main(): 'reference_point': top_section_swiss, # always in WGS84 'local_proj_string': proj_string, 'epoch_duration': 600, - 'n_phases':8 + 'n_phases': 8 } } diff --git a/setup.cfg b/setup.cfg index 5de7d5d889c11702b743b023a8d68562e36f1cd5..3c898943f26233f2291c2dd6f69048e3068fb3d0 100644 --- a/setup.cfg +++ b/setup.cfg @@ -26,7 +26,7 @@ install_requires = pandas pyarrow pydantic - ramsis-model[hydws] @ git+https://gitlab.seismo.ethz.ch/indu/ramsis-model.git + hermes-model[hydws] @ git+https://gitlab.seismo.ethz.ch/indu/hermes-model.git scikit-learn seismostats @ git+https://github.com/swiss-seismological-service/SeismoStats.git tables