Browse Source

Added interactive GM-map; Fixed bug in DCSource generation

GM
lukas_dell 1 year ago
parent
commit
abf3625a1b
  1. 2
      setup.py
  2. 31
      src/apps/gm_mapping.py
  3. 316
      src/apps/interactive_gm_plot.py
  4. 10
      src/config/gm_mapping.cfg
  5. 0
      src/gm/__init__.py
  6. 48
      src/gm/misc.py
  7. 114
      src/gm/networks.py
  8. 301
      src/gm/plot.py
  9. 41
      src/gm/sources.py
  10. BIN
      src/mapping_functions/NN/DCSource_sklearn.bin
  11. BIN
      src/mapping_functions/NN/DCSource_sklearn_model.bin
  12. BIN
      src/mapping_functions/NN/DCSource_tensorflow.bin
  13. BIN
      src/mapping_functions/NN/DCSource_tensorflow_model.h5

2
setup.py

@ -44,7 +44,7 @@ packages = [
'%s' % packname,
'%s.dataset' % packname,
'%s.apps' % packname,
'%s.GM' % packname]
'%s.gm' % packname]
# entry_points = {

31
src/apps/gm_mapping.py

@ -11,10 +11,11 @@ from multiprocessing import Pool, Manager
from pyrocko import gf
import ewrica.sources as ews
import ewrica.GM.networks as GMn
import ewrica.GM.plot as GMp
import ewrica.GM.misc as GMm
import ewrica.GM.sources as GMs
import ewrica.gm.networks as GMn
import ewrica.gm.plot as GMp
import ewrica.gm.misc as GMm
import ewrica.gm.sources as GMs
def parser_func():
@ -133,7 +134,7 @@ def parser_func():
parser.add_argument('--nn-directory', dest='nn_dir',
nargs=1, metavar='dirPATH',
default=['../data/mapping_functions/NN/'],
default=['../mapping_functions/NN/'],
help='''
Defines the directory where all the NN are saved,
according to the specific sources [..]
@ -146,7 +147,7 @@ def parser_func():
''')
args = parser.parse_args()
print(args)
datafile = args.datadir[0] + 'data.yaml'
args.datafile = datafile
if os.path.exists(datafile):
@ -254,7 +255,7 @@ def main(argv):
elif args.method[0] == 'Pyrocko':
engine = gf.LocalEngine(store_superdirs=[args.gf[0]],
store_dirs=[args.gf[0]])
store_dirs=[args.gf[0]])
### Create Waveform targets
waveform_targets = [gf.Target(quantity='displacement',
@ -294,8 +295,11 @@ def main(argv):
### Create dataframe with all values
#############################
chagms = []
weights = []
allGDF = gpd.GeoDataFrame()
for srccnt, gdf in enumerate(srcsDict.values()):
for srccnt, gdf in srcsDict.items():
weights.append(srcs[srccnt].probability)
for im in gdf.columns:
if im == 'geometry':
if srccnt == 0:
@ -309,8 +313,7 @@ def main(argv):
#############################
### Get weighting
#############################
weights = [src.probability for src in srcs]
probSrcCnt = weights.index(max(weights))
probSrcCnt = [src.probability for src in srcs].index(max(weights))
#############################
### Calc (weighted) mean and std per comp and gm/im
@ -322,7 +325,8 @@ def main(argv):
'''
for chagm in chagms:
cols = [col for col in allGDF.columns if chagm in col]
df = allGDF[cols]
df = allGDF[cols].copy(deep=True)
allGDF['%s_mean' % (chagm)] = df.mean(axis=1)
allGDF['%s_std' % (chagm)] = df.std(axis=1)
@ -332,6 +336,11 @@ def main(argv):
allGDF['%s_wmean' % (chagm)] = wmean
allGDF['%s_wstd' % (chagm)] = num.sqrt(wvar)
# ls = df.values.tolist()
# allGDF[chagm] = ls
print(allGDF)
# exit()
#############################
### Save gdf
#############################

316
src/apps/interactive_gm_plot.py

@ -0,0 +1,316 @@
from bokeh.plotting import figure
from bokeh.models import ColorBar, ColumnDataSource, HoverTool,\
TapTool, Select, LinearColorMapper
from bokeh.io import curdoc
from bokeh.layouts import column, row
from bokeh.tile_providers import get_provider, Vendors
from bokeh.colors import RGB
import numpy as num
import pandas as pd
import geopandas as gpd
from matplotlib import cm
import argparse
def parser_func():
parser = argparse.ArgumentParser(description='''
Reads a geojson file and converts it to an interactive map
''',
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('-d', '--datadir', dest='datadir',
nargs=1, metavar='dirPATH',
help='''
Path to the directory where the geojson can be found.
''')
args = parser.parse_args()
print(args)
return args
def wgs84_to_web_mercator(lons, lats):
k = 6378137
x = lons * (k * num.pi / 180.0)
y = num.log(num.tan((90 + lats) * num.pi / 360.0)) * k
return x, y
def get_dataset(src, gm, comp, target='mean'):
rawlons = num.array(gdf.geometry.x)
rawlats = num.array(gdf.geometry.y)
lons, lats = wgs84_to_web_mercator(rawlons, rawlats)
lonsSet = num.sort(list(set(lons)))
latsSet = num.sort(list(set(lats)))
dlon = abs(lonsSet[1] - lonsSet[0]) / 2
dlat = abs(latsSet[1] - latsSet[0]) / 2
xlons = list(num.array([lons + dlon, lons + dlon,
lons - dlon, lons - dlon]).T)
ylats = list(num.array([lats + dlat, lats - dlat,
lats - dlat, lats + dlat]).T)
select_cols = [col for col in src.columns if gm in col]
selectSide = [col for col in select_cols if comp in col]
appendix = ['mean', 'std']
selectMain = [col for col in selectSide for app in appendix if app in col]
dfMain = src[selectMain].copy(deep=True)
dfMain = dfMain.rename(columns={
'%s_%s_mean' % (comp, gm): 'mean',
'%s_%s_std' % (comp, gm): 'std',
'%s_%s_wmean' % (comp, gm): 'wmean',
'%s_%s_wstd' % (comp, gm): 'wstd'})
# dfMain = 10**(dfMain)
dfMain['show'] = dfMain[target]
dfMain['lon'] = xlons
dfMain['lat'] = ylats
dfMain['rawlon'] = rawlons
dfMain['rawlat'] = rawlats
dfSide = src[selectSide].copy(deep=True)
return ColumnDataSource(data=dfMain), ColumnDataSource(data=dfSide)
def make_main_plot(source, title):
tile_provider = get_provider(Vendors.CARTODBPOSITRON)
maxlon = max([max(sublist) for sublist in source.data['lon']])
minlon = min([min(sublist) for sublist in source.data['lon']])
maxlat = max([max(sublist) for sublist in source.data['lat']])
minlat = min([min(sublist) for sublist in source.data['lat']])
# range bounds supplied in web mercator coordinates
plot = figure(plot_width=800, plot_height=800,
tools="pan, wheel_zoom, box_zoom, reset",
toolbar_location='right',
x_range=(minlon, maxlon), y_range=(minlat, maxlat),
x_axis_type="mercator", y_axis_type="mercator")
plot.add_tile(tile_provider)
plot.title.text = title
plot.patches(xs='lon', ys='lat', source=source,
color={'field': 'show', 'transform': color_mapper},
hover_line_color="black", fill_alpha=0.7,
hover_fill_alpha=1.0,)
plot.add_layout(color_bar, 'below')
my_hover = HoverTool()
my_hover.tooltips = [
("index", "$index"),
('Lon', '@rawlon'),
('Lat', '@rawlat'),
('mean', '@mean'),
('std', '@std'),
('wmean', '@wmean'),
('wstd', '@wstd')]
# my_hover.point_policy = 'snap_to_data'
plot.add_tools(my_hover, TapTool())
plot.xaxis.axis_label = 'Longitude'
plot.yaxis.axis_label = "Latitude"
plot.axis.axis_label_text_font_style = "bold"
plot.grid.grid_line_alpha = 0.3
return plot
def update_plot(attrname, old, new):
gm = gmSelect.value
comp = compSelect.value
target = targetSelect.value
srcMain, srcSide = get_dataset(gdf, gm, comp, target)
sourceMain.data.update(srcMain.data)
plot.title.text = '%s_%s_%s' % (comp, gm, target)
if target in ['std', 'wstd']:
color_bar.color_mapper.palette = coolwarm_palette
color_bar.color_mapper.low = 0
color_bar.color_mapper.high = 1
else:
color_bar.color_mapper.palette = 'Viridis256'
color_bar.color_mapper.low = -1
color_bar.color_mapper.high = 2.3
srcHist = get_source_side(srcSide, gidx)
sourceHist.data.update(srcHist.data)
###
def get_source_side(source, idx):
d = pd.DataFrame(source.data).drop(columns='index')
appendix = ['mean', 'std']
statCols = [col for col in d.columns for app in appendix if app in col]
statData = d[statCols]
histData = d.drop(columns=statCols)
histData = list(histData.iloc[idx])
statData = statData.iloc[idx]
hist, edges = num.histogram(histData, bins=35, range=(-1, 2.5))
histDf = pd.DataFrame({'column': hist,
"left": edges[:-1],
"right": edges[1:]})
histDf["interval"] = ["%0.2f to %0.2f" % (left, right)
for left, right in zip(histDf["left"], histDf["right"])]
histDf = pd.DataFrame(histDf)
statList = list(statData.index)
lslen = len(histDf['column'])
histDf['mean'] = [statData[stat] for stat in statList
if stat.endswith('_mean')] * lslen
histDf['wmean'] = [statData[stat] for stat in statList
if stat.endswith('_wmean')] * lslen
histDf['std'] = [statData[stat] for stat in statList
if stat.endswith('_std')] * lslen
histDf['wstd'] = [statData[stat] for stat in statList
if stat.endswith('_wstd')] * lslen
histDf['ys'] = num.linspace(-2, 0, lslen)
histDf['mean-std'] = histDf['mean'] - histDf['std']
histDf['mean+std'] = histDf['mean'] + histDf['std']
histDf['wmean-wstd'] = histDf['wmean'] - histDf['wstd']
histDf['wmean+wstd'] = histDf['wmean'] + histDf['wstd']
return ColumnDataSource(data=histDf)
def make_hist_plot(sourceFin, title, index):
plotHist = figure(plot_width=600, plot_height=800, title=title,
x_axis_label="Val", y_axis_label="Count")
plotHist.quad(bottom=0, top='column', left="left", fill_color='black',
line_color='black',
right="right", source=sourceFin, fill_alpha=0.7, hover_fill_alpha=1.0,
name='hist')
plotHist.line(x='mean', y='ys', source=sourceFin,
line_width=3, color='red', legend_label='Mean')
plotHist.harea(x1='mean+std', x2='mean-std', y='ys', source=sourceFin,
color='red', alpha=0.2, legend_label='Std')
plotHist.line(x='wmean', y='ys', source=sourceFin,
line_width=3, color='green', legend_label='WMean')
plotHist.harea(x1='wmean+wstd', x2='wmean-wstd', y='ys', source=sourceFin,
color='green', alpha=0.2, legend_label='WStd')
hover = HoverTool(names=['hist'],
tooltips=[('Interval', '@interval'),
('Count', str("@" + 'column'))])
plotHist.add_tools(hover)
plotHist.legend.location = "top_left"
plotHist.legend.click_policy="hide"
return plotHist
def update_hist_plot(attr, old, new):
try:
index = new[0]
flag = True
except IndexError:
sourceHist.data.update([])
flag = False
plotHist.visible = False
if flag:
plotHist.visible = True
plotHist.title.text = str(index)
gm = gmSelect.value
comp = compSelect.value
target = targetSelect.value
srcMain, srcSide = get_dataset(gdf, gm, comp, target)
srcHist = get_source_side(srcSide, index)
sourceHist.data.update(srcHist.data)
global gidx
gidx = index
######
args = parser_func()
filename = 'predicted_data.geojson'
file = args.datadir[0] + '/' + filename
m_coolwarm_rgb = (255 * cm.Reds(range(256))).astype('int')
coolwarm_palette = [RGB(*tuple(rgb)).to_hex() for rgb in m_coolwarm_rgb]
gm = 'pgv'
comp = 'Z'
target = 'mean'
global gidx
gidx = 0
gms = {'pga': 'pga',
'pgv': 'pgv'}
comps = {'Z': 'Z',
'N': 'N',
'E': 'E',
'H': 'H'}
targets = {'mean': 'mean',
'wmean': 'wmean',
'std': 'std',
'wstd': 'wstd'}
gmSelect = Select(value=gm, title='GM',
options=list(gms.keys()))
compSelect = Select(value=comp, title='Component',
options=list(comps.keys()))
targetSelect = Select(value=target, title='Targets',
options=list(targets.keys()))
gdf = gpd.read_file(file)
color_mapper = LinearColorMapper(palette='Viridis256', low=-1, high=2.3)
color_bar = ColorBar(color_mapper=color_mapper,
title='GM in log10', orientation='horizontal',
label_standoff=6, border_line_color=None, location='bottom_center')
sourceMain, sourceSide = get_dataset(gdf, gm, comp, target)
plot = make_main_plot(sourceMain, '%s_%s_%s' % (comp, gm, target))
sourceHist = get_source_side(sourceSide, gidx)
plotHist = make_hist_plot(sourceHist, '', gidx)
plotHist.visible = False
gmSelect.on_change('value', update_plot)
compSelect.on_change('value', update_plot)
targetSelect.on_change('value', update_plot)
sourceMain.selected.on_change('indices', update_hist_plot)
controls = column(gmSelect, compSelect, targetSelect)
curdoc().add_root(row(plot, controls, plotHist))
curdoc().title = "Interactive Ground-Motion Map"

10
src/config/gm_mapping.cfg

@ -1,13 +1,13 @@
-d /home/lehmann/dr/yaml_test/Norcia/DC_nearest_neighbor/
--method Pyrocko
--method NN
--mp 10
-s Z E N
-s Z
--gm pga pgv
--gf-database /home/lehmann/dr/pyrocko/own_2hz/
--H2 True
--delete-comp True
--nn-method sklearn
-m 3 1.5 1.5 rectangular
--nn-method tensorflow
-m 30 1.5 1.5 rectangular

0
src/GM/__init__.py → src/gm/__init__.py

48
src/GM/misc.py → src/gm/misc.py

@ -1,6 +1,7 @@
import random
import argparse
import sys
import h5py
import numpy as num
@ -131,6 +132,52 @@ def create_locationDict(fileLocation):
return locDict
def get_vs30(file, coords):
coords = num.array(coords)
f = h5py.File(file, 'r')
dCord = 0.1
flonMask = (min(coords.T[0]) - dCord < f['lon']) &\
(f['lon'] < max(coords.T[0]) + dCord)
flatMask = (min(coords.T[1]) - dCord < f['lat']) &\
(f['lat'] < max(coords.T[1]) + dCord)
flonIdx = num.where(flonMask)[0]
flatIdx = num.where(flatMask)[0]
flon = f['lon'][flonMask]
flat = f['lat'][flatMask]
fz = f['z'][min(flatIdx):max(flatIdx) + 1]
fnew = num.zeros((len(fz), len(flonIdx)))
for ii in range(len(fz)):
fnew[ii] = fz[ii][flonMask]
fz = fnew
dLon = flon[1] - flon[0]
dLat = flat[1] - flat[0]
vals = []
for refLon, refLat in coords:
for lon in flon:
if dLon > abs(refLon - lon):
idxLon = (num.where(flon == lon)[0][0])
break
for lat in flat:
if dLat > abs(refLat - lat):
idxLat = (num.where(flat == lat)[0][0])
break
vals.append(fz[idxLat][idxLon])
f.close()
del fz
del flon
del flat
return vals
#############################
### Waveform functions
#############################
@ -207,6 +254,7 @@ def create_synthetic_waveform(gfPath, source, coords, stf=None,
pyrockoSource = gf.DCSource(lat=source.lat, lon=source.lon,
depth=source.depth * 1000.,
magnitude=source.magnitude,
time=source.time,
dip=source.dip, strike=source.strike,
rake=source.rake)

114
src/GM/networks.py → src/gm/networks.py

@ -13,17 +13,19 @@ from pyrocko import orthodrome
from openquake.hazardlib.geo import geodetic
import ewrica.GM.sources as GMs
import ewrica.GM.plot as GMp
import ewrica.gm.sources as GMs
import ewrica.gm.plot as GMp
import matplotlib.backends.backend_pdf
from tensorflow import keras
import tensorflow as tf
#####
### General
#####
def get_gdf_NN(ii, args, src, lons, lats, srcsDict):
print('Starting: %s' % (ii))
@ -76,7 +78,6 @@ def get_gdf_NN(ii, args, src, lons, lats, srcsDict):
## rearranging the columns, otherwise the NN does not work
data = data[scalingDict.keys()]
### fix it multiple gms are the output
if args.nnmethod[0] == 'tensorflow':
test_dataset = tf.data.Dataset.from_tensor_slices((data.values))
data = test_dataset.batch(len(data))
@ -86,10 +87,10 @@ def get_gdf_NN(ii, args, src, lons, lats, srcsDict):
staDict = {}
for chagm, pr in zip(targets, pred):
staDict[chagm] = 10 ** pr
staDict[chagm] = pr
gdf = gpd.GeoDataFrame(staDict, geometry=gpd.points_from_xy(lons, lats))
srcsDict[src] = gdf
srcsDict[ii] = gdf
print('Finished: %s in %s' % (ii, time.time() - absRefTime))
@ -361,10 +362,8 @@ def NN_saving(hiddenLayers, args, xTrain, yTrain, scalingDict, targets):
return
def NN_testing(hiddenLayers, args,
xTrain, yTrain, xTest, yTest, evlData, targets, appendix):
evalData = evlData.copy(deep=True)
def NN_testing(hiddenLayers, args, xTrain, yTrain, xTest, yTest,
evaldata, scalingDict, targets, appendix):
hiddenLayers = ast.literal_eval(hiddenLayers)
if args.module[0] == 'sklearn':
@ -373,49 +372,52 @@ def NN_testing(hiddenLayers, args,
appendix += '_%s_%s' % (args.sks[0], args.ska[0])
pred = num.array(pred).T
#############################
### Evaluation
#############################
diffs = []
evalData = pd.DataFrame(evalData)
for target, pr in zip(targets, pred):
evalData[target] = 10 ** evalData[target]
evalData['pred%s' % (target)] = 10 ** pr
diff = num.log10(evalData[target] / evalData['pred%s' % (target)])
diffs.append(diff)
print('Finished: %s_%s_%s with mean:%0.3f & std:%0.3f' % (hiddenLayers,
args.sks[0], args.ska[0],
num.mean(num.mean(diffs)),
num.mean(num.std(diff))))
# output = '%0.3f, %0.3f | %s | %s' \
# % (num.mean(diff), num.std(diff), appendix, hiddenLayers)
# print(output, file=open('%s/results_%s.txt' % (args.outputfile,
# args.module[0]), "a"))
history = []
elif args.module[0] == 'tensorflow':
model, test_dataset, history = tensorflow_fit(args, hiddenLayers,
xTrain, yTrain, xTest, yTest)
pred = model.predict(test_dataset)
appendix += '_%s_%s_%s' % (args.ska[0], args.batchsize[0], args.epochs[0])
appendix += '_%s_%s_%s' % (args.ska[0], args.batchsize[0],
args.epochs[0])
pred = num.array(pred).T
#############################
### Evaluation
#############################
evalData = pd.DataFrame(evalData)
for target, pr in zip(targets, pred):
evalData[target] = 10 ** evalData[target]
evalData['pred%s' % (target)] = 10 ** pr
NN_evaluation(args, pred, targets, hiddenLayers, evaldata, scalingDict,
appendix, history=history)
return
def NN_evaluation(args, pred, targets, hiddenLayers, evaldata, scalingDict,
appendix, history=[]):
evlData = evaldata.copy(deep=True)
evlData = pd.DataFrame(evlData)
for target, pr in zip(targets, pred):
evlData['pred%s' % (target)] = pr
diff = evlData[target] - evlData['pred%s' % (target)]
print('Finished: %s_%s_%s %s with mean:%0.3f & std:%0.3f' % (hiddenLayers,
args.sks[0], args.ska[0],
target,
num.mean(diff),
num.std(diff)))
if args.module[0] == 'sklearn':
pass
# output = '%0.3f, %0.3f | %s | %s' \
# % (num.mean(diff), num.std(diff), appendix, hiddenLayers)
# print(output, file=open('%s/results_%s.txt' % (args.outputfile,
# args.module[0]), "a"))
if args.module[0] == 'tensorflow':
print('\nHistory:')
hist = pd.DataFrame(history.history)
hist['epoch'] = history.epoch
print(hist.tail(10))
print(hist.tail(11))
outputtrain = hist.tail(1)[['loss', 'mean_absolute_error',
'mean_squared_error']].values[0]
@ -427,24 +429,22 @@ def NN_testing(hiddenLayers, args,
outputval[2])
output = '%s | %s | %s %s | %s' \
% (outputtrain, outputval, appendix, hist.tail(1)[['epoch']].values[0], hiddenLayers)
% (outputtrain, outputval, appendix,
hist.tail(1)[['epoch']].values[0], hiddenLayers)
print(output, file=open('%s/results_%s.txt' % (args.outputfile,
args.module[0]), "a"))
#############################
### Plotting
#############################
# random.seed(2)
# evtNames = random.sample(list(set(evalData['eventname'])), 2)
cond = 'magnitude'
mini = evalData[evalData[cond] == min(evalData[cond])]
maxi = evalData[evalData[cond] == max(evalData[cond])]
mini = evlData[evlData[cond] == min(evlData[cond])]
maxi = evlData[evlData[cond] == max(evlData[cond])]
evtNames = [list(set(mini['eventname']))[0], list(set(maxi['eventname']))[0]]
print(evtNames)
# evtNames = ['Synthetic_1053', 'Synthetic_1127']
for evtName in evtNames:
evtData = evalData[evalData['eventname'] == evtName]
evtData = evlData[evlData['eventname'] == evtName]
source = GMs.SourceClass(
lat=float(evtData.hypolat.iloc[0]),
@ -479,11 +479,6 @@ def NN_testing(hiddenLayers, args,
GMp.plot_1d_gdf(gdf, refgdf, cols + targets, source, mode='dist')
GMp.plot_1d_gdf(gdf, refgdf, cols + targets, source, mode='azimuth')
# GMp.plot_1d(obsCont=obsCont, resCont=predCont, mode='distance',
# distType='hypo', savename=[])
# GMp.plot_1d(obsCont=obsCont, resCont=predCont, mode='azimuth',
# savename=[])
pdf = matplotlib.backends.backend_pdf.PdfPages('%s/%s_%s_%s.pdf'
% (args.outputfile, evtName, hiddenLayers, appendix))
for fig in range(1, plt.gcf().number + 1):
@ -492,15 +487,13 @@ def NN_testing(hiddenLayers, args,
plt.close('all')
print('Finished Plotting %s' % (evtName))
return
### sklearn
def sklearn_fit(args, layers, xTrain, yTrain):
from sklearn.neural_network import MLPRegressor
print('Starting:', layers)
model = MLPRegressor(solver=args.sks[0], alpha=1e-5,
model = MLPRegressor(solver=args.sks[0], alpha=1e-5, verbose=True,
activation=args.ska[0], max_iter=args.skn[0],
hidden_layer_sizes=layers)
model.fit(xTrain, yTrain)
@ -620,6 +613,11 @@ class EarlyStoppingAtMinLoss(keras.callbacks.Callback):
def adapt_learning_rate(epoch):
# ‘adaptive’ keeps the learning rate constant to ‘learning_rate_init’
# as long as training loss keeps decreasing. Each time two consecutive
# epochs fail to decrease training loss by at least tol, or fail to increase
# validation score by at least tol if ‘early_stopping’ is on, the current
# learning rate is divided by 5.
return 0.05 / ((epoch * 2) + 1)
@ -632,9 +630,9 @@ def tensorflow_fit(args, layers, xTrain, yTrain, xTest=[], yTest=[]):
yTrain.values))
train_dataset = train_dataset.shuffle(len(xTrain)).batch(args.batchsize[0])
my_lr_scheduler = keras.callbacks.LearningRateScheduler(adapt_learning_rate)
callbacks = [my_lr_scheduler]
# callbacks = []
# my_lr_scheduler = keras.callbacks.LearningRateScheduler(adapt_learning_rate)
# callbacks = [my_lr_scheduler]
callbacks = []
if len(xTest) == 0 and len(yTest) == 0:
callbacks.append(EarlyStoppingAtMinLoss(10))
@ -652,7 +650,7 @@ def tensorflow_fit(args, layers, xTrain, yTrain, xTest=[], yTest=[]):
test_dataset = tf.data.Dataset.from_tensor_slices((xTest.values,
yTest.values))
test_dataset = test_dataset.batch(args.batchsize[0])
callbacks.append(EarlyStoppingAtMinLoss_validate(5))
callbacks.append(EarlyStoppingAtMinLoss_validate(10))
history = model.fit(train_dataset,
shuffle=True,
verbose=2,

301
src/GM/plot.py → src/gm/plot.py

@ -9,7 +9,7 @@ from matplotlib.colors import BoundaryNorm
from mpl_toolkits import basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
import ewrica.GM.misc as GMm
import ewrica.gm.misc as GMm
#############################
### Plot functions
@ -436,148 +436,6 @@ def plot_1d(obsCont, resCont, mode='dist', distType='hypo',
plt.savefig('%s.png' % (savename))
def plot_1d_gdf(gdf, refgdf, chagms, source, add='', mode='dist', distType='hypo',
savename='gm_diagram'):
chas, ims = num.array([ii.rsplit('_') for ii in chagms]).T
ims = list(set(ims))
chas = list(set(chas))
if mode in ['dist', 'distance']:
if hasattr(source, 'surface'):
if distType == 'hypo':
distStr = 'Hypocentral'
elif distType == 'rrup':
distStr = 'Minimal Rupture'
elif distType == 'rx':
distStr = 'Perpendicular-to-Strike (rx)'
elif distType == 'rjb':
distStr = 'Joyner-Boore'
elif distType == 'ry0':
distStr = 'Parallel-to-Strike (ry0)'
else:
distType = 'hypo'
distStr = 'Hypocentral'
valMax = max(max(gdf[chagms].max()), max(refgdf[chagms].max()))
valMin = min(min(gdf[chagms].min()), min(refgdf[chagms].min()))
outCol = len(chas)
outRow = len(ims)
figWidth = outCol * 10
figLength = outRow * 10
fig = plt.figure(figsize=(figWidth, figLength))
outer = gridspec.GridSpec(outRow, outCol, wspace=0., hspace=0.)
compCnt = -1
for im in ims:
compCnt += 1
n = -1
for cha in chas:
n += 1
### data
lons = list(gdf.geometry.x)
lats = list(gdf.geometry.y)
reflons = list(refgdf.geometry.x)
reflats = list(refgdf.geometry.y)
if add == '':
data = num.array(gdf['%s_%s' % (cha, im)])
refdata = num.array(refgdf['%s_%s' % (cha, im)])
else:
data = num.array(gdf['%s_%s_%s' % (cha, im, add)])
refdata = num.array(refgdf['%s_%s_%s' % (cha, im, add)])
data = num.log10(data)
refdata = num.log10(refdata)
if type(data) == float:
continue
inner = gridspec.GridSpecFromSubplotSpec(2, 1,
subplot_spec=outer[compCnt * len(chas) + (n)],
wspace=0., hspace=0.,
width_ratios=[1], height_ratios=[3, 1])
ax1 = plt.Subplot(fig, inner[0])
ax2 = plt.Subplot(fig, inner[1])
ax2color = 'b'
if mode in ['dist', 'distance']:
'''
Distance
'''
ydata = GMm.get_distances(lons, lats, source, distType=distType)
ydataref = GMm.get_distances(reflons, reflats, source, distType=distType)
if compCnt + 1 == len(ims):
ax2.set_xlabel('%s - Distance [km]' % (distStr))
elif mode in ['azi', 'azimuth']:
'''
Azimuth
'''
ydata = GMm.get_azimuths(lons, lats, source)
ydataref = GMm.get_azimuths(reflons, reflats, source)
ax1.set_xlim((-180., 180))
ax2.set_xlim((-180., 180))
if compCnt + 1 == len(ims):
ax2.set_xlabel('Azimuth [Degree]')
else:
print('Wrong mode')
exit()
if n == 0:
ax1.set_ylabel('Amplitude [log]')
ax2.set_ylabel('Residuum [log]', color=ax2color)
'''
Plotting
'''
## Primary plot, e.g. Azimuth or distance
ax1.plot(ydata, data, '*r', label='GM_Pre')
ax1.plot(ydataref, refdata, 'og', label='GM_Obs', fillstyle='none')
ax1.set_ylim((num.log10(valMin / 2.), num.log10(valMax * 5.)))
## Secondary plot, residual
residuum = refdata - data
ax2.plot(ydata, residuum, '+', color=ax2color, label='log10-Res')
coef = num.polyfit(ydata, residuum, 1)
poly1d_fn = num.poly1d(coef)
dummyX = num.linspace(min(ydata), max(ydata))
ax2.plot(dummyX, poly1d_fn(dummyX), '--', color=ax2color)
ax2.set_ylim((min(-1.1, min(residuum)), max(1.1, max(residuum))))
ax2.axhline(y=0, color=ax2color, linestyle='-', alpha=0.25)
ax2.axhline(y=1, color=ax2color, linestyle=':', alpha=0.5)
ax2.axhline(y=-1, color=ax2color, linestyle=':', alpha=0.5)
ax2.tick_params(axis='y', colors=ax2color)
# ax1.set_title('%s\n%s vs. %s\nres: %0.2f, std: %0.2f, slope:%0.1E'
# % (gm, comp, obsComp, num.nanmean(residuum),
# num.nanstd(residuum), coef[0]))
fig.add_subplot(ax1)
fig.add_subplot(ax2)
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines + lines2, labels + labels2)
if mode in ['dist', 'distance']:
plt.suptitle('Distance')
else:
plt.suptitle('Azimuth')
outer.tight_layout(fig, rect=[0., 0., 1., 0.98])
if savename != []:
plt.savefig('%s.png' % (savename))
def plot_gdf(gdf, chagms, refgdf=False, source=False, mode='default', add='',
savename='gm_map',
cmapname='viridis', predPlotMode='scatter', smoothfac=False):
@ -617,9 +475,8 @@ def plot_gdf(gdf, chagms, refgdf=False, source=False, mode='default', add='',
else:
refdata = num.array(refgdf['%s_%s_%s' % (cha, im, add)])
data = num.log10(refdata) - num.log10(data)
data = refdata - data
shmLevels = num.linspace(-1, 1, 21)
mode = 'diff'
elif add in ['std', 'stds', 'wstd', 'wstds']:
shmLevels = num.array(num.linspace(-2, 2, 21))
@ -630,9 +487,6 @@ def plot_gdf(gdf, chagms, refgdf=False, source=False, mode='default', add='',
10, 20, 50,
100, 200]))
if mode == 'default':
data = num.log10(data)
if type(data) == float:
continue
@ -749,13 +603,162 @@ def plot_gdf(gdf, chagms, refgdf=False, source=False, mode='default', add='',
ticks=shmLevels,
label=label)
titlestr = '%s\n' % (add)
ax.set_title('%s\n%s' % (im, cha))
fig.add_subplot(ax)
if source is not False:
titlestr = '%s\nMag: %0.1f, Depth: %0.1f' \
% (source.name, source.magnitude, source.depth)
# titlestr += ' Strike: %0.1f, Dip: %0.1f, Rake: %0.1f' \
# % (source.strike, source.dip, source.rake)
# titlestr += '\nNucX: %.1f, NucY: %.1f' \
# % (source.nucleation_X, source.nucleation_Y)
titlestr += '%s\n' % (add)
else:
titlestr = '%s\n' % (add)
plt.suptitle(titlestr)
outer.tight_layout(fig, rect=[0., 0., 1., 0.95])
if savename != []:
fig.savefig('%s.png' % (savename))
def plot_1d_gdf(gdf, refgdf, chagms, source, add='', mode='dist', distType='hypo',
savename='gm_diagram'):
chas, ims = num.array([ii.rsplit('_') for ii in chagms]).T
ims = list(set(ims))
chas = list(set(chas))
if mode in ['dist', 'distance']:
if hasattr(source, 'surface'):
if distType == 'hypo':
distStr = 'Hypocentral'
elif distType == 'rrup':
distStr = 'Minimal Rupture'
elif distType == 'rx':
distStr = 'Perpendicular-to-Strike (rx)'
elif distType == 'rjb':
distStr = 'Joyner-Boore'
elif distType == 'ry0':
distStr = 'Parallel-to-Strike (ry0)'
else:
distType = 'hypo'
distStr = 'Hypocentral'
valMax = max(max(gdf[chagms].max()), max(refgdf[chagms].max()))
valMin = min(min(gdf[chagms].min()), min(refgdf[chagms].min()))
outCol = len(chas)
outRow = len(ims)
figWidth = outCol * 10
figLength = outRow * 10
fig = plt.figure(figsize=(figWidth, figLength))
outer = gridspec.GridSpec(outRow, outCol, wspace=0., hspace=0.)
compCnt = -1
for im in ims:
compCnt += 1
n = -1
for cha in chas:
n += 1
### data
lons = list(gdf.geometry.x)
lats = list(gdf.geometry.y)
reflons = list(refgdf.geometry.x)
reflats = list(refgdf.geometry.y)
if add == '':
data = num.array(gdf['%s_%s' % (cha, im)])
refdata = num.array(refgdf['%s_%s' % (cha, im)])
else:
data = num.array(gdf['%s_%s_%s' % (cha, im, add)])
refdata = num.array(refgdf['%s_%s_%s' % (cha, im, add)])
if type(data) == float:
continue
inner = gridspec.GridSpecFromSubplotSpec(2, 1,
subplot_spec=outer[compCnt * len(chas) + (n)],
wspace=0., hspace=0.,
width_ratios=[1], height_ratios=[3, 1])
ax1 = plt.Subplot(fig, inner[0])
ax2 = plt.Subplot(fig, inner[1])
ax2color = 'b'
if mode in ['dist', 'distance']:
'''
Distance
'''
ydata = GMm.get_distances(lons, lats, source, distType=distType)
ydataref = GMm.get_distances(reflons, reflats, source, distType=distType)
if compCnt + 1 == len(ims):
ax2.set_xlabel('%s - Distance [km]' % (distStr))
elif mode in ['azi', 'azimuth']:
'''
Azimuth
'''
ydata = GMm.get_azimuths(lons, lats, source)
ydataref = GMm.get_azimuths(reflons, reflats, source)
ax1.set_xlim((-180., 180))
ax2.set_xlim((-180., 180))
if compCnt + 1 == len(ims):
ax2.set_xlabel('Azimuth [Degree]')
else:
print('Wrong mode')
exit()
if n == 0:
ax1.set_ylabel('Amplitude [log]')
ax2.set_ylabel('Residuum [log]', color=ax2color)
'''
Plotting
'''
## Primary plot, e.g. Azimuth or distance
ax1.plot(ydata, data, '*r', label='GM_Pre')
ax1.plot(ydataref, refdata, 'og', label='GM_Obs', fillstyle='none')
ax1.set_ylim((valMin - num.log10(2.), valMax + num.log10(5.)))
## Secondary plot, residual
residuum = refdata - data
ax2.plot(ydata, residuum, '+', color=ax2color, label='log10-Res')
coef = num.polyfit(ydata, residuum, 1)
poly1d_fn = num.poly1d(coef)
dummyX = num.linspace(min(ydata), max(ydata))
ax2.plot(dummyX, poly1d_fn(dummyX), '--', color=ax2color)
ax2.set_ylim((min(-1.1, min(residuum)), max(1.1, max(residuum))))
ax2.axhline(y=0, color=ax2color, linestyle='-', alpha=0.25)
ax2.axhline(y=1, color=ax2color, linestyle=':', alpha=0.5)
ax2.axhline(y=-1, color=ax2color, linestyle=':', alpha=0.5)
ax2.tick_params(axis='y', colors=ax2color)
# ax1.set_title('%s\n%s vs. %s\nres: %0.2f, std: %0.2f, slope:%0.1E'
# % (gm, comp, obsComp, num.nanmean(residuum),
# num.nanstd(residuum), coef[0]))
fig.add_subplot(ax1)
fig.add_subplot(ax2)
lines, labels = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines + lines2, labels + labels2)
if mode in ['dist', 'distance']:
plt.suptitle('Distance')
else:
plt.suptitle('Azimuth')
outer.tight_layout(fig, rect=[0., 0., 1., 0.98])
if savename != []:
plt.savefig('%s.png' % (savename))

41
src/GM/sources.py → src/gm/sources.py

@ -4,14 +4,14 @@ import geopandas as gpd
import geojson
from pyrocko import orthodrome, gf, trace, util
from pyrocko.guts import Object, StringChoice, Float, Int, \
String, List, Any, Dict, Choice
from pyrocko.guts import Object, StringChoice, Float, String, List, Dict, Choice
from openquake.hazardlib.geo import geodetic, Point, Mesh, PlanarSurface
import ewrica.sources as ews
import ewrica.GM.misc as GMm
from ewrica.GM.misc import G, deg
from ewrica.gm.misc import G, deg
import ewrica.gm.misc as GMm
#############################
### Sources, probably move them to ew.sources
@ -429,7 +429,6 @@ class StationContainer(Object, Cloneable):
return
def calc_azimuths(self):
lons = []
@ -445,10 +444,8 @@ class StationContainer(Object, Cloneable):
for nn, ii in enumerate(self.stations):
self.stations[ii].azimuth = float(azimuth[nn])
return
def calc_h2_from_gm(self, delete=True, imts=['pgv', 'pga'],
SDOFs=[0.3, 1, 3]):
'''
@ -473,7 +470,6 @@ class StationContainer(Object, Cloneable):
for gm in delGM:
del self.stations[sta].components[comp].gms[gm]
if 'E' in self.stations[sta].components.keys() and \
'N' in self.stations[sta].components.keys():
@ -498,7 +494,6 @@ class StationContainer(Object, Cloneable):
del self.stations[sta].components['N']
del self.stations[sta].components['E']
if '1' in self.stations[sta].components.keys() and \
'2' in self.stations[sta].components.keys():
@ -519,14 +514,12 @@ class StationContainer(Object, Cloneable):
COMP.gms[gm] = GM
self.stations[sta].components['H'] = COMP
if delete:
del self.stations[sta].components['1']
del self.stations[sta].components['2']
return
def get_gm_from_wv(self, imts=['pga', 'pgv'],
SDOFs=[0.3, 1.0, 3.0],
H2=False, delete=False, deleteWvData=True):
@ -539,7 +532,7 @@ class StationContainer(Object, Cloneable):
intV = trace.IntegrationResponse(1)
for sta in self.stations:
if H2 == True:
if H2:
cha = 'H'
flag = False
@ -587,7 +580,7 @@ class StationContainer(Object, Cloneable):
GM = GMClass(
name=gm,
value=float(val),
value=float(num.log10(val)),
unit=unit)
COMP.gms[GM.name] = GM
@ -630,7 +623,7 @@ class StationContainer(Object, Cloneable):
GM = GMClass(
name=gm,
value=float(val),
value=float(num.log10(val)),
unit=unit)
COMP.gms[GM.name] = GM
@ -659,6 +652,22 @@ class StationContainer(Object, Cloneable):
for sta in rmSta:
del self.stations[sta]
def remove_distant_stations_radius(self, refdist):
rmSta = []
source = self.refSource
for sta in self.stations:
lon = self.stations[sta].lon
lat = self.stations[sta].lat
dist = orthodrome.distance_accurate50m_numpy(source.lat, source.lon,
lat, lon)[0]
if dist > refdist:
rmSta.append(sta)
for sta in rmSta:
del self.stations[sta]
def remove_stations_without_components(self, componentList):
rmSta = []
for sta in self.stations:
@ -904,8 +913,8 @@ def get_gdf_Pyrocko(ii, args, src, engine, waveform_targets, srcsDict):
staCont.get_gm_from_wv(imts=args.imts, SDOFs=args.sdofs,
H2=args.H2, delete=args.delete, deleteWvData=True)
staCont.validate()
# staCont.validate()
gdf = staCont.to_geodataframe()
srcsDict[src] = gdf
srcsDict[ii] = gdf
print('Finished: %s' % (ii))

BIN
src/mapping_functions/NN/DCSource_sklearn.bin

Binary file not shown.

BIN
src/mapping_functions/NN/DCSource_sklearn_model.bin

Binary file not shown.

BIN
src/mapping_functions/NN/DCSource_tensorflow.bin

Binary file not shown.

BIN
src/mapping_functions/NN/DCSource_tensorflow_model.h5

Binary file not shown.
Loading…
Cancel
Save