dt-geo/DeterministicMmaxEstimates
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
d9ededf5bc
DeterministicMmaxEstimates/Mmax.py

714 lines
35 KiB
Python
Raw Normal View History

Populated app repository with code and licence
2024-10-22 13:46:54 +02:00
import sys
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# import os
import scipy.io
import logging
from geopy.distance import geodesic
from geopy.point import Point
def main(Input_catalog, Input_injection_rate, time_win_in_hours, time_step_in_hour, time_win_type,
End_time, ev_limit, Inpar, time_inj, time_shut_in, time_big_ev, Model_index,
Mc, Mu, G, ssd, C, b_value_type, cl, mag_name):
"""
Main function for application: MAXIMUM_MAGNITUDE_DETERMINISTIC_MODELS
Arguments:
Input catalog: path to input file of type 'catalog'
Input injection rate: path to input file of type 'injection_rate'
**kwargs: Model paramters.
Returns:
'PLOT_Mmax_param': Plot of all results (also check 'application.log' for more info)
'DATA_Mmax_param': Results with 'csv' format
'application.log': Logging file
"""
f_indx = Model_index
b_method = b_value_type
Plot_flag = 1
# Setting up the logfile--------------------------------
logging.basicConfig(filename="application.log",
filemode='a',
format='%(asctime)s,%(msecs)d %(name)s %(levelname)s %(message)s',
datefmt='%H:%M:%S',
level=logging.INFO)
logger = logging.getLogger(__name__)
# Importing utilities ----------------------
logger.info("Import utilities")
from util.Find_idx4Time import Find_idx4Time
from util.CandidateEventsTS import CandidateEventsTS
from util.M_max_models import M_max_models
def latlon_to_enu(lat, lon, alt):
lat0 = np.mean(lat)
lon0 = np.mean(lon)
alt0 = 0
# Reference point (origin)
origin = Point(lat0, lon0, alt0)
east = np.zeros_like(lon)
north = np.zeros_like(lat)
up = np.zeros_like(alt)
for i in range(len(lon)):
# Target point
target = Point(lat[i], lon[i], alt[i])
# Calculate East-North-Up
east[i] = geodesic((lat0, lon0), (lat0, lon[i])).meters
if lon[i] < lon0:
east[i] = -east[i]
north[i] = geodesic((lat0, lon0), (lat[i], lon0)).meters
if lat[i] < lat0:
north[i] = -north[i]
up = alt - alt0
return east, north, up
# Importing data
logger.info("Import data")
mat = scipy.io.loadmat(Input_catalog)
Cat_structure = mat['Catalog']
Cat_id, Cat_t, Cat_m = [], [], []
Cat_x, Cat_y, Cat_z = [], [], []
Cat_lat, Cat_lon, Cat_elv, Cat_depth = [], [], [], []
for i in range(1,Cat_structure.shape[1]):
if Cat_structure.T[i][0][0][0] == 'X':
Cat_x = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_x)):
raise ValueError("Catalog-X contains infinite value")
if Cat_structure.T[i][0][3][0] == 'km':
Cat_x *= 1000
if Cat_structure.T[i][0][0][0] == 'Lat':
Cat_lat = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_lat)):
raise ValueError("Catalog-Lat contains infinite value")
if Cat_structure.T[i][0][0][0] == 'Y':
Cat_y = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_y)):
raise ValueError("Catalog-Y contains infinite value")
if Cat_structure.T[i][0][3][0] == 'km':
Cat_y *= 1000
if Cat_structure.T[i][0][0][0] == 'Long':
Cat_lon = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_lon)):
raise ValueError("Catalog-Long contains infinite value")
if Cat_structure.T[i][0][0][0] == 'Z':
Cat_z = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_z)):
raise ValueError("Catalog-Z contains infinite value")
if Cat_structure.T[i][0][3][0] == 'km':
Cat_z *= 1000
if Cat_structure.T[i][0][0][0] == 'Depth':
Cat_depth = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_depth)):
raise ValueError("Catalog-Depth contains infinite value")
if Cat_structure.T[i][0][3][0] == 'km':
Cat_depth *= 1000
if Cat_structure.T[i][0][0][0] == 'Elevation':
Cat_elv = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_elv)):
raise ValueError("Catalog-Elevation contains infinite value")
if Cat_structure.T[i][0][3][0] == 'km':
Cat_elv *= 1000
if Cat_structure.T[i][0][0][0] == 'Time':
Cat_t = Cat_structure.T[i][0][2]
if not all(np.isfinite(Cat_t)):
raise ValueError("Catalog-Time contains infinite value")
if Cat_structure.T[i][0][0][0] == mag_name:
Cat_m = Cat_structure.T[i][0][2]
# if not all(np.isfinite(Cat_m)):
if np.argwhere(all(np.isfinite(Cat_m))):
raise ValueError("Catalog-Magnitude contains infinite value")
if any(Cat_x):
Cat_id = np.linspace(0,Cat_x.shape[0],Cat_x.shape[0]).reshape((Cat_x.shape[0],1))
arg = (Cat_id, Cat_x, Cat_y, Cat_z, Cat_t, Cat_m)
Cat = np.concatenate(arg, axis=1)
elif any(Cat_lat):
if any(Cat_elv):
Cat_x, Cat_y, Cat_z = latlon_to_enu(Cat_lat, Cat_lon, Cat_elv)
elif any(Cat_depth):
Cat_x, Cat_y, Cat_z = latlon_to_enu(Cat_lat, Cat_lon, Cat_depth)
else:
raise ValueError("Catalog Depth or Elevation is not available")
Cat_id = np.linspace(0,Cat_x.shape[0],Cat_x.shape[0]).reshape((Cat_x.shape[0],1))
arg = (Cat_id, Cat_x, Cat_y, Cat_z, Cat_t, Cat_m)
Cat = np.concatenate(arg, axis=1)
else:
raise ValueError("Catalog data are not available")
mat = scipy.io.loadmat(Input_injection_rate)
Inj_structure = mat['d']
if 'Date' in Inj_structure.dtype.names:
inj_date = Inj_structure['Date'][0,0]
inj_rate = Inj_structure['Injection_rate'][0,0]
Hyd = np.concatenate((inj_date,inj_rate), axis=1)
else:
raise ValueError("Injection data are not available")
if Cat[0,4]>np.max(Hyd[:,0]) or Hyd[0,0]>np.max(Cat[:,4]):
raise ValueError('Catalog and injection data do not have time coverage!')
if Hyd[0,0] < Cat[0,4]:
same_time_idx = Find_idx4Time(Hyd[:,0], Cat[0,4])
Hyd = Hyd[same_time_idx:,:]
Hyd[:,0] = (Hyd[:,0] - Cat[0,4])*24*3600
Cat[:,4] = (Cat[:,4] - Cat[0,4])*24*3600
logger.info('Start of the computations is based on the time of catalog data.')
# Model dictionary
Feat_dic = {
'indx' :[0 ,1 ,2 ,3 ,4 ,5 ],
'Full_names':['All_M_max' ,'McGarr' ,'Hallo' ,'Li' ,'van-der-Elst','Shapiro' ],
'Short_name':['max_all' , 'max_mcg', 'max_hlo', 'max_li', 'max_vde' , 'max_shp'],
'f_num' :[5 ,4 ,4 ,1 ,6 ,5 ],
'Param' :[{'Mc': 0.8, 'b_method': ['b'], 'cl': [0.37], 'Inpar': ['dv','d_num'], 'Mu':0.6, 'G': 35*10**(9), 'ssd': 3*10**6, 'C': 0.95, 'ev_limit': 20, 'num_bootstraps': 100}]
}
Feat_dic['Param'][0]['Inpar'] = Inpar
Feat_dic['Param'][0]['ev_limit'] = ev_limit
num_bootstraps = 100 # Number of bootstraping for standard error computation
Feat_dic['Param'][0]['num_bootstraps'] = num_bootstraps
if f_indx in [0,1,2,4]:
if Mc < np.min(Cat[:,-1]) or Mc > np.max(Cat[:,-1]):
raise ValueError("Completeness magnitude (Mc) is out of magnitude range")
Feat_dic['Param'][0]['Mc'] = Mc
if f_indx in [0,1,2]:
if Mu < 0.2 or Mu > 0.8:
raise ValueError("Friction coefficient (Mu) must be between [0.2, 0.8]")
Feat_dic['Param'][0]['Mu'] = Mu
if f_indx in [0,1,2,3]:
if G < 1*10**(9) or G > 100*10**(9):
raise ValueError("Shear modulus of reservoir rock (G) must be between [1, 100] GPa")
Feat_dic['Param'][0]['G'] = G
if f_indx in [0,5]:
if ssd < 0.1*10**(6) or ssd > 100*10**(6):
raise ValueError("Static stress drop (ssd) must be between [0.1, 100] MPa")
Feat_dic['Param'][0]['ssd'] = ssd
if f_indx in [0,5]:
if C < 0.5 or C > 5:
raise ValueError("Geometrical constant (C) of Shaprio's model must be between [0.5, 5]")
Feat_dic['Param'][0]['C'] = C
if f_indx in [0,1,2,4]:
Feat_dic['Param'][0]['b_method'] = b_method
if f_indx in [0,4]:
for cl_i in cl:
if cl_i < 0 or cl_i > 1:
raise ValueError("Confidence level (cl) of van der Elst model must be between [0, 1]")
Feat_dic['Param'][0]['cl'] = cl
# Setting up based on the config and model dic --------------
ModelClass = M_max_models()
Model_name = Feat_dic['Full_names'][f_indx]
ModelClass.f_name = Feat_dic['Short_name'][f_indx]
f_num = Feat_dic['f_num'][f_indx]
if Feat_dic['Param'][0]['Mc']:
ModelClass.Mc = Mc
else:
Mc = np.min(Cat[:,-1])
ModelClass.Mc = Mc
time_win = time_win_in_hours*3600 # in sec
ModelClass.time_win = time_win
if f_indx == 0:
# Only first b_methods is used
Feat_dic['Param'][0]['b_method'] = b_method[0]
ModelClass.b_method = Feat_dic['Param'][0]['b_method']
logger.info(f"All models are based on b_method: { b_method[0]}")
Feat_dic['Param'][0]['cl'] = [cl[0]]
ModelClass.cl = Feat_dic['Param'][0]['cl']
logger.info(f"All models are based on cl: { cl[0]}")
ModelClass.Mu = Feat_dic['Param'][0]['Mu']
ModelClass.G = Feat_dic['Param'][0]['G']
ModelClass.ssd = Feat_dic['Param'][0]['ssd']
ModelClass.C = Feat_dic['Param'][0]['C']
ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps']
if f_indx == 1 or f_indx == 2:
ModelClass.b_method = Feat_dic['Param'][0]['b_method']
ModelClass.Mu = Feat_dic['Param'][0]['Mu']
ModelClass.G = Feat_dic['Param'][0]['G']
ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps']
Feat_dic['Param'][0]['cl'] = [None]
if f_indx == 3:
Feat_dic['Param'][0]['b_method'] = [None]
ModelClass.G = Feat_dic['Param'][0]['G']
Feat_dic['Param'][0]['cl'] = [None]
if f_indx == 4:
ModelClass.b_method = Feat_dic['Param'][0]['b_method']
ModelClass.num_bootstraps = Feat_dic['Param'][0]['num_bootstraps']
ModelClass.G = Feat_dic['Param'][0]['G']
ModelClass.cl = Feat_dic['Param'][0]['cl']
if f_indx == 5:
ModelClass.ssd = Feat_dic['Param'][0]['ssd']
ModelClass.C = Feat_dic['Param'][0]['C']
Feat_dic['Param'][0]['b_method'] = [None]
Feat_dic['Param'][0]['cl'] = [None]
# Output dictionary --------------------------------
Output_dict = {
'Type' :['idx', 'Time[day]'],
'label' :[None, None],
'b_method' :[None, None],
'cl' :[None, None],
}
c_out = 2
if any(Feat_dic['Param'][0]['b_method']) and any(Feat_dic['Param'][0]['cl']):
for i in range(len(Feat_dic['Param'][0]['b_method'])):
for j in range(len(Feat_dic['Param'][0]['cl'])):
if f_indx == 0: # f_index == 0
for i in Feat_dic['Full_names'][1:]:
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(i)
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
c_out += 1
elif j == 0: # f_index == 4
Output_dict['Type'].append('b_value')
Output_dict['label'].append('b_value')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('b_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Seismogenic Index')
Output_dict['label'].append('SI')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('si_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(Model_name)
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(Feat_dic['Param'][0]['cl'][j])
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('M_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
c_out += 6
else: # f_index == 4
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(Model_name)
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(Feat_dic['Param'][0]['cl'][j])
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('M_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
c_out += 2
elif any(Feat_dic['Param'][0]['b_method']): # f_index == 1, 2
for i in range(len(Feat_dic['Param'][0]['b_method'])):
Output_dict['Type'].append('b_value')
Output_dict['label'].append('b_value')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('b_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(Model_name)
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
Output_dict['Type'].append('Standard Error')
Output_dict['label'].append('M_std_err')
Output_dict['b_method'].append(Feat_dic['Param'][0]['b_method'][i])
Output_dict['cl'].append(None)
c_out += 4
elif f_indx == 5: # f_index == 5
for i in ['L(max)','L(int)','L(min)', 'L(avg)']:
Output_dict['Type'].append('Length[m]')
Output_dict['label'].append(i)
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(Model_name)
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
c_out += 5
else: # f_index == 3
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append(Model_name)
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
c_out += 1
# Add True maximum magnitude to the outputs
Output_dict['Type'].append('Maximum magnitude')
Output_dict['label'].append('True Max-Mag')
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
c_out += 1
# if any(Feat_dic['Param'][0]['Inpar']): # Input parameters
if Inpar:
for i in Feat_dic['Param'][0]['Inpar']:
if i == 'd_num':
Output_dict['Type'].append('Number of events')
Output_dict['label'].append('Ev_Num')
elif i == 'dv':
Output_dict['Type'].append('Volume[m3]')
Output_dict['label'].append('Vol')
elif i == 'Mo':
Output_dict['Type'].append('Seismic moment')
Output_dict['label'].append('Mo')
elif i == 'SER':
Output_dict['Type'].append('Seismic efficiency ratio')
Output_dict['label'].append('SER')
Output_dict['b_method'].append(None)
Output_dict['cl'].append(None)
c_out += 1
# Functions ------------------------------
# Computing in extending time or time window
def Feature_in_Time(In_arr):
SER_l = 0
SER_c = 0
i_row = 0
for i in np.arange(1,Times):
# Defining time window and data
if time_win_type == 0:
ModelClass.time_win = i*time_step
candidate_events = CandidateEventsTS(Cat, i*time_step, None, i*time_step)
else:
candidate_events = CandidateEventsTS(Cat, i*time_step, None, time_win)
data = candidate_events.filter_data()
# USER: Check if cumulative dv is correctly computed !!!
end_time_indx = Find_idx4Time(Hyd[:,0], i*time_step)
# ModelClass.dv = np.sum(Hyd[:end_time_indx,1])
ModelClass.dv = np.trapz(Hyd[:end_time_indx,1], Hyd[:end_time_indx,0])/60
if len(data) > Feat_dic['Param'][0]['ev_limit'] and ModelClass.dv > 0:
ModelClass.Mo = np.sum(10**(1.5*data[:end_time_indx,5]+9.1))
if f_indx != 5: # Parameters have not been assinged for f_index == 5
SER_c = ModelClass.Mo/(2.0*ModelClass.G*ModelClass.dv)
SER_l = np.max((SER_c, SER_l))
ModelClass.SER = SER_l
ModelClass.data = data
In_arr[i_row,0] = i # index
In_arr[i_row,1] = i*time_step # Currecnt time
c = 2
if any(Feat_dic['Param'][0]['b_method']) and any(Feat_dic['Param'][0]['cl']):
for ii in range(len(Feat_dic['Param'][0]['b_method'])):
ModelClass.b_method = Feat_dic['Param'][0]['b_method'][ii]
for jj in range(len(Feat_dic['Param'][0]['cl'])): # f_index == 4
ModelClass.cl = Feat_dic['Param'][0]['cl'][jj]
Out = ModelClass.ComputeModel()
if jj == 0:
In_arr[i_row,c:c+f_num] = Out
c += f_num
else:
In_arr[i_row,c:c+2] = Out[-2:]
c += 2
elif any(Feat_dic['Param'][0]['b_method']): # f_index == 1, 2
for ii in range(len(Feat_dic['Param'][0]['b_method'])):
ModelClass.b_method = Feat_dic['Param'][0]['b_method'][ii]
In_arr[i_row,c:c+f_num] = ModelClass.ComputeModel()
c += f_num
else: # f_index = 0, 3, 5
In_arr[i_row,c:c+f_num] = ModelClass.ComputeModel()
c += f_num
# Compute true maximum magnitude in the data
In_arr[i_row,c] = np.max(data[:end_time_indx,5])
c += 1
if Inpar:
for ii in range(len(Feat_dic['Param'][0]['Inpar'])): # Add input parameters
if Feat_dic['Param'][0]['Inpar'][ii] == 'd_num':
In_arr[i_row,c+ii] = data.shape[0]
elif Feat_dic['Param'][0]['Inpar'][ii] == 'dv':
In_arr[i_row,c+ii] = ModelClass.dv
elif Feat_dic['Param'][0]['Inpar'][ii] == 'Mo':
In_arr[i_row,c+ii] = ModelClass.Mo
elif Feat_dic['Param'][0]['Inpar'][ii] == 'SER':
if ModelClass.SER:
In_arr[i_row,c+ii] = ModelClass.SER
i_row += 1
return In_arr[:i_row,:]
# Plotting models and parameters
def Plot_feature(Model_Param_array,Output_dict):
myVars = locals()
# Function for findin min-max of all similar parameters
def Extermom4All(Model_Param_array, itr_loc):
Mat1D = np.reshape(Model_Param_array[:,itr_loc], -1)
NotNone = np.isfinite(Mat1D)
if np.min(Mat1D[NotNone])>0:
return [np.min(Mat1D[NotNone])*0.95, np.max(Mat1D[NotNone])*1.05]
elif np.min(Mat1D[NotNone])<0 and np.max(Mat1D[NotNone])>0:
return [np.min(Mat1D[NotNone])*1.05, np.max(Mat1D[NotNone])*1.05]
elif np.max(Mat1D[NotNone])<0:
return [np.min(Mat1D[NotNone])*1.05, np.max(Mat1D[NotNone])*0.95]
# Function for setting relevant lagends in the plot
def Legend_label(loc):
l = Output_dict_c['label'][loc]
if Output_dict_c['b_method'][loc]:
if Output_dict_c['cl'][loc]:
l+='('+Output_dict_c['b_method'][loc]+', cl='+str(Output_dict_c['cl'][loc])+')'
else:
l+='('+Output_dict_c['b_method'][loc]+')'
return l
c_NotNone = [] # Removing all parameters with None or constant value
for i in range(Model_Param_array.shape[1]):
NotNone = np.isfinite(Model_Param_array[:,i])
Eq_value = np.mean(Model_Param_array[:,i])
if any(NotNone) and Eq_value != Model_Param_array[0,i]:
c_NotNone.append(i)
else:
logger.info(f"No-PLOT: All values of {Output_dict['Type'][i]} are {Model_Param_array[0,i]}!")
if len(c_NotNone) > 1:
Model_Param_array = Model_Param_array[:,c_NotNone]
# New output dictionary based on valid parameters for plotting
Output_dict_c = {'Type':[], 'label':[], 'b_method':[], 'cl':[]}
for i in range(len(c_NotNone)):
Output_dict_c['Type'].append(Output_dict['Type'][c_NotNone[i]])
Output_dict_c['label'].append(Output_dict['label'][c_NotNone[i]])
Output_dict_c['b_method'].append(Output_dict['b_method'][c_NotNone[i]])
Output_dict_c['cl'].append(Output_dict['cl'][c_NotNone[i]])
coloring=['blue','g','r','c','m','y',
'brown', 'darkolivegreen', 'teal', 'steelblue', 'slateblue',
'purple', 'darksalmon', '#c5b0d5', '#c49c94',
'#e377c2', '#f7b6d2', '#7f7f7f', '#c7c7c7', '#bcbd22', '#dbdb8d',
'#17becf', '#9edae5']
# All parameters to be plotted
All_vars = Output_dict_c['Type'][2:]
Uniqe_var = list(dict.fromkeys([s for s in All_vars if 'Standard Error' not in s])) #list(set(All_vars))
# defining handels and labels to make final legend
All_handels = ['p0']
for i in range(1,len(All_vars)):
All_handels.append('p'+str(i))
handels = []
labels = []
# itr_loc: location of paramteres with similar type
itr_loc = np.where(np.array(All_vars) == Uniqe_var[0])[0]+2
fig, myVars[Output_dict_c['Type'][0]] = plt.subplots(1,1,figsize=(8+int(len(All_vars)/3),6))
fig.subplots_adjust(right=1-len(Uniqe_var)*0.09)
if Output_dict_c['label'][itr_loc[0]] == 'True Max-Mag': # plot with dash-line
myVars[All_handels[itr_loc[0]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]], ls='--')
else:
myVars[All_handels[itr_loc[0]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]])
handels.append(All_handels[itr_loc[0]-2])
labels.append(Legend_label(itr_loc[0]))
myVars[Output_dict_c['Type'][0]].set_ylabel(Output_dict_c['Type'][itr_loc[0]])
myVars[Output_dict_c['Type'][0]].set_ylim(Extermom4All(Model_Param_array, itr_loc)[0], Extermom4All(Model_Param_array, itr_loc)[1])
myVars[Output_dict_c['Type'][0]].set_xlabel('Day (From start of the recording)')
if End_time:
myVars[Output_dict_c['Type'][0]].set_xlim(0,End_time)
# Plotting statndard error (if exists)
if Output_dict_c['Type'][itr_loc[0]+1] == 'Standard Error':
myVars[Output_dict_c['Type'][0]].fill_between(Model_Param_array[:,1]/24/3600,
Model_Param_array[:,itr_loc[0]] - Model_Param_array[:,itr_loc[0]+1],
Model_Param_array[:,itr_loc[0]] + Model_Param_array[:,itr_loc[0]+1], color= coloring[itr_loc[0]], alpha=0.1)
# Plotting similar parameters on one axis
for j in range(1,len(itr_loc)):
if Output_dict_c['label'][itr_loc[j]] == 'True Max-Mag': # plot with dash-line
myVars[All_handels[itr_loc[j]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]], ls='--')
else:
myVars[All_handels[itr_loc[j]-2]], = myVars[Output_dict_c['Type'][0]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]])
handels.append(All_handels[itr_loc[j]-2])
labels.append(Legend_label(itr_loc[j]))
# Plotting statndard error (if exists)
if itr_loc[j]+1 <= len(itr_loc) and Output_dict_c['Type'][itr_loc[j]+1] == 'Standard Error':
myVars[Output_dict_c['Type'][0]].fill_between(Model_Param_array[:,1]/24/3600,
Model_Param_array[:,itr_loc[j]] - Model_Param_array[:,itr_loc[j]+1],
Model_Param_array[:,itr_loc[j]] + Model_Param_array[:,itr_loc[j]+1], color= coloring[itr_loc[j]], alpha=0.1)
first_itr = 0
# Check if there is any more parameter to be plotted in second axes
# The procedure is similar to last plots.
if len(Uniqe_var) > 1:
for i in range(1,len(Uniqe_var)):
itr_loc = np.where(np.array(All_vars) == Uniqe_var[i])[0]+2
myVars[Uniqe_var[i]] = myVars[Output_dict_c['Type'][0]].twinx()
# if it is third or more axis, make a distance between them
if first_itr == 0:
first_itr += 1
set_right = 1
else:
set_right = 1 + first_itr*0.2
first_itr += 1
myVars[Uniqe_var[i]].spines.right.set_position(("axes", set_right))
if Output_dict_c['label'][itr_loc[0]] == 'True Max-Mag': # plot with dash-line
myVars[All_handels[itr_loc[0]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]], ls='--')
else:
myVars[All_handels[itr_loc[0]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[0]], c= coloring[itr_loc[0]])
handels.append(All_handels[itr_loc[0]-2])
labels.append(Legend_label(itr_loc[0]))
myVars[Uniqe_var[i]].set_ylabel(Output_dict_c['Type'][itr_loc[0]])
myVars[Uniqe_var[i]].yaxis.label.set_color(coloring[itr_loc[0]])
myVars[Uniqe_var[i]].spines["right"].set_edgecolor(coloring[itr_loc[0]])
myVars[Uniqe_var[i]].tick_params(axis='y', colors= coloring[itr_loc[0]])
myVars[Uniqe_var[i]].set_ylim(Extermom4All(Model_Param_array, itr_loc)[0], Extermom4All(Model_Param_array, itr_loc)[1])
if itr_loc[0]+1 < len(Output_dict_c['Type']) and Output_dict_c['Type'][itr_loc[0]+1] == 'Standard Error':
myVars[Uniqe_var[i]].fill_between(Model_Param_array[:,1]/24/3600,
Model_Param_array[:,itr_loc[0]] - Model_Param_array[:,itr_loc[0]+1],
Model_Param_array[:,itr_loc[0]] + Model_Param_array[:,itr_loc[0]+1], color= coloring[itr_loc[0]], alpha=0.1)
for j in range(1,len(itr_loc)):
if Output_dict_c['label'][itr_loc[j]] == 'True Max-Mag': # plot with dash-line
myVars[All_handels[itr_loc[j]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]], ls = '--')
else:
myVars[All_handels[itr_loc[j]-2]], = myVars[Uniqe_var[i]].plot(Model_Param_array[:,1]/24/3600, Model_Param_array[:,itr_loc[j]], c= coloring[itr_loc[j]])
handels.append(All_handels[itr_loc[j]-2])
labels.append(Legend_label(itr_loc[j]))
if Output_dict_c['Type'][itr_loc[j]+1] == 'Standard Error':
myVars[Uniqe_var[i]].fill_between(Model_Param_array[:,1]/24/3600,
Model_Param_array[:,itr_loc[j]] - Model_Param_array[:,itr_loc[j]+1],
Model_Param_array[:,itr_loc[j]] + Model_Param_array[:,itr_loc[j]+1], color= coloring[itr_loc[j]], alpha=0.1)
# If there are timing, plot them as vertical lines
if time_inj:
myVars['l1'], = plt.plot([time_inj,time_inj], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='--', c='k')
handels.append('l1')
labels.append('Start-inj')
if time_shut_in:
myVars['l2'], = plt.plot([time_shut_in,time_shut_in], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='-.', c='k')
handels.append('l2')
labels.append('Shut-in')
if time_big_ev:
myVars['l3'], = plt.plot([time_big_ev,time_big_ev], [Extermom4All(Model_Param_array, itr_loc)[0],Extermom4All(Model_Param_array, itr_loc)[1]], ls='dotted', c='k')
handels.append('l3')
labels.append('Large-Ev')
box = myVars[Output_dict_c['Type'][0]].get_position()
if len(handels) < 6:
myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])
plt.legend([myVars[ii] for ii in handels], labels, loc='upper center',
bbox_to_anchor=(0.5+0.06*first_itr, -0.15), fancybox=True, shadow=True, ncol=len(handels))
elif len(handels) < 13:
myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.04*int(len(handels)/2),
box.width, box.height * (1 - 0.04*int(len(handels)/2))])
plt.legend([myVars[ii] for ii in handels], labels, loc='upper center',
bbox_to_anchor=(0.5+0.1*first_itr, -0.04*int(len(handels)/2)), fancybox=True, shadow=True, ncol=int(len(handels)/2)+1, handleheight=2)
else:
myVars[Output_dict_c['Type'][0]].set_position([box.x0, box.y0 + box.height * 0.04*int(len(handels)/2),
box.width, box.height * (1 - 0.04*int(len(handels)/2))])
plt.legend([myVars[ii] for ii in handels], labels, loc='upper center',
bbox_to_anchor=(0.6+0.1*first_itr, -0.04*int(len(handels)/2)), fancybox=True, shadow=True, ncol=int(len(handels)/2)+1, handleheight=2)
plt.title(Model_name)
# plt.savefig(cwd+'/Results/'+Model_name+'.png', dpi=300)
plt.savefig('PLOT_Mmax_param.png', dpi=300)
# plt.show()
# Run functions based on the configurations -------------------
# Computing model
if time_step_in_hour > time_win_in_hours:
raise ValueError('Time steps should be <= time window.')
if (time_win_in_hours+3*time_step_in_hour)*3600 > np.max(Cat[:,4]):
raise ValueError('Time window and time steps are like that no more than three computations is possible. Use smaller value for either time window or time step.')
time_step = time_step_in_hour*3600
if End_time:
if End_time*24*3600 > np.max(Cat[:,4])+24*3600:
raise ValueError(f'End_time is longer than the maximum time of catalog, which is {np.ceil(np.max(Cat[:,4])/24/3600)} day!')
Times = int(np.floor((End_time*24*3600 - Cat[0,4]) / time_step))
else:
Times = int(np.floor((np.max(Cat[:,4]) - Cat[0,4]) / time_step))
Model_Param_array = np.zeros((Times,c_out))
logger.info("Computing input parameter(s) and the model(s)")
Model_Param_array = Feature_in_Time(Model_Param_array)
# Plotting
if Plot_flag > 0:
if Model_Param_array.any():
logger.info("Plotting results")
Plot_feature(Model_Param_array, Output_dict)
else:
logger.info("Model_Param_array is empty.")
raise ValueError("No model is generated since no event in all time windows is available due to imporoper values of either time window, minimum limit of events or completeness magnitude!")
# Saving
logger.info("Saving results")
# Save models and parameters in
def OneLineHeader(loc):
l = Output_dict['Type'][loc]
if Output_dict['b_method'][loc]:
if Output_dict['label'][loc] != 'b_value' and Output_dict['label'][loc] != 'b_std_err' and Output_dict['label'][loc] != 'M_std_err':
l += '_'+Output_dict['label'][loc]
else:
l = Output_dict['label'][loc]
if Output_dict['cl'][loc]:
l += '(b_method: '+Output_dict['b_method'][loc]+', q='+str(Output_dict['cl'][loc])+')'
else:
l += '(b_method: '+Output_dict['b_method'][loc]+')'
elif Output_dict['label'][loc]:
l += '('+ Output_dict['label'][loc] +')'
return l
db_df = pd.DataFrame(data = Model_Param_array,
columns = [OneLineHeader(i) for i in range(len(Output_dict['Type']))])
db_df = db_df.round(4)
# db_df.to_excel(cwd+'/Results/'+Full_feat_name+'.xlsx')
db_df.to_csv('DATA_Mmax_param.csv', sep=';', index=False)
# np.savez_compressed(cwd+'/Results/'+Full_feat_name+'.npz', Model_Param_array = Model_Param_array) # change the name!
if __name__=='__main__':
# Input_catalog = 'Data/Cooper_Basin_Catalog_HAB_1_2003_Reprocessed.mat'
# Input_injection_rate = 'Data/Cooper_Basin_HAB_1_2003_Injection_Rate.mat'
# Model_index = 4
# b_value_type = 'TGR'
main(Input_catalog, Input_injection_rate, time_win_in_hours, time_step_in_hour, time_win_type,
End_time, ev_limit, Inpar, time_inj, time_shut_in, time_big_ev, Model_index,
Mc, Mu, G, ssd, C, b_value_type, cl, mag_name)
Reference in New Issue Copy Permalink