pylot/map_projection.py

177 lines
4.9 KiB
Python
Raw Normal View History

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import obspy
2017-04-11 14:32:29 +02:00
from matplotlib import cm
2017-04-11 15:16:53 +02:00
from PySide import QtCore, QtGui
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT as NavigationToolbar
#import QtPyLoT
from pylot.core.util.dataprocessing import read_metadata
2017-04-11 14:32:29 +02:00
from scipy.interpolate import griddata
#pf=QtPyLoT.main()
def onpick(event):
ind = event.ind
2017-04-11 15:16:53 +02:00
for index in ind:
print(station_names[index])
2017-04-11 14:32:29 +02:00
def get_metadata(path):
metadata=read_metadata(path)
parser=metadata[1]
return parser
2017-04-11 15:16:53 +02:00
def get_station_names_lat_lon(parser):
2017-04-11 14:32:29 +02:00
station_names=[]
lat=[]
lon=[]
2017-04-11 15:16:53 +02:00
for station in parser.stations:
station_name=station[0].station_call_letters
if not station_name in station_names:
station_names.append(station_name)
lat.append(station[0].latitude)
lon.append(station[0].longitude)
return station_names, lat, lon
2017-04-11 14:32:29 +02:00
def get_picks(pf, station_names):
picks=[]
for station in station_names:
try:
picks.append(pf.autopicks[station]['P']['mpp'])
except:
picks.append(np.nan)
return picks
def get_picks_rel(picks):
picks_rel=[]
minp = min(picks)
for pick in picks:
if type(pick) is obspy.core.utcdatetime.UTCDateTime:
pick -= minp
picks_rel.append(pick)
return picks_rel
def remove_nan_picks(picks):
picks_no_nan=[]
for pick in picks:
if not np.isnan(pick):
picks_no_nan.append(pick)
return picks_no_nan
def remove_nan_lat_lon(picks, lat, lon):
lat_no_nan=[]
lon_no_nan=[]
for index, pick in enumerate(picks):
if not np.isnan(pick):
lat_no_nan.append(lat[index])
lon_no_nan.append(lon[index])
return lat_no_nan, lon_no_nan
def get_lon_lat_dim(lon, lat):
londim = max(lon) - min(lon)
latdim = max(lat) - min(lat)
return londim, latdim
def get_x_y_dim(x, y):
xdim = max(x) - min(x)
ydim = max(y) - min(y)
return xdim, ydim
2017-04-11 15:16:53 +02:00
def init_map(projection, ax, resolution='l'):
m = Basemap(projection=projection, resolution = resolution, ax=ax)
2017-04-11 14:32:29 +02:00
m.drawmapboundary(fill_color='darkblue')
m.drawcountries()
m.drawstates()
m.fillcontinents(color='grey', lake_color='aqua')
m.drawcoastlines()
return m
def get_lat_lon_axis(lat, lon):
steplat = (max(lat)-min(lat))/250
steplon = (max(lon)-min(lon))/250
lataxis = np.arange(min(lat), max(lat), steplat)
lonaxis = np.arange(min(lon), max(lon), steplon)
return lataxis, lonaxis
def get_lat_lon_grid(lataxis, lonaxis):
longrid, latgrid = np.meshgrid(lonaxis, lataxis)
return latgrid, longrid
def draw_contour_filled(picks, longrid, latgrid, picksgrid, levels='50'):
levels = np.linspace(min(picks), max(picks), 50)
contourf = m.contourf(longrid, latgrid, picksgrid, levels, latlon=True, zorder=9)
return contourf
def annotate_ax(ax, x, y, station_names):
for index, name in enumerate(station_names):
ax.annotate(' %s' % name, xy=(x[index], y[index]), fontsize='x-small', zorder=12)
2017-04-11 14:32:29 +02:00
def connect_pick(ax, onpick):
ax.figure.canvas.mpl_connect('pick_event', onpick)
2017-04-11 14:32:29 +02:00
def add_cbar(ax, scatter, label):
cbar = ax.figure.colorbar(scatter)
cbar.set_label(label)
return cbar
2017-04-11 15:16:53 +02:00
qwidget = QtGui.QWidget()
QtGui.QVBoxLayout()
main_box = QtGui.QVBoxLayout()
qwidget.setLayout(main_box)
fig = plt.figure()
main_ax = fig.add_subplot(111)
canvas = main_ax.figure.canvas
main_box.addWidget(canvas)
toolbar = NavigationToolbar(canvas, qwidget)
main_box.addWidget(toolbar)
2017-04-11 14:32:29 +02:00
parser = get_metadata('/data/Geothermie/Insheim/STAT_INFO/MAGS2_net.dless')
2017-04-11 15:16:53 +02:00
station_names, lat, lon = get_station_names_lat_lon(parser)
2017-04-11 14:32:29 +02:00
picks = get_picks(pf, station_names)
picks_rel = get_picks_rel(picks)
2017-04-11 14:32:29 +02:00
picks_no_nan = remove_nan_picks(picks_rel)
lat_no_nan, lon_no_nan = remove_nan_lat_lon(picks_rel, lat, lon)
2017-04-11 14:32:29 +02:00
londim, latdim = get_lon_lat_dim(lon, lat)
x, y = m(lon, lat)
xdim, ydim = get_x_y_dim(x, y)
2017-04-11 15:16:53 +02:00
m = init_map(projection='mill', ax=main_ax, resolution='l')
2017-04-11 14:32:29 +02:00
lataxis, lonaxis = get_lat_lon_axis(lat, lon)
latgrid, longrid = get_lat_lon_grid(lataxis, lonaxis)
picksgrid_no_nan = griddata((lat_no_nan, lon_no_nan), picks_no_nan, (latgrid, longrid), method='linear')
contourf = draw_contour_filled(picks_no_nan, longrid, latgrid, picksgrid_no_nan)
sc = m.scatter(lon, lat, s=50, facecolor='none', latlon=True, zorder=10, picker=True, edgecolor='m', label='Not Picked')
sc_picked = m.scatter(lon_no_nan, lat_no_nan, s=50, c=picks_no_nan, latlon=True, zorder=11, label='Picked')
2017-04-11 15:16:53 +02:00
annotate_ax(main_ax, x, y, station_names)
main_ax.legend()
2017-04-11 15:16:53 +02:00
connect_pick(main_ax, onpick)
2017-04-11 14:32:29 +02:00
2017-04-11 15:16:53 +02:00
cbar = add_cbar(main_ax, sc_picked, label='Time relative to first onset [s]')
2017-04-11 14:32:29 +02:00
2017-04-11 15:16:53 +02:00
canvas.draw()
2017-04-11 14:32:29 +02:00
2017-04-11 15:16:53 +02:00
qwidget.show()
2017-04-11 14:32:29 +02:00
# ax.set_xlim(min(x)-0.5*xdim, max(x)+0.5*xdim)
# ax.set_ylim(min(y)-0.5*ydim, max(y)+0.5*ydim)
2017-04-11 15:16:53 +02:00
#plt.show()