|
|
|
|
@ -1,68 +1,67 @@
|
|
|
|
|
'''This module provides basic cluster processing for seismic events.'''
|
|
|
|
|
|
|
|
|
|
import sys, os, datetime, math
|
|
|
|
|
import sys
|
|
|
|
|
import os
|
|
|
|
|
import datetime
|
|
|
|
|
import math
|
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
import numpy as num
|
|
|
|
|
from pyrocko import util, model, gmtpy
|
|
|
|
|
from pyrocko import model, gmtpy
|
|
|
|
|
from pyrocko.automap import Map
|
|
|
|
|
from pyrocko import orthodrome as od
|
|
|
|
|
from pyrocko import moment_tensor as pmt
|
|
|
|
|
from matplotlib import collections as mc
|
|
|
|
|
from matplotlib import dates, colors
|
|
|
|
|
|
|
|
|
|
epsilon = 1e-6
|
|
|
|
|
km = 1000.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def CartesianToLambert(e,n,u):
|
|
|
|
|
def CartesianToLambert(e, n, u):
|
|
|
|
|
'''
|
|
|
|
|
Convert axis orientation to coordinate for axis plot
|
|
|
|
|
(using Lambert azimuthal equal-area projection)
|
|
|
|
|
'''
|
|
|
|
|
# Using Lambert azimuthal equal-area projection
|
|
|
|
|
x=(1/math.sqrt(2.))*e*math.sqrt(2/(1-u))
|
|
|
|
|
y=(1/math.sqrt(2.))*n*math.sqrt(2/(1-u))
|
|
|
|
|
return x,y
|
|
|
|
|
|
|
|
|
|
x = (1/math.sqrt(2.))*e*math.sqrt(2/(1-u))
|
|
|
|
|
y = (1/math.sqrt(2.))*n*math.sqrt(2/(1-u))
|
|
|
|
|
return x, y
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def get_axis_coords(events):
|
|
|
|
|
'''
|
|
|
|
|
Get axis plot coordinates for all events
|
|
|
|
|
'''
|
|
|
|
|
xs,ys=[],[]
|
|
|
|
|
xs, ys = [], []
|
|
|
|
|
|
|
|
|
|
for ev in events:
|
|
|
|
|
pax=ev.moment_tensor.p_axis()
|
|
|
|
|
tax=ev.moment_tensor.t_axis()
|
|
|
|
|
bax=ev.moment_tensor.null_axis()
|
|
|
|
|
p=[pax[0,0],pax[0,1],pax[0,2]]
|
|
|
|
|
t=[tax[0,0],tax[0,1],tax[0,2]]
|
|
|
|
|
b=[bax[0,0],bax[0,1],bax[0,2]]
|
|
|
|
|
|
|
|
|
|
if p[2]<0:
|
|
|
|
|
p=[-p[0],-p[1],-p[2]]
|
|
|
|
|
if t[2]<0:
|
|
|
|
|
t=[-t[0],-t[1],-t[2]]
|
|
|
|
|
if b[2]<0:
|
|
|
|
|
b=[-b[0],-b[1],-b[2]]
|
|
|
|
|
pax = ev.moment_tensor.p_axis()
|
|
|
|
|
tax = ev.moment_tensor.t_axis()
|
|
|
|
|
bax = ev.moment_tensor.null_axis()
|
|
|
|
|
p = [pax[0, 0], pax[0, 1], pax[0, 2]]
|
|
|
|
|
t = [tax[0, 0], tax[0, 1], tax[0, 2]]
|
|
|
|
|
b = [bax[0, 0], bax[0, 1], bax[0, 2]]
|
|
|
|
|
|
|
|
|
|
px,py=CartesianToLambert(p[1],p[0],-p[2])
|
|
|
|
|
tx,ty=CartesianToLambert(t[1],t[0],-t[2])
|
|
|
|
|
bx,by=CartesianToLambert(b[1],b[0],-b[2])
|
|
|
|
|
if p[2] < 0:
|
|
|
|
|
p = [-p[0], -p[1], -p[2]]
|
|
|
|
|
if t[2] < 0:
|
|
|
|
|
t = [-t[0], -t[1], -t[2]]
|
|
|
|
|
if b[2] < 0:
|
|
|
|
|
b = [-b[0], -b[1], -b[2]]
|
|
|
|
|
|
|
|
|
|
xs.append([px,tx,bx])
|
|
|
|
|
ys.append([py,ty,by])
|
|
|
|
|
px, py = CartesianToLambert(p[1], p[0], -p[2])
|
|
|
|
|
tx, ty = CartesianToLambert(t[1], t[0], -t[2])
|
|
|
|
|
bx, by = CartesianToLambert(b[1], b[0], -b[2])
|
|
|
|
|
|
|
|
|
|
return xs,ys
|
|
|
|
|
xs.append([px, tx, bx])
|
|
|
|
|
ys.append([py, ty, by])
|
|
|
|
|
|
|
|
|
|
return xs, ys
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def getCoordinatesHudsonPlot(mts):
|
|
|
|
|
'''
|
|
|
|
|
Get Hudson plot coordinates for all events
|
|
|
|
|
'''
|
|
|
|
|
us,vs=[],[]
|
|
|
|
|
us, vs = [], []
|
|
|
|
|
for mt in mts:
|
|
|
|
|
mt = pmt.values_to_matrix(mt)
|
|
|
|
|
eig_m = pmt.eigh_check(mt)[0]
|
|
|
|
|
@ -74,61 +73,63 @@ def getCoordinatesHudsonPlot(mts):
|
|
|
|
|
return us, vs
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def get_triangle_coords(events):
|
|
|
|
|
'''
|
|
|
|
|
Get triangle plot coordinates for all events
|
|
|
|
|
'''
|
|
|
|
|
xs,ys,cs=[],[],[]
|
|
|
|
|
xs, ys, cs = [], [], []
|
|
|
|
|
|
|
|
|
|
for ev in events:
|
|
|
|
|
pax=ev.moment_tensor.p_axis()
|
|
|
|
|
tax=ev.moment_tensor.t_axis()
|
|
|
|
|
bax=ev.moment_tensor.null_axis()
|
|
|
|
|
p=[pax[0,0],pax[0,1],pax[0,2]]
|
|
|
|
|
t=[tax[0,0],tax[0,1],tax[0,2]]
|
|
|
|
|
b=[bax[0,0],bax[0,1],bax[0,2]]
|
|
|
|
|
pax = ev.moment_tensor.p_axis()
|
|
|
|
|
tax = ev.moment_tensor.t_axis()
|
|
|
|
|
bax = ev.moment_tensor.null_axis()
|
|
|
|
|
p = [pax[0, 0], pax[0, 1], pax[0, 2]]
|
|
|
|
|
t = [tax[0, 0], tax[0, 1], tax[0, 2]]
|
|
|
|
|
b = [bax[0, 0], bax[0, 1], bax[0, 2]]
|
|
|
|
|
|
|
|
|
|
# t,p,b=sds2tpb(ev.strike,ev.dip,ev.rake)
|
|
|
|
|
a3526=35.26*math.pi/180.
|
|
|
|
|
a3526 = 35.26*math.pi/180.
|
|
|
|
|
|
|
|
|
|
deltab=math.asin(abs(b[2]))
|
|
|
|
|
deltat=math.asin(abs(t[2]))
|
|
|
|
|
deltap=math.asin(abs(p[2]))
|
|
|
|
|
deltab = math.asin(abs(b[2]))
|
|
|
|
|
deltat = math.asin(abs(t[2]))
|
|
|
|
|
deltap = math.asin(abs(p[2]))
|
|
|
|
|
|
|
|
|
|
psi=math.atan(math.sin(deltat)/math.sin(deltap))-0.25*math.pi
|
|
|
|
|
den=math.sin(a3526)*math.sin(deltab)+math.cos(a3526)*math.cos(deltab)*math.cos(psi)
|
|
|
|
|
psi = math.atan(math.sin(deltat)/math.sin(deltap))-0.25*math.pi
|
|
|
|
|
den = math.sin(a3526)*math.sin(deltab) +\
|
|
|
|
|
math.cos(a3526)*math.cos(deltab)*math.cos(psi)
|
|
|
|
|
|
|
|
|
|
if abs(den)<=epsilon:
|
|
|
|
|
den=epsilon
|
|
|
|
|
if abs(den) <= epsilon:
|
|
|
|
|
den = epsilon
|
|
|
|
|
|
|
|
|
|
x=math.cos(deltab)*math.sin(psi)/den
|
|
|
|
|
y=(math.cos(a3526)*math.sin(deltab)-math.sin(a3526)*math.cos(deltab)*math.cos(psi))/den
|
|
|
|
|
s2p,s2t,s2b=(math.sin(deltap))**2,(math.sin(deltat))**2,(math.sin(deltab))**2
|
|
|
|
|
c=(s2t,s2b,s2p)
|
|
|
|
|
x = math.cos(deltab)*math.sin(psi)/den
|
|
|
|
|
y = (math.cos(a3526)*math.sin(deltab) -
|
|
|
|
|
math.sin(a3526)*math.cos(deltab)*math.cos(psi))/den
|
|
|
|
|
s2p = (math.sin(deltap))**2
|
|
|
|
|
s2t = (math.sin(deltat))**2
|
|
|
|
|
s2b = (math.sin(deltab))**2
|
|
|
|
|
c = (s2t, s2b, s2p)
|
|
|
|
|
|
|
|
|
|
xs.append(x)
|
|
|
|
|
ys.append(y)
|
|
|
|
|
cs.append(c)
|
|
|
|
|
|
|
|
|
|
return xs,ys,cs
|
|
|
|
|
|
|
|
|
|
return xs, ys, cs
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def view_and_savefig_similarity_matrix(simmat,figname,title):
|
|
|
|
|
def view_and_savefig_similarity_matrix(simmat, figname, title):
|
|
|
|
|
'''
|
|
|
|
|
View and save a similarity matrix
|
|
|
|
|
'''
|
|
|
|
|
nev=len(simmat)
|
|
|
|
|
nev = len(simmat)
|
|
|
|
|
|
|
|
|
|
f = plt.figure(1,figsize=(12,10), facecolor='w', edgecolor='k')
|
|
|
|
|
f = plt.figure(1, figsize=(12, 10), facecolor='w', edgecolor='k')
|
|
|
|
|
grid = plt.GridSpec(10, 10, wspace=0.2, hspace=0.2)
|
|
|
|
|
|
|
|
|
|
plt.subplot(grid[0:9, 0:9])
|
|
|
|
|
|
|
|
|
|
plt.imshow(simmat,interpolation='none',cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax = nev, xmin = 0)
|
|
|
|
|
plt.ylim(ymax = nev, ymin = 0)
|
|
|
|
|
plt.imshow(simmat, interpolation='none', cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax=nev, xmin=0)
|
|
|
|
|
plt.ylim(ymax=nev, ymin=0)
|
|
|
|
|
|
|
|
|
|
plt.xlabel("Event number")
|
|
|
|
|
plt.ylabel("Event number")
|
|
|
|
|
@ -137,21 +138,20 @@ def view_and_savefig_similarity_matrix(simmat,figname,title):
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def view_similarity_matrix(simmat,title):
|
|
|
|
|
def view_similarity_matrix(simmat, title):
|
|
|
|
|
'''
|
|
|
|
|
View (only) a similarity matrix
|
|
|
|
|
'''
|
|
|
|
|
nev=len(simmat)
|
|
|
|
|
nev = len(simmat)
|
|
|
|
|
|
|
|
|
|
f = plt.figure(1,figsize=(12,10), facecolor='w', edgecolor='k')
|
|
|
|
|
plt.figure(1, figsize=(12, 10), facecolor='w', edgecolor='k')
|
|
|
|
|
grid = plt.GridSpec(10, 10, wspace=0.2, hspace=0.2)
|
|
|
|
|
|
|
|
|
|
plt.subplot(grid[0:9, 0:9])
|
|
|
|
|
|
|
|
|
|
plt.imshow(simmat,interpolation='none',cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax = nev, xmin = 0)
|
|
|
|
|
plt.ylim(ymax = nev, ymin = 0)
|
|
|
|
|
plt.imshow(simmat, interpolation='none', cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax=nev, xmin=0)
|
|
|
|
|
plt.ylim(ymax=nev, ymin=0)
|
|
|
|
|
|
|
|
|
|
plt.xlabel("Event number")
|
|
|
|
|
plt.ylabel("Event number")
|
|
|
|
|
@ -159,21 +159,20 @@ def view_similarity_matrix(simmat,title):
|
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def savefig_similarity_matrix(simmat,figname,title):
|
|
|
|
|
def savefig_similarity_matrix(simmat, figname, title):
|
|
|
|
|
'''
|
|
|
|
|
Save (only) a similarity matrix
|
|
|
|
|
'''
|
|
|
|
|
nev=len(simmat)
|
|
|
|
|
nev = len(simmat)
|
|
|
|
|
|
|
|
|
|
f = plt.figure(1,figsize=(12,10), facecolor='w', edgecolor='k')
|
|
|
|
|
f = plt.figure(1, figsize=(12, 10), facecolor='w', edgecolor='k')
|
|
|
|
|
grid = plt.GridSpec(10, 10, wspace=0.2, hspace=0.2)
|
|
|
|
|
|
|
|
|
|
plt.subplot(grid[0:9, 0:9])
|
|
|
|
|
|
|
|
|
|
plt.imshow(simmat,interpolation='none',cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax = nev, xmin = 0)
|
|
|
|
|
plt.ylim(ymax = nev, ymin = 0)
|
|
|
|
|
plt.imshow(simmat, interpolation='none', cmap='GnBu_r')
|
|
|
|
|
plt.xlim(xmax=nev, xmin=0)
|
|
|
|
|
plt.ylim(ymax=nev, ymin=0)
|
|
|
|
|
|
|
|
|
|
plt.xlabel("Event number")
|
|
|
|
|
plt.ylabel("Event number")
|
|
|
|
|
@ -181,68 +180,64 @@ def savefig_similarity_matrix(simmat,figname,title):
|
|
|
|
|
f.savefig(figname)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def cluster_to_color(cluster_id):
|
|
|
|
|
'''
|
|
|
|
|
Given a cluster id provide the corresponding color
|
|
|
|
|
'''
|
|
|
|
|
mycolors=['black','red','blue','green','darkviolet','gold','darkorange',\
|
|
|
|
|
'dodgerblue','brown','lightseagreen','plum','lawngreen',\
|
|
|
|
|
'palevioletred','royalblue','limegreen','indigo']
|
|
|
|
|
my_def_color='gainsboro'
|
|
|
|
|
if cluster_id>(len(mycolors)-2):
|
|
|
|
|
color=def_color
|
|
|
|
|
mycolors = ['black', 'red', 'blue', 'green', 'darkviolet', 'gold',
|
|
|
|
|
'darkorange', 'dodgerblue', 'brown', 'lightseagreen', 'plum',
|
|
|
|
|
'lawngreen', 'palevioletred', 'royalblue', 'limegreen',
|
|
|
|
|
'indigo']
|
|
|
|
|
my_def_color = 'gainsboro'
|
|
|
|
|
if cluster_id > (len(mycolors)-2):
|
|
|
|
|
color = my_def_color
|
|
|
|
|
else:
|
|
|
|
|
color=mycolors[cluster_id+1]
|
|
|
|
|
color = mycolors[cluster_id+1]
|
|
|
|
|
return color
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def color2rgb(col):
|
|
|
|
|
'''
|
|
|
|
|
Return a red/green/blue string for GMT for a given color name
|
|
|
|
|
'''
|
|
|
|
|
colarray = colors.to_rgb(col)
|
|
|
|
|
r,g,b=int(255*colarray[0]),int(255*colarray[1]),int(255*colarray[2])
|
|
|
|
|
rgb=str(r)+'/'+str(g)+'/'+str(b)
|
|
|
|
|
r, g, b = int(255*colarray[0]), int(255*colarray[1]), int(255*colarray[2])
|
|
|
|
|
rgb = str(r)+'/'+str(g)+'/'+str(b)
|
|
|
|
|
return rgb
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def plot_spatial(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
def plot_spatial(events, eventsclusters, clusters, conf, plotdir):
|
|
|
|
|
'''
|
|
|
|
|
Plot map of seismicity clusters
|
|
|
|
|
'''
|
|
|
|
|
lats=[ev.lat for ev in events]
|
|
|
|
|
lons=[ev.lon for ev in events]
|
|
|
|
|
deps=[ev.depth for ev in events]
|
|
|
|
|
colors = [cluster_to_color(clid) for clid in eventsclusters]
|
|
|
|
|
lats = [ev.lat for ev in events]
|
|
|
|
|
lons = [ev.lon for ev in events]
|
|
|
|
|
# deps = [ev.depth for ev in events]
|
|
|
|
|
# colors = [cluster_to_color(clid) for clid in eventsclusters]
|
|
|
|
|
|
|
|
|
|
if conf.sw_filterevent:
|
|
|
|
|
latmin,latmax=cong.latmin,conf.latmax
|
|
|
|
|
lonmin,lonmax=cong.lonmin,conf.lonmax
|
|
|
|
|
depmin,depmax=cong.depthmin,conf.depthmax
|
|
|
|
|
if latmin>latmax:
|
|
|
|
|
print ('cases over lon +-180 still to be implemented')
|
|
|
|
|
latmin, latmax = conf.latmin, conf.latmax
|
|
|
|
|
lonmin, lonmax = conf.lonmin, conf.lonmax
|
|
|
|
|
# depmin, depmax = conf.depthmin, conf.depthmax
|
|
|
|
|
if latmin > latmax:
|
|
|
|
|
print('cases over lon +-180 still to be implemented')
|
|
|
|
|
sys.exit()
|
|
|
|
|
center=model.event(0.5*(latmin+latmax),0.5*(lonmin+lonmax))
|
|
|
|
|
center = model.event(0.5*(latmin+latmax), 0.5*(lonmin+lonmax))
|
|
|
|
|
else:
|
|
|
|
|
latmin=min([ev.lat for ev in events])-0.1
|
|
|
|
|
latmax=max([ev.lat for ev in events])+0.1
|
|
|
|
|
lonmin=min([ev.lon for ev in events])-0.1
|
|
|
|
|
lonmax=max([ev.lon for ev in events])+0.1
|
|
|
|
|
depmin=0.*km
|
|
|
|
|
depmax=1.05*max([ev.depth for ev in events])
|
|
|
|
|
center=od.Loc(0.5*(latmin+latmax),0.5*(lonmin+lonmax))
|
|
|
|
|
latmin, latmax = min(lats)-0.1, max(lats)+0.1
|
|
|
|
|
lonmin, lonmax = min(lons)-0.1, max(lons)+0.1
|
|
|
|
|
# depmin = 0.*km
|
|
|
|
|
# depmax = 1.05*max(deps)
|
|
|
|
|
center = od.Loc(0.5*(latmin+latmax), 0.5*(lonmin+lonmax))
|
|
|
|
|
|
|
|
|
|
# Map size
|
|
|
|
|
if conf.map_radius is not None:
|
|
|
|
|
safe_radius = conf.map_radius
|
|
|
|
|
else:
|
|
|
|
|
corners=[od.Loc(latmin,lonmin),od.Loc(latmin,lonmax)]
|
|
|
|
|
corners = [od.Loc(latmin, lonmin), od.Loc(latmin, lonmax)]
|
|
|
|
|
dist1 = od.distance_accurate50m(center, corners[0])
|
|
|
|
|
dist2 = od.distance_accurate50m(center, corners[1])
|
|
|
|
|
safe_radius = max(dist1,dist2)
|
|
|
|
|
safe_radius = max(dist1, dist2)
|
|
|
|
|
|
|
|
|
|
# Generate the basic map
|
|
|
|
|
m = Map(
|
|
|
|
|
@ -261,64 +256,62 @@ def plot_spatial(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
show_plates=False)
|
|
|
|
|
|
|
|
|
|
if conf.sw_filterevent:
|
|
|
|
|
rectlons=[lonmin,lonmin,lonmax,lonmax,lonmin]
|
|
|
|
|
rectlats=[latmin,latmax,latmax,latmin,latmin]
|
|
|
|
|
m.gmt.psxy(in_columns=(rectlons, rectlats), W='thin,0/0/0,dashed', *m.jxyr)
|
|
|
|
|
rectlons = [lonmin, lonmin, lonmax, lonmax, lonmin]
|
|
|
|
|
rectlats = [latmin, latmax, latmax, latmin, latmin]
|
|
|
|
|
m.gmt.psxy(in_columns=(rectlons, rectlats),
|
|
|
|
|
W='thin,0/0/0,dashed', *m.jxyr)
|
|
|
|
|
|
|
|
|
|
# Draw some larger cities covered by the map area
|
|
|
|
|
m.draw_cities()
|
|
|
|
|
|
|
|
|
|
# Events in clusters
|
|
|
|
|
for id_cluster in clusters:
|
|
|
|
|
col=cluster_to_color(id_cluster)
|
|
|
|
|
mylats,mylons=[],[]
|
|
|
|
|
for iev,evcl in enumerate(eventsclusters):
|
|
|
|
|
if evcl==id_cluster:
|
|
|
|
|
col = cluster_to_color(id_cluster)
|
|
|
|
|
mylats, mylons = [], []
|
|
|
|
|
for iev, evcl in enumerate(eventsclusters):
|
|
|
|
|
if evcl == id_cluster:
|
|
|
|
|
mylats.append(events[iev].lat)
|
|
|
|
|
mylons.append(events[iev].lon)
|
|
|
|
|
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',\
|
|
|
|
|
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',
|
|
|
|
|
G=color2rgb(col), *m.jxyr)
|
|
|
|
|
|
|
|
|
|
figname=os.path.join(plotdir,'plot_map.'+conf.figure_format)
|
|
|
|
|
figname = os.path.join(plotdir, 'plot_map.'+conf.figure_format)
|
|
|
|
|
m.save(figname)
|
|
|
|
|
# m.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def plot_spatial_with_dcs(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
def plot_spatial_with_dcs(events, eventsclusters, clusters, conf, plotdir):
|
|
|
|
|
'''
|
|
|
|
|
Plot map of seismicity clusters, with focal mechanisms
|
|
|
|
|
'''
|
|
|
|
|
lats=[ev.lat for ev in events]
|
|
|
|
|
lons=[ev.lon for ev in events]
|
|
|
|
|
deps=[ev.depth for ev in events]
|
|
|
|
|
colors = [cluster_to_color(clid) for clid in eventsclusters]
|
|
|
|
|
lats = [ev.lat for ev in events]
|
|
|
|
|
lons = [ev.lon for ev in events]
|
|
|
|
|
# deps = [ev.depth for ev in events]
|
|
|
|
|
# colors = [cluster_to_color(clid) for clid in eventsclusters]
|
|
|
|
|
|
|
|
|
|
if conf.sw_filterevent:
|
|
|
|
|
latmin,latmax=cong.latmin,conf.latmax
|
|
|
|
|
lonmin,lonmax=cong.lonmin,conf.lonmax
|
|
|
|
|
depmin,depmax=cong.depthmin,conf.depthmax
|
|
|
|
|
if latmin>latmax:
|
|
|
|
|
print ('cases over lon +-180 still to be implemented')
|
|
|
|
|
latmin, latmax = conf.latmin, conf.latmax
|
|
|
|
|
lonmin, lonmax = conf.lonmin, conf.lonmax
|
|
|
|
|
# depmin, depmax = conf.depthmin, conf.depthmax
|
|
|
|
|
if latmin > latmax:
|
|
|
|
|
print('cases over lon +-180 still to be implemented')
|
|
|
|
|
sys.exit()
|
|
|
|
|
center=model.event(0.5*(latmin+latmax),0.5*(lonmin+lonmax))
|
|
|
|
|
center = model.event(0.5*(latmin+latmax), 0.5*(lonmin+lonmax))
|
|
|
|
|
else:
|
|
|
|
|
latmin=min([ev.lat for ev in events])-0.1
|
|
|
|
|
latmax=max([ev.lat for ev in events])+0.1
|
|
|
|
|
lonmin=min([ev.lon for ev in events])-0.1
|
|
|
|
|
lonmax=max([ev.lon for ev in events])+0.1
|
|
|
|
|
depmin=0.*km
|
|
|
|
|
depmax=1.05*max([ev.depth for ev in events])
|
|
|
|
|
center=od.Loc(0.5*(latmin+latmax),0.5*(lonmin+lonmax))
|
|
|
|
|
latmin, latmax = min(lats)-0.1, max(lats)+0.1
|
|
|
|
|
lonmin, lonmax = min(lons)-0.1, max(lons)+0.1
|
|
|
|
|
# depmin = 0.*km
|
|
|
|
|
# depmax = 1.05*max(deps)
|
|
|
|
|
center = od.Loc(0.5*(latmin+latmax), 0.5*(lonmin+lonmax))
|
|
|
|
|
|
|
|
|
|
# Map size
|
|
|
|
|
if conf.map_radius is not None:
|
|
|
|
|
safe_radius = conf.map_radius
|
|
|
|
|
else:
|
|
|
|
|
corners=[od.Loc(latmin,lonmin),od.Loc(latmin,lonmax)]
|
|
|
|
|
corners = [od.Loc(latmin, lonmin), od.Loc(latmin, lonmax)]
|
|
|
|
|
dist1 = od.distance_accurate50m(center, corners[0])
|
|
|
|
|
dist2 = od.distance_accurate50m(center, corners[1])
|
|
|
|
|
safe_radius = max(dist1,dist2)
|
|
|
|
|
safe_radius = max(dist1, dist2)
|
|
|
|
|
|
|
|
|
|
# Generate the basic map
|
|
|
|
|
m = Map(
|
|
|
|
|
@ -337,9 +330,10 @@ def plot_spatial_with_dcs(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
show_plates=False)
|
|
|
|
|
|
|
|
|
|
if conf.sw_filterevent:
|
|
|
|
|
rectlons=[lonmin,lonmin,lonmax,lonmax,lonmin]
|
|
|
|
|
rectlats=[latmin,latmax,latmax,latmin,latmin]
|
|
|
|
|
m.gmt.psxy(in_columns=(rectlons, rectlats), W='thin,0/0/0,dashed', *m.jxyr)
|
|
|
|
|
rectlons = [lonmin, lonmin, lonmax, lonmax, lonmin]
|
|
|
|
|
rectlats = [latmin, latmax, latmax, latmin, latmin]
|
|
|
|
|
m.gmt.psxy(in_columns=(rectlons, rectlats),
|
|
|
|
|
W='thin,0/0/0,dashed', *m.jxyr)
|
|
|
|
|
|
|
|
|
|
# Draw some larger cities covered by the map area
|
|
|
|
|
m.draw_cities()
|
|
|
|
|
@ -349,11 +343,11 @@ def plot_spatial_with_dcs(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
beachball_symbol = 'd'
|
|
|
|
|
|
|
|
|
|
for id_cluster in clusters:
|
|
|
|
|
col=cluster_to_color(id_cluster)
|
|
|
|
|
g_col=color2rgb(col)
|
|
|
|
|
for iev,evcl in enumerate(eventsclusters):
|
|
|
|
|
if evcl==id_cluster:
|
|
|
|
|
ev=events[iev]
|
|
|
|
|
col = cluster_to_color(id_cluster)
|
|
|
|
|
g_col = color2rgb(col)
|
|
|
|
|
for iev, evcl in enumerate(eventsclusters):
|
|
|
|
|
if evcl == id_cluster:
|
|
|
|
|
ev = events[iev]
|
|
|
|
|
devi = ev.moment_tensor.deviatoric()
|
|
|
|
|
beachball_size = 3.*factor_symbl_size
|
|
|
|
|
mt = devi.m_up_south_east()
|
|
|
|
|
@ -362,16 +356,15 @@ def plot_spatial_with_dcs(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
m6 = pmt.to6(mt)
|
|
|
|
|
|
|
|
|
|
if m.gmt.is_gmt5():
|
|
|
|
|
kwargs = dict(
|
|
|
|
|
M=True,
|
|
|
|
|
S='%s%g' % (beachball_symbol[0],
|
|
|
|
|
(beachball_size) / gmtpy.cm))
|
|
|
|
|
kwargs = dict(M=True, S='%s%g' %
|
|
|
|
|
(beachball_symbol[0],
|
|
|
|
|
(beachball_size)/gmtpy.cm))
|
|
|
|
|
else:
|
|
|
|
|
kwargs = dict(
|
|
|
|
|
S='%s%g' % (beachball_symbol[0],
|
|
|
|
|
(beachball_size)*2 / gmtpy.cm))
|
|
|
|
|
kwargs = dict(S='%s%g' %
|
|
|
|
|
(beachball_symbol[0],
|
|
|
|
|
(beachball_size)*2 / gmtpy.cm))
|
|
|
|
|
|
|
|
|
|
data=(ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0)
|
|
|
|
|
data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0)
|
|
|
|
|
|
|
|
|
|
m.gmt.psmeca(
|
|
|
|
|
in_rows=[data],
|
|
|
|
|
@ -381,222 +374,218 @@ def plot_spatial_with_dcs(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
*m.jxyr,
|
|
|
|
|
**kwargs)
|
|
|
|
|
|
|
|
|
|
figname=os.path.join(plotdir,'plot_map_with_dcs.'+conf.figure_format)
|
|
|
|
|
figname = os.path.join(plotdir, 'plot_map_with_dcs.'+conf.figure_format)
|
|
|
|
|
m.save(figname)
|
|
|
|
|
# m.show()
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
def plot_tm(events,eventsclusters,clusters,conf,plotdir):
|
|
|
|
|
def plot_tm(events, eventsclusters, clusters, conf, plotdir):
|
|
|
|
|
'''
|
|
|
|
|
Plot magnitude vs time for the eismicity clusters
|
|
|
|
|
'''
|
|
|
|
|
times=[ev.time for ev in events]
|
|
|
|
|
times = [ev.time for ev in events]
|
|
|
|
|
orig_dates = [datetime.datetime.fromtimestamp(ev.time) for ev in events]
|
|
|
|
|
mpl_dates=dates.date2num(orig_dates)
|
|
|
|
|
mags=[ev.magnitude for ev in events]
|
|
|
|
|
mpl_dates = dates.date2num(orig_dates)
|
|
|
|
|
mags = [ev.magnitude for ev in events]
|
|
|
|
|
colors = [cluster_to_color(clid) for clid in eventsclusters]
|
|
|
|
|
|
|
|
|
|
dates_format = dates.DateFormatter('%Y-%m-%d')
|
|
|
|
|
|
|
|
|
|
if conf.sw_filterevent:
|
|
|
|
|
tmin,tmax=cong.tmin,conf.tmax
|
|
|
|
|
magmin,magmax=cong.magmin,conf.magmax
|
|
|
|
|
tmin, tmax = conf.tmin, conf.tmax
|
|
|
|
|
magmin, magmax = conf.magmin, conf.magmax
|
|
|
|
|
else:
|
|
|
|
|
if (max(times)-min(times))/86400. > 720.:
|
|
|
|
|
dt=86400.
|
|
|
|
|
dates_loc=dates.YearLocator()
|
|
|
|
|
dt = 86400.
|
|
|
|
|
dates_loc = dates.YearLocator()
|
|
|
|
|
elif (max(times)-min(times))/86400. > 10.:
|
|
|
|
|
dt=86400.
|
|
|
|
|
dates_loc=dates.MonthLocator()
|
|
|
|
|
dt = 86400.
|
|
|
|
|
dates_loc = dates.MonthLocator()
|
|
|
|
|
elif (max(times)-min(times))/3600. > 10.:
|
|
|
|
|
dt=3600.
|
|
|
|
|
dates_loc=dates.DayLocator()
|
|
|
|
|
dt = 3600.
|
|
|
|
|
dates_loc = dates.DayLocator()
|
|
|
|
|
else:
|
|
|
|
|
dt=1.
|
|
|
|
|
dates_loc=dates.HourLocator()
|
|
|
|
|
dt = 1.
|
|
|
|
|
dates_loc = dates.HourLocator()
|
|
|
|
|
dates_format = dates.DateFormatter('%Y-%m-%d %h:%m:%s')
|
|
|
|
|
tmin,tmax=min(times)-dt,max(times)+dt
|
|
|
|
|
magmin,magmax=min(mags)-0.1,max(mags)+0.1
|
|
|
|
|
dmin=dates.date2num(datetime.datetime.fromtimestamp(tmin))
|
|
|
|
|
dmax=dates.date2num(datetime.datetime.fromtimestamp(tmax))
|
|
|
|
|
tmin, tmax = min(times)-dt, max(times)+dt |
Simone Cesca
5 years ago