In [1]:
%matplotlib inline
import os
import re
import sys
import numpy
import matplotlib.pyplot as plt
try:
import cPickle as pickle
except:
import pickle
import matplotlib.patheffects as PathEffects
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection
sys.path.append('/Users/mpagani/Projects/hmtk/')
sys.path.append('/Users/mpagani/Projects/original/oq-engine/')
from hmtk.subduction.cross_sections import CrossSection, Trench
from hmtk.seismicity.selector import CatalogueSelector
from mpl_toolkits.basemap import Basemap
In [2]:
fin = open('./../data/trench.xy', 'r')
trench = []
for line in fin:
aa = re.split('\s+', re.sub('^\s+', '', line))
trench.append((float(aa[0]), float(aa[1])))
fin.close()
trench = Trench(numpy.array(trench))
In [3]:
cat = pickle.load(open("./../data/catalogue_ext_cac.p", "rb" ))
In [4]:
minlo = -110
minla = 5
maxlo = -75
maxla = 25
midlo = -100
midla = 20
In [5]:
from hmtk.parsers.catalogue.gcmt_ndk_parser import ParseNDKtoGCMT
from obspy.imaging.beachball import beach
gcmt_filename = '/Users/mpagani/Data/catalogues/gcmt/jan76_dec13.ndk'
gcmtc = ParseNDKtoGCMT(gcmt_filename)
gcmtc.read_file()
def plot_nodal_planes(catalogue, ax, minlo, minla, maxlo, maxla):
beach1 = beach(np1, xy=(-70, 80), width=30)
beach2 = beach(mt, xy=(50, 50), width=50)
ax.add_collection(beach1)
ax.add_collection(beach2)
In [8]:
fig = plt.figure(figsize=(12,9))
#
# Plot the basemap
m = Basemap(llcrnrlon=minlo, llcrnrlat=minla,
urcrnrlon=maxlo, urcrnrlat=maxla,
resolution='i', projection='tmerc',
lon_0=midlo, lat_0=midla)
#
# Draw paralleles and meridians with labels
# labels = [left,right,top,bottom]
m.drawcoastlines()
m.drawmeridians(numpy.arange(numpy.floor(minlo/10.)*10,
numpy.ceil(maxlo/10.)*10,5.),
labels=[False, False, False, True])
m.drawparallels(numpy.arange(numpy.floor(minla/10.)*10,
numpy.ceil(maxla/10.)*10,5.),
labels=[True, False, False, False])
#
# Plot the instrumental catalogue
xa, ya = m(cat.data['longitude'], cat.data['latitude'])
szea = (cat.data['magnitude']*100)**1.5
patches = []
for x, y, sze in zip(list(xa), list(ya), szea):
circle = Circle((x, y), sze, ec='white')
patches.append(circle)
print ('depths: %f %f ' % (min(cat.data['depth']), max(cat.data['depth'])))
colors = cat.data['depth']
p = PatchCollection(patches, zorder=6, edgecolors='white')
p.set_alpha(0.5)
p.set_clim([0, 200])
p.set_array(numpy.array(colors))
plt.gca().add_collection(p)
plt.colorbar(p,fraction=0.02, pad=0.04, extend='max')
#
# GCMT
x, y = m(gcmtc.catalogue.data['longitude'],
gcmtc.catalogue.data['latitude'])
#plt.plot(x, y, 'sr', zorder=10, alpha=.5)
#
# Plot the traces of cross-sections
distance = 100
cs_len = 400
ts = trench.resample(distance)
fou = open('cs_traces.csv', 'w')
x, y = m(trench.axis[:, 0], trench.axis[:, 1])
plt.plot(x, y, '-g', linewidth=2, zorder=10)
x, y = m(ts.axis[:, 0], ts.axis[:, 1])
plt.plot(x, y, '--y', linewidth=4, zorder=20)
for idx, cs in enumerate(trench.iterate_cross_sections(distance, cs_len)):
if cs is not None:
x, y = m(cs.plo, cs.pla)
plt.plot(x, y, ':r', linewidth=2, zorder=20)
text = plt.text(x[-1], y[-1], '%d' % idx, ha='center', va='center', size=10, zorder=30)
text.set_path_effects([PathEffects.withStroke(linewidth=3, foreground="w")])
tmps = '%f %f %f %f %d\n' % (cs.plo[0], cs.pla[0], cs_len, cs.strike[0], idx)
print (tmps.rstrip())
fou.write(tmps)
fou.close()
In [ ]: