Skip to content

Commit cbc3dbe

Browse files
Diego MelgarDiego Melgar
Diego Melgar
authored and
Diego Melgar
committed
Updatesd
1 parent e3ecb74 commit cbc3dbe

4 files changed

Lines changed: 181 additions & 44 deletions

File tree

‎scripts/quadtree_los_ixtepec.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,16 @@
11
format compact
22

3-
fout='/Users/dmelgar/Ixtepec2017/InSAR/T107_Ascending/ascending_quadtree.los'
4-
insar=textread('/Users/dmelgar/Ixtepec2017/InSAR/T107_Ascending/los.lltnde');
3+
fout='/Users/dmelgar/Ixtepec2017/InSAR/T107_Ascending/asc_quadtree.los'
4+
insar=textread('/Users/dmelgarm/Turkey2020/InSAR/asc_los_detrended.txt');
55
mindim=32
66
maxdim=256
77
thresh=0.04
88

9-
plot_lims=[-0.1,0.1]
9+
plot_lims=[-0.3,0.3]
1010

1111
lon=insar(:,1);
1212
lat=insar(:,2);
13-
los=insar(:,7)/1000;
13+
los=insar(:,3);
1414
lookE=insar(:,4);
1515
lookN=insar(:,5);
1616
lookU=insar(:,6);

‎scripts/quadtree_los_turkey.m‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
format compact
2+
3+
fout='/Users/dmelgarm/Turkey2020/InSAR/desc_quadtree.los'
4+
insar=textread('/Users/dmelgarm/Turkey2020/InSAR/desc_los_detrended.txt');
5+
mindim=32
6+
maxdim=256
7+
thresh=0.1
8+
9+
plot_lims=[-0.3,0.3]
10+
11+
lon=insar(:,1);
12+
lat=insar(:,2);
13+
los=insar(:,3);
14+
lookE=insar(:,4);
15+
lookN=insar(:,5);
16+
lookU=insar(:,6);
17+
18+
%Prepare grid
19+
dlon=8.34e-4;
20+
dlat=8.34e-4;
21+
22+
lon_i=linspace(min(lon),max(lon),2048);
23+
lat_i=linspace(min(lat),max(lat),2048);
24+
25+
26+
[X,Y]=meshgrid(lon_i,lat_i);
27+
los_interp_sharp = griddata(lon,lat,los,X,Y);
28+
29+
los_interp=los_interp_sharp;
30+
31+
display('QT decomp')
32+
s=qtdecomp(los_interp,thresh,[mindim,maxdim]);
33+
%Calcualte maximum dimension
34+
maxdim=max(max(full(s)));
35+
%get power of 2 of possible block values
36+
idim=(log(mindim)/log(2)):1:(log(maxdim)/log(2));
37+
%Now loop through
38+
los_out=[];
39+
lon_out=[];
40+
lat_out=[];
41+
for k=1:length(idim)
42+
[vals, r, c] = qtgetblk(los_interp_sharp, s,2^idim(k));
43+
icurrent=find(s==2^idim(k));
44+
%Now get mean and cellc enter of each grid
45+
for kgrid=1:length(icurrent)
46+
%Get values
47+
new_los=mean(mean(vals(:,:,kgrid)));
48+
if ~isnan(new_los) %Add to the list
49+
los_out=[los_out,mean(mean(vals(:,:,kgrid)))];
50+
%get indices of center as upepr left plus half the size of the
51+
%block
52+
r1=r(kgrid)+(2^idim(k))/2;
53+
r2=r(kgrid)+(2^idim(k))/2+1;
54+
c1=c(kgrid)+(2^idim(k))/2;
55+
c2=c(kgrid)+(2^idim(k))/2+1;
56+
%Now figure out coordiantes of the center
57+
lon_out=[lon_out,0.5*(X(r1,c1)+X(r2,c2))];
58+
lat_out=[lat_out,0.5*(Y(r1,c1)+Y(r2,c2))];
59+
end
60+
end
61+
end
62+
63+
64+
65+
%Find closest look direction vector
66+
lookE_out=[];
67+
lookN_out=[];
68+
lookU_out=[];
69+
for k=1:length(lon_out)
70+
d=sqrt((lon_out(k)-lon).^2+(lat_out(k)-lat).^2);
71+
[a,i]=min(d);
72+
lookE_out=[lookE_out,lookE(i)];
73+
lookN_out=[lookN_out,lookN(i)];
74+
lookU_out=[lookU_out,lookU(i)];
75+
end
76+
77+
%other random filters
78+
i=find(lon_out>38.68);
79+
lon_out=lon_out(i);
80+
lat_out=lat_out(i);
81+
los_out=los_out(i);
82+
lookE_out=lookE_out(i);
83+
lookN_out=lookN_out(i);
84+
lookU_out=lookU_out(i);
85+
i=find(lat_out>38.11);
86+
lon_out=lon_out(i);
87+
lat_out=lat_out(i);
88+
los_out=los_out(i);
89+
lookE_out=lookE_out(i);
90+
lookN_out=lookN_out(i);
91+
lookU_out=lookU_out(i);
92+
93+
%Write to file
94+
out=[lon_out' lat_out' los_out' lookE_out' lookN_out' lookU_out'];
95+
dlmwrite(fout,out,'delimiter','\t','precision',6)
96+
97+
%Plot
98+
figure
99+
subplot(1,2,1)
100+
mesh(X,Y,los_interp_sharp)
101+
colorbar
102+
view([0,90])
103+
hold on
104+
scatter3(lon_out,lat_out,100000*ones(size(lat_out)),'kx')
105+
xlabel('Longitude','FontSize',16)
106+
ylabel('Latitude','FontSize',16)
107+
108+
title('LOS (m)','FontSize',16)
109+
caxis(plot_lims)
110+
colormap('jet')
111+
set(gca,'FontSize',14)
112+
subplot(1,2,2)
113+
scatter(lon_out,lat_out,40,los_out,'filled')
114+
colorbar
115+
view([0,90])
116+
caxis(plot_lims)
117+
colormap('jet')
118+
title('QuadTree Resampled LOS (m)','FontSize',16)
119+
xlabel('Longitude','FontSize',16)
120+
ylabel('Latitude','FontSize',16)
121+
grid on
122+
set(gca,'FontSize',14)
123+
display([int2str(length(lat_out)) ' grids created'])
124+
125+

‎src/python/mudpy/insartools.py‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,7 @@
11
'''
2-
InSAR tools gop in here
2+
InSAR tools go in here
33
Diego Melgar
4-
UC Berkeley
5-
05/2015
4+
University of Oregon
65
'''
76

87
def quadtree2mudpy(home,project_name,quadtree_file,gflist_file,prefix):
@@ -31,13 +30,13 @@ def quadtree2mudpy(home,project_name,quadtree_file,gflist_file,prefix):
3130
E=insar[k,3]
3231
U=insar[k,5]
3332
los=c_[los_c,N,E,U]
34-
#Make plot as you go to verify
35-
plt.scatter(insar[k,0],insar[k,1],c=los_c,cmap=cm.magma,s=80,lw=0,vmin=vmin,vmax=vmax)
36-
plt.legend(['LOS'])
3733
savetxt(home+project_name+'/data/statics/'+sta+'.los',los,header='LOS(m),los unit vector (positive towards satellite) n,e,u')
3834
#Generate gflist file as well
3935
gflist.write('%s\t%10.6f\t%10.6f\t0\t0\t0\t0\t1\t/foo/bar\t/foo/bar\t/foo/bar\t/foo/bar\t%s\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n'%(sta,insar[k,0],insar[k,1],home+project_name+'/data/statics/'+sta+'.los'))
4036

37+
#Make plot to verify
38+
plt.scatter(insar[:,0],insar[:,1],c=insar[:,2],cmap=cm.magma,s=80,lw=0,vmin=vmin,vmax=vmax)
39+
plt.legend(['LOS'])
4140
plt.colorbar()
4241
gflist.close()
4342
plt.show()

‎src/python/mudpy/view.py‎

Lines changed: 47 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -479,7 +479,7 @@ def tile_slip(rupt,nstrike,ndip,slip_bounds,geographic=False,epicenter=0,epicent
479479

480480
from numpy import genfromtxt,unique,where,zeros,tile,sqrt
481481
import matplotlib.pyplot as plt
482-
from obspy.core.util.geodetics import gps2DistAzimuth
482+
from obspy.geodetics import gps2dist_azimuth
483483

484484
f=genfromtxt(rupt)
485485
num=f[:,0]
@@ -504,19 +504,19 @@ def tile_slip(rupt,nstrike,ndip,slip_bounds,geographic=False,epicenter=0,epicent
504504
slip_min=slip_bounds[0]
505505
slip_max=slip_bounds[1]
506506
#Aftershocks
507-
af=genfromtxt('/Users/dmelgar/Chiapas2017/afters/afters_10days_close2fault.txt')
508-
lon_afters=af[:,0]
509-
lat_afters=af[:,1]
510-
depth_afters=af[:,2]
511-
print(af)
507+
# af=genfromtxt('/Users/dmelgar/Chiapas2017/afters/afters_10days_close2fault.txt')
508+
# lon_afters=af[:,0]
509+
# lat_afters=af[:,1]
510+
# depth_afters=af[:,2]
511+
# print(af)
512512
if geographic==True: #Get geographic coordinates to compute along strike and along dip distance
513513
lon=f[(epicenter_line-1)*nstrike:epicenter_line*nstrike,1] #Only compute line at the epicenter depth
514514
lat=f[(epicenter_line-1)*nstrike:epicenter_line*nstrike,2]
515515
depth=-f[:,3]
516516
depth=depth[0:len(unum)]
517517
along_strike=zeros(nstrike)
518518
for k in range(len(lat)):
519-
out=gps2DistAzimuth(epicenter[1],epicenter[0],lat[k],lon[k])
519+
out=gps2dist_azimuth(epicenter[1],epicenter[0],lat[k],lon[k])
520520
if lat[k]<epicenter[1]: #It's to the south
521521
#along_strike[k]=-out[0]/1000
522522
along_strike[k]=out[0]/1000
@@ -525,15 +525,15 @@ def tile_slip(rupt,nstrike,ndip,slip_bounds,geographic=False,epicenter=0,epicent
525525
along_strike[k]=-out[0]/1000
526526
#Now tile
527527
along_strike=tile(along_strike,ndip)
528-
#Process the aftershocks
529-
along_strike_afters=zeros(len(lon_afters))
530-
for k in range(len(lat_afters)):
531-
out=gps2DistAzimuth(epicenter[1],epicenter[0],lat_afters[k],lon_afters[k])
532-
if lat_afters[k]<epicenter[1]: #It's to the south
533-
#along_strike_afters[k]=-out[0]/1000
534-
along_strike_afters[k]=out[0]/1000
535-
else:
536-
along_strike_afters[k]=-out[0]/1000
528+
# #Process the aftershocks
529+
# along_strike_afters=zeros(len(lon_afters))
530+
# for k in range(len(lat_afters)):
531+
# out=gps2dist_azimuth(epicenter[1],epicenter[0],lat_afters[k],lon_afters[k])
532+
# if lat_afters[k]<epicenter[1]: #It's to the south
533+
# #along_strike_afters[k]=-out[0]/1000
534+
# along_strike_afters[k]=out[0]/1000
535+
# else:
536+
# along_strike_afters[k]=-out[0]/1000
537537
#Get indices for plot
538538
istrike=zeros(nstrike*ndip)
539539
idip=zeros(nstrike*ndip)
@@ -1576,7 +1576,8 @@ def plot_data(home,project_name,gflist,vord,decimate,lowpass,t_lim,sort,scale,k_
15761576
#plt.subplots_adjust(left=0.2, bottom=0.05, right=0.8, top=0.95, wspace=0, hspace=0)
15771577
plt.subplots_adjust(left=0.2, bottom=0.15, right=0.8, top=0.85, wspace=0, hspace=0)
15781578

1579-
def synthetics(home,project_name,run_name,run_number,gflist,vord,decimate,lowpass,t_lim,sort,scale,k_or_g,uncert=False,waveforms_as_accel=False):
1579+
def synthetics(home,project_name,run_name,run_number,gflist,vord,decimate,lowpass,t_lim,
1580+
sort,scale,k_or_g,uncert=False,waveforms_as_accel=False,units='m',uncerth=0.01,uncertv=0.03):
15801581
'''
15811582
Plot synthetics vs real data
15821583
@@ -1630,12 +1631,14 @@ def synthetics(home,project_name,run_name,run_number,gflist,vord,decimate,lowpas
16301631
es=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.e.sac')
16311632
us=read(synthpath+run_name+'.'+run_number+'.'+sta[i[k]]+'.'+synthsuffix+'.u.sac')
16321633

1633-
#Zero out a stations
1634-
if sta[i[k]]=="NILT":
1635-
n[0].data=zeros(len(n[0].data))
1636-
e[0].data=zeros(len(e[0].data))
1637-
ns[0].data=zeros(len(ns[0].data))
1638-
es[0].data=zeros(len(es[0].data))
1634+
if units=='cm':
1635+
n[0].data*=100
1636+
e[0].data*=100
1637+
u[0].data*=100
1638+
ns[0].data*=100
1639+
es[0].data*=100
1640+
us[0].data*=100
1641+
16391642

16401643
if lowpass!=None:
16411644
print('Lowpassing')
@@ -1674,9 +1677,9 @@ def synthetics(home,project_name,run_name,run_number,gflist,vord,decimate,lowpas
16741677
axe.grid(which='both')
16751678
axu.plot(u[0].times(),u[0].data,'k',us[0].times(),us[0].data,'r')
16761679
if uncert==True:
1677-
axn.fill_between(n[0].times(),n[0].data-0.01,n[0].data+0.01,alpha=0.2,color='k')
1678-
axe.fill_between(e[0].times(),e[0].data-0.01,e[0].data+0.01,alpha=0.2,color='k')
1679-
axu.fill_between(u[0].times(),u[0].data-0.03,u[0].data+0.03,alpha=0.2,color='k')
1680+
axn.fill_between(n[0].times(),n[0].data-uncerth,n[0].data+uncerth,alpha=0.2,color='k')
1681+
axe.fill_between(e[0].times(),e[0].data-uncerth,e[0].data+uncerth,alpha=0.2,color='k')
1682+
axu.fill_between(u[0].times(),u[0].data-uncertv,u[0].data+uncertv,alpha=0.2,color='k')
16801683
axu.grid(which='both')
16811684
axe.yaxis.set_ticklabels([])
16821685
axu.yaxis.set_ticklabels([])
@@ -1746,13 +1749,23 @@ def synthetics(home,project_name,run_name,run_number,gflist,vord,decimate,lowpas
17461749
axn.set_ylabel(sta[i[k]],rotation=0,labelpad=20)
17471750
if k==0:
17481751
if vord.lower()=='d':
1749-
axn.set_title('North (m)')
1750-
axe.set_title('East (m)')
1751-
axu.set_title('Up (m)')
1752+
if units=='cm':
1753+
axn.set_title('North (cm)')
1754+
axe.set_title('East (cm)')
1755+
axu.set_title('Up (cm)')
1756+
else:
1757+
axn.set_title('North (m)')
1758+
axe.set_title('East (m)')
1759+
axu.set_title('Up (m)')
17521760
else:
1753-
axn.set_title('North (m/s)')
1754-
axe.set_title('East (m/s)')
1755-
axu.set_title('Up (m/s)')
1761+
if units=='cm':
1762+
axn.set_title('North (cm/s)')
1763+
axe.set_title('East (cm/s)')
1764+
axu.set_title('Up (cm/s)')
1765+
else:
1766+
axn.set_title('North (m/s)')
1767+
axe.set_title('East (m/s)')
1768+
axu.set_title('Up (m/s)')
17561769
if k!=len(i)-1:
17571770
axn.xaxis.set_ticklabels([])
17581771
axe.xaxis.set_ticklabels([])
@@ -1866,7 +1879,7 @@ def insar_residual(home,project_name,run_name,run_number,gflist,zlims):
18661879

18671880
matplotlib.rcParams.update({'font.size': 14})
18681881
#Decide what to plot
1869-
sta=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=0,dtype='S')
1882+
sta=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=0,dtype='U')
18701883
lon_all=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[1],dtype='f')
18711884
lat_all=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[2],dtype='f')
18721885
gf=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[7],dtype='f')
@@ -1909,7 +1922,7 @@ def insar_results(home,project_name,run_name,run_number,gflist,zlims,cmap,figsiz
19091922

19101923
matplotlib.rcParams.update({'font.size': 14})
19111924
#Decide what to plot
1912-
sta=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=0,dtype='S')
1925+
sta=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=0,dtype='U')
19131926
lon_all=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[1],dtype='f')
19141927
lat_all=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[2],dtype='f')
19151928
gf=genfromtxt(home+project_name+'/data/station_info/'+gflist,usecols=[7],dtype='f')

0 commit comments

Comments
 (0)