cesca
/
seiscloud
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added option for global plots, instead of regional. Default is regional. Code and manual updated.

Browse Source
master
cesca 3 years ago
parent 6c66963607
commit e3db1be89c
3 changed files with 90 additions and 154 deletions
Show Stats Download Patch File Download Diff File

BIN
gx/seiscloud_manual.odt
View File

Binary file not shown.

BIN
gx/seiscloud_manual.pdf
View File

Binary file not shown.

244
src/plot.py
Unescape Escape View File

@ -11,9 +11,11 @@ from pyrocko.automap import Map
from pyrocko import orthodrome as od
from pyrocko import moment_tensor as pmt
from pyrocko.plot import beachball
from pyrocko.plot.gmtpy import GMT, cm
from matplotlib import dates, colors
from seiscloud import cluster as scc
epsilon = 1e-6
km = 1000.
@ -214,56 +216,29 @@ def color2rgb(col):
return rgb
def global_spatial(events, eventsclusters, clusters, conf, plotdir):
def plot_global(events, eventsclusters, clusters, conf, plotdir):
'''
Plot map of seismicity clusters
'''
# lats = [ev.lat for ev in events]
# lons = [ev.lon for ev in events]
latmean = 0.5*(conf.latmin+conf.latmax)
if conf.lonmin <= conf.lonmax:
lonmean = 0.5*(conf.lonmin+conf.lonmax)
else:
lonmean = 0.5*(conf.lonmin+(360.+conf.lonmax))
if lonmean > 180.:
lonmean = lonmean - 360.
center = od.Loc(latmean, lonmean)
# Map size
if conf.map_radius is not None:
safe_radius = conf.map_radius
else:
corners = [od.Loc(conf.latmin, conf.lonmin),
od.Loc(conf.latmin, conf.lonmax)]
dist1 = od.distance_accurate50m(center, corners[0])
dist2 = od.distance_accurate50m(center, corners[1])
safe_radius = max(dist1, dist2)
# Generate the basic map
m = Map(
lat=center.lat,
lon=center.lon,
radius=safe_radius,
width=30., height=30.,
show_grid=False,
show_topo=conf.sw_topography,
color_dry=(238, 236, 230),
topo_cpt_wet='light_sea_uniform',
topo_cpt_dry='light_land_uniform',
illuminate=True,
illuminate_factor_ocean=0.15,
show_rivers=False,
show_plates=False)
# if conf.sw_filter_events:
# 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()
gmt = GMT(config={
'MAP_FRAME_TYPE': 'fancy',
'GRID_CROSS_SIZE_PRIMARY': '0',
'PS_MEDIA': 'Custom_%ix%i' % (21*cm, 13*cm)})
gmt.psbasemap(
R='g', # global
J='H0/16c',
B='0')
gmt.pscoast(
R='g', # region
J='H0/16c', # projection
B='g30', # grid
D='c', # resolution
S=(220, 220, 250), # wet fill color
G=(250, 250, 220), # dry fill color
W='thinnest') # shoreline pen
# Events in clusters
for id_cluster in clusters:
@ -273,12 +248,15 @@ def global_spatial(events, eventsclusters, clusters, conf, plotdir):
if evcl == id_cluster:
mylats.append(events[iev].lat)
mylons.append(events[iev].lon)
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',
G=color2rgb(col), *m.jxyr)
gmt.psxy(
R='g',
J='H0/6i',
in_columns=(mylons, mylats),
S='c3p',
G=color2rgb(col))
figname = os.path.join(plotdir, 'plot_map.'+conf.figure_format)
m.save(figname)
# m.show()
gmt.save(figname)
def plot_global_dbscan(events, eventsclusters, clusters, conf, plotdir):
@ -286,43 +264,24 @@ def plot_global_dbscan(events, eventsclusters, clusters, conf, plotdir):
Plot map of seismicity dbscan types
'''
latmean = 0.5*(conf.latmin+conf.latmax)
if conf.lonmin <= conf.lonmax:
lonmean = 0.5*(conf.lonmin+conf.lonmax)
else:
lonmean = 0.5*(conf.lonmin+(360.+conf.lonmax))
if lonmean > 180.:
lonmean = lonmean - 360.
center = od.Loc(latmean, lonmean)
# Map size
if conf.map_radius is not None:
safe_radius = conf.map_radius
else:
corners = [od.Loc(conf.latmin, conf.lonmin),
od.Loc(conf.latmin, conf.lonmax)]
dist1 = od.distance_accurate50m(center, corners[0])
dist2 = od.distance_accurate50m(center, corners[1])
safe_radius = max(dist1, dist2)
# Generate the basic map
m = Map(
lat=center.lat,
lon=center.lon,
radius=safe_radius,
width=30., height=30.,
show_grid=False,
show_topo=conf.sw_topography,
color_dry=(238, 236, 230),
topo_cpt_wet='light_sea_uniform',
topo_cpt_dry='light_land_uniform',
illuminate=True,
illuminate_factor_ocean=0.15,
show_rivers=False,
show_plates=False)
# Draw some larger cities covered by the map area
m.draw_cities()
gmt = GMT(config={
'MAP_FRAME_TYPE': 'fancy',
'GRID_CROSS_SIZE_PRIMARY': '0',
'PS_MEDIA': 'Custom_%ix%i' % (21*cm, 13*cm)})
gmt.psbasemap(
R='g', # global
J='H0/16c',
B='0')
gmt.pscoast(
R='g', # region
J='H0/16c', # projection
B='g30', # grid
D='c', # resolution
S=(220, 220, 250), # wet fill color
G=(250, 250, 220), # dry fill color
W='thinnest') # shoreline pen
# Events in clusters
mylats, mylons = [], []
@ -331,8 +290,12 @@ def plot_global_dbscan(events, eventsclusters, clusters, conf, plotdir):
if 'dbtype:edge' in events[iev].tags:
mylats.append(events[iev].lat)
mylons.append(events[iev].lon)
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',
G=color2rgb(col), *m.jxyr)
gmt.psxy(
R='g',
J='H0/6i',
in_columns=(mylons, mylats),
S='c3p',
G=color2rgb(col))
mylats, mylons = [], []
col = 'dodgerblue'
@ -340,8 +303,12 @@ def plot_global_dbscan(events, eventsclusters, clusters, conf, plotdir):
if 'dbtype:core' in events[iev].tags:
mylats.append(events[iev].lat)
mylons.append(events[iev].lon)
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',
G=color2rgb(col), *m.jxyr)
gmt.psxy(
R='g',
J='H0/6i',
in_columns=(mylons, mylats),
S='c3p',
G=color2rgb(col))
mylats, mylons = [], []
col = 'gray'
@ -349,67 +316,42 @@ def plot_global_dbscan(events, eventsclusters, clusters, conf, plotdir):
if 'dbtype:isle' in events[iev].tags:
mylats.append(events[iev].lat)
mylons.append(events[iev].lon)
m.gmt.psxy(in_columns=(mylons, mylats), S='c7p',
G=color2rgb(col), *m.jxyr)
gmt.psxy(
R='g',
J='H0/6i',
in_columns=(mylons, mylats),
S='c3p',
G=color2rgb(col))
figname = os.path.join(plotdir, 'plot_map_dbscan.'+conf.figure_format)
m.save(figname)
# m.show()
gmt.save(figname)
def plot_global_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]
latmean = 0.5*(conf.latmin+conf.latmax)
if conf.lonmin <= conf.lonmax:
lonmean = 0.5*(conf.lonmin+conf.lonmax)
else:
lonmean = 0.5*(conf.lonmin+(360.+conf.lonmax))
if lonmean > 180.:
lonmean = lonmean - 360.
center = od.Loc(latmean, lonmean)
# Map size
if conf.sw_manual_radius:
safe_radius = conf.map_radius
else:
corners = [od.Loc(conf.latmin, conf.lonmin),
od.Loc(conf.latmin, conf.lonmax)]
dist1 = od.distance_accurate50m(center, corners[0])
dist2 = od.distance_accurate50m(center, corners[1])
safe_radius = max(dist1, dist2)
# Generate the basic map
m = Map(
lat=center.lat,
lon=center.lon,
radius=safe_radius,
width=30., height=30.,
show_grid=False,
show_topo=conf.sw_topography,
color_dry=(238, 236, 230),
topo_cpt_wet='light_sea_uniform',
topo_cpt_dry='light_land_uniform',
illuminate=True,
illuminate_factor_ocean=0.15,
show_rivers=False,
show_plates=False)
# if conf.sw_filter_events:
# 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()
gmt = GMT(config={
'MAP_FRAME_TYPE': 'fancy',
'GRID_CROSS_SIZE_PRIMARY': '0',
'PS_MEDIA': 'Custom_%ix%i' % (21*cm, 13*cm)})
gmt.psbasemap(
R='g', # global
J='H0/16c',
B='0')
gmt.pscoast(
R='g', # region
J='H0/16c', # projection
B='g30', # grid
D='c', # resolution
S=(220, 220, 250), # wet fill color
G=(250, 250, 220), # dry fill color
W='thinnest') # shoreline pen
# Events in clusters
factor_symbl_size = 5.
factor_symbl_size = 2.5
beachball_symbol = 'd'
for id_cluster in clusters:
@ -421,34 +363,27 @@ def plot_global_with_dcs(events, eventsclusters, clusters, conf, plotdir):
if ev.moment_tensor is not None:
factor_symbl_size = ev.magnitude
devi = ev.moment_tensor.deviatoric()
beachball_size = 3.*factor_symbl_size
beachball_size = 1.*factor_symbl_size
mt = devi.m_up_south_east()
mt = mt / ev.moment_tensor.scalar_moment() \
* pmt.magnitude_to_moment(5.0)
m6 = pmt.to6(mt)
if m.gmt.is_gmt5():
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))
data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0)
m.gmt.psmeca(
gmt.psmeca(
R='g',
J='H0/6i',
in_rows=[data],
G=g_col,
E='white',
W='1p,%s' % g_col,
*m.jxyr,
**kwargs)
M=True,
S='%s%g' % (beachball_symbol[0],
(beachball_size)/gmtpy.cm))
figname = os.path.join(plotdir, 'plot_map_with_dcs.'+conf.figure_format)
m.save(figname)
# m.show()
gmt.save(figname)
def plot_spatial(events, eventsclusters, clusters, conf, plotdir):
@ -1237,6 +1172,7 @@ def plot_all(events, eventsclusters, clusters, conf, resdir, plotdir):
except:
print('skipped')
plot_global(events, eventsclusters, clusters, conf, plotdir)
if conf.sw_global_plots:
print("spatial plot...")
try:

Write Preview
Loading…
Cancel
Save