-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathshakemap_sampler_test.py
71 lines (63 loc) · 2.9 KB
/
shakemap_sampler_test.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 28 10:01:50 2022
@author: Hugo Rosero
"""
#import shakemap
#import shakeml
import quakeml
import numpy as np
import random
import os
import lxml.etree as le
import pandas
import io
import datetime
import sampler
from openquake.hazardlib.imt import PGA, PGV, IA, SA
import correlation as crl
from openquake.hazardlib.site import Site, SiteCollection
from openquake.hazardlib.geo import Mesh, Point
import pylab
import scipy.stats as stats
import scipy.linalg as spla
"""
Reads a normal shakemap (as it is the output of shakyground for earthquake scenarios)
:return: The shakemap with the random residuals (separated and mixed)
"""
#local testing
#change this directory where the test file shakemap.xml is located
file_name="C:\\Users\\Public\\JournalPaperGFZ\\PythonCodesBD\\deus-master\\testinputs\\shakemap.xml"
shakemap_outfile="C:\\Users\\Public\\JournalPaperGFZ\\PythonCodesBD\\shakyground-master\\out_shakemap_uncorr.xml"
shakemap_outfile_corr="C:\\Users\\Public\\JournalPaperGFZ\\PythonCodesBD\\shakyground-master\\out_shakemap_corr.xml"
random_seed=123
event,units,grid_data, event_specific_uncertainties, regular_grid = sampler.extract_shakemap_data(file_name)
grid_data,units = sampler.create_uncorrelated_residuals(grid_data,units,random_seed)
sampler.save_random_shakemap(shakemap_outfile,event,units,grid_data, event_specific_uncertainties,regular_grid,random_seed)
#L=sampler.build_sparse_correlation_matrix(grid_data,PGA,tol=1e-2)
grid_data,units = sampler.create_correlated_residuals(grid_data,units,random_seed)
sampler.save_random_shakemap(shakemap_outfile_corr,event,units,grid_data, event_specific_uncertainties,regular_grid,random_seed)
#with open("C:\\Users\\Public\\JournalPaperGFZ\\PythonCodesBD\\shakyground-master\\out_shakemap.xml", 'w') as f:
#f.write(eq_with_residuals)
#test correlation model
print('test')
#rand_seed=1234578
rand_idx=np.random.permutation(grid_data.index)
npts=500
first_idx=[rand_idx[i] for i in range(npts)]
gx=grid_data['LAT'][first_idx]
gy=grid_data['LON'][first_idx]
rpga=grid_data['RESPGA'][first_idx]#obtain the residuals
sites = SiteCollection([Site(Point(x,y), 1, vs30measured=True, z1pt0=3.4, z2pt5=5.6, backarc=False) for x, y in zip(gx, gy)])#dummy vs30. not needed in this test
corrmat=crl.jbcorrelation(sites, PGA, True)
print('corr mat done')
Lcorrmat=spla.sqrtm(corrmat)#np.linalg.cholesky(corrmat)
print('cholesky done')
Lcorrmatinv=np.linalg.inv(Lcorrmat)
print('inverse done')
Lcorrmatinv=[[Lcorrmatinv[i,j] for j in range(npts)] for i in range(npts)]
print('rebuild inverse done')
std_residuals=np.dot(Lcorrmatinv,list(rpga.array))#these should be independent, standard normal distributed random variables
#print('mean: '+str(np.mean(std_residuals))+' std dev: '+str(np.std(std_residuals)))
stats.probplot(std_residuals, dist="norm", plot=pylab)
pylab.show()