In [6]:
from __future__ import division
import netCDF4 as NC
import numpy as np
import os
import scipy.io as sio
%matplotlib inline
import matplotlib.pyplot as plt
from salishsea_tools import viz_tools
from matplotlib.path import Path
import matplotlib.patches as patches
import pandas as pd
Coastline information:
/ocean/rich/home/matlab/m_mapWK/pb/britishcolumbia2pts.txt
Washington coastline
/ocean/rich/home/matlab/m_mapWK/pb/washington2pts.txt
/ocean/rich/more/mmapbase/bcgeo/PNW.mat
has a 2-column variable 'ncst' with NaN-separated line segments that make polygons for the BC/WA coast.
Coastline¶
In [7]:
topo_datastruct = sio.loadmat('/ocean/rich/more/mmapbase/bcgeo/PNW.mat')
coast={}
coast['lat'] = topo_datastruct['ncst'][:,1]
coast['lon'] = topo_datastruct['ncst'][:,0]
def draw_coast(ax):
#ax.set_aspect(5/4.4)
#Plot coast line
ax.plot(coast['lon'],coast['lat'],'-k',rasterized=True,markersize=1)
ax.set_xlim([-128.5, -121])
ax.set_ylim([46, 52])
#ax.set_xticks([])
#ax.set_yticks([])
Bathy¶
In [8]:
#bathymetry
fB = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc','r')
lats = fB.variables['nav_lat']
lons = fB.variables['nav_lon']
D = fB.variables['Bathymetry']
#rivers info
test = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/rivers/rivers_cnst.nc','r')
plotting = test.variables['rorunoff'][0,:,:]
def draw_bathy(ax,grid):
if grid=='grid':
mesh=ax.pcolormesh(D[:],cmap=cmap)
elif grid=='map':
mesh=ax.pcolormesh(lons[:],lats[:],D[:],cmap=cmap)
mesh=ax.contourf(lons[:],lats[:],D[:],15,cmap=cmap)
viz_tools.plot_coastline(ax,fB,coords=grid)
cbar=fig.colorbar(mesh,ax=ax)
ax.set_aspect(5/4.4)
cbar.set_label('{depth.long_name} [{depth.units}]'.format(depth=D), **{'size': '16'})
Rivers¶
In [9]:
def draw_rivers(ax,grid):
for i in range(0,898):
for j in range(0,398):
if plotting[i,j] > 0:
if grid=='grid':
ax.plot(j,i,'ob',markersize=3)
elif grid=='map':
ax.plot(lons[i,j],lats[i,j],'og',markersize=5,rasterized=True)
Storms¶
In [10]:
#Storm surge points
istorms=[328,195,213,124]
jstorms=[467,298,351,747]
ls=['Point Atkinson','Victoria','Patricia Bay', 'Campbell River']
def draw_storms(ax,grid):
for k in range(0,len(istorms)):
if grid=='grid':
ax.plot(istorms[k],jstorms[k],'*r',markersize=15)
elif grid=='map':
x=lons[jstorms[k],istorms[k]]; y=lats[jstorms[k],istorms[k]]
label=ls[k]
ax.plot(x,y,'*r',markersize=15,rasterized=True)
if label == 'Victoria':
ax.annotate(label, xy=(x,y), xytext = (-65, 8),textcoords = 'offset points',
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=11)
elif label =='Point Atkinson':
ax.annotate(label, xy=(x,y), xytext = (10, 30),textcoords = 'offset points',
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=11)
elif label =='Campbell River':
ax.annotate('Campbell \n River', xy=(x,y), xytext = (-70, -30),textcoords = 'offset points',
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=11)
if label == 'Patricia Bay':
ax.annotate(label, xy=(x,y), xytext = (-90, -5),textcoords = 'offset points',
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=11)
Tides¶
In [11]:
filename = '/data/nsoontie/MEOPAR/analysis/compare_tides/obs_tidal_wlev_const_all.csv'
harm_obs = pd.read_csv(filename,sep=';',header=0)
harm_obs = harm_obs.rename(columns={'Site': 'site', 'Lat': 'lat', 'Lon': 'lon',
'M2 amp': 'M2_amp', 'M2 phase (deg UT)': 'M2_pha',
'K1 amp': 'K1_amp', 'K1 phase (deg UT)': 'K1_pha'})
def draw_tidestations(ax):
stations_obs = ['Port Renfrew','Sheringham Point','Pedder Bay', 'Esquimalt',
'Victoria','Clover Point','Finnerty Cove', 'Fulford Harbour',
'Tumbo Channel','Patos Island','Whaler Bay', 'Tsawwassen',
'Sandheads', 'Point Grey','Point Atkinson','Gibsons Landing',
'Halfmoon Bay','Irvines Landing','Powell River', 'Lund',
'Twin Islets','Campbell River','Maude Island E', 'Nymphe Cove',
'Seymour Narrows','Brown Bay','Chatham Point','Kelsey Bay','Yorke Island']
numsta=len(stations_obs)
lon=[]; lat=[]
for stn in range(numsta):
location = stations_obs[stn]
lon.append(-harm_obs.lon[harm_obs.site==location])
lat.append(harm_obs.lat[harm_obs.site==location])
print stn+1, location
ax.annotate(stn+1, xy = (-harm_obs.lon[harm_obs.site==location],harm_obs.lat[harm_obs.site==location]),
xytext = (-5,5),ha = 'left', va = 'bottom',textcoords = 'offset points',fontsize=10)
ax.scatter(lon,lat,marker='^',s=80,edgecolor='black',linewidth='1',facecolor='green',rasterized=True)
In [12]:
cmap = plt.get_cmap('Blues')
cmap.set_bad('gainsboro')
fig, ax = plt.subplots(1, 1, figsize=(20,15))
g = 'map'
draw_coast(ax)
draw_bathy(ax,g)
ax.set_ylabel('{latitude.long_name} [{latitude.units}]'.format(latitude=lats), **{'size': '16'})
ax.set_xlabel('{longitude.long_name} [{longitude.units}]'.format(longitude=lons), **{'size': '16'})
ax.set_xlim([-126.4,-121.3])
ax.set_ylim([46.8,51.1])
ax.set_axis_bgcolor('darkgray')
ax.text(-122.6,50.2,' British\nColumbia',fontsize=20)
ax.text(-124.6,48.8,'Vancouver\n Island', fontsize=20)
ax.text(-125.8,48,'Pacific\nOcean', fontsize=20)
ax.text(-123.9,47.5,'Washington\n State', fontsize=20)
ax.text(-123.57,49,' Strait\n of\nGeorgia', fontsize=16)
ax.text(-125.9,50.1,'Johnstone \n Strait', fontsize=16)
ax.text(-124.64,48.43,'Strait of Juan de Fuca',fontsize=16,rotation=-16)
ax.text(-122.3,47.7,'Puget\nSound', fontsize=16)
ax.plot(-123.1,49.23,'sw',markersize=12,rasterized=True)
ax.text(-123.03,49.21,'Vancouver',fontsize=16)
ax.annotate('Fraser River', xy=(-123.06,49.13), xytext = (-122.8,49.13),textcoords = 'data',
arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=16)
fig.patch.set_facecolor('#8BB0CC')
fig.show()
/home/imachuca/anaconda/lib/python2.7/site-packages/matplotlib/figure.py:387: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure "matplotlib is currently using a non-GUI backend, "
In [ ]:
In [ ]: