In [8]:
"""
Description:
Draft code for demonstrating data-extraction and visualization of an
ASDF file contianing MT data
References:
CreationDate: 31/07/18
Developer: rakib.hassan@ga.gov.au
Revision History:
LastUpdate: 31/07/18 RH
LastUpdate: dd/mm/yyyy Who Optional description
"""
import sys, os, math
import numpy as np
import scipy
import matplotlib.pyplot as plt
import glob
from mpl_toolkits.basemap import Basemap
from descartes import PolygonPatch
from shapely.geometry import Polygon
from collections import defaultdict
from obspy.core import inventory, Stream, UTCDateTime, read
import pyasdf
import configparser, fnmatch
import pandas as pd
import mpld3
In [9]:
# Define file-path
#ASDF_PATH = 'U:/Software/mtpy/development/asdf/au.vic2.h5'
ASDF_PATH = '/g/data/ha3/rakib/ausLAMP/Data/Output/au.vic.h5'
In [10]:
ds = pyasdf.ASDFDataSet(ASDF_PATH, mode='r')
In [11]:
coordsDict = ds.get_all_coordinates()
In [12]:
def drawBBox( minLon, minLat, maxLon, maxLat, bm, **kwargs):
bblons = np.array([minLon, maxLon, maxLon, minLon, minLon])
bblats = np.array([minLat, minLat, maxLat, maxLat, minLat])
x, y = bm( bblons, bblats )
xy = zip(x,y)
poly = Polygon(xy)
bm.ax.add_patch(PolygonPatch(poly, **kwargs))
#end func
fig = plt.figure(figsize=(10,10))
minLon = 1e32
maxLon = -1e32
minLat = 1e32
maxLat = -1e32
for k in coordsDict.keys():
coords = coordsDict[k]
lon,lat = coords['longitude'], coords['latitude'],
if(minLon > lon): minLon = lon
if(maxLon < lon): maxLon = lon
if(minLat > lat): minLat = lat
if(maxLat < lat): maxLat = lat
# end for
m = Basemap(width=800000,height=800000,projection='lcc',
resolution='l',lat_1=minLat,lat_2=maxLat,
lat_0=(minLat+maxLat)/2., lon_0=(minLon + maxLon)/2.)
# draw coastlines.
m.drawcoastlines()
#draw grid
parallels = np.linspace(np.floor(minLat)-5, np.ceil(maxLat)+5, \
int((np.ceil(maxLat)+5) - (np.floor(minLat)-5))+1)
m.drawparallels(parallels,labels=[True,True,False,False])
meridians = np.linspace(np.floor(minLon)-5, np.ceil(maxLon)+5, \
int((np.ceil(maxLon)+5) - (np.floor(minLon)-5))+1)
m.drawmeridians(meridians,labels=[False,False,True,True])
# plot stations
for k in coordsDict.keys():
coords = coordsDict[k]
lon,lat = coords['longitude'], coords['latitude'],
px,py = m(lon, lat)
pxl,pyl = m(lon, lat-0.1)
m.scatter(px, py, 50, marker='v', c='g', edgecolor='none')
plt.annotate(k, xy=(px, py), fontsize=5)
# end for
insetAx = fig.add_axes([0.75,0.75,0.125,0.125])
mInset = Basemap(resolution='c', # c, l, i, h, f or None
ax=insetAx,
projection='merc',
lat_0=-20, lon_0=132,
llcrnrlon=110, llcrnrlat=-40, urcrnrlon=155, urcrnrlat=-10)
#mInset.drawcoastlines()
mInset.fillcontinents(color='lightgray')
mInset.drawstates(color="grey")
drawBBox(minLon, minLat, maxLon, maxLat, mInset, fill='True', facecolor='k')
plt.savefig('stations.pdf')