from __future__ import division
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
from salishsea_tools import nc_tools
%matplotlib inline
def find_points(flow):
for i in range(390,435):
for j in range(280,398):
if flow1[0,i,j] > 0:
print i,j, lat[i,j], lon[i,j], flow[0,i,j]
grid = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
lat = grid.variables['nav_lat'][:,:]
lon = grid.variables['nav_lon'][:,:]
depth = grid.variables['Bathymetry'][:]
river1 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m05d15.nc')
river2 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m07d01.nc')
river3 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m07d02.nc')
river4 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m07d03.nc')
river5 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m07d04.nc')
river6 = nc.Dataset('/ocean/sallen/allen/research/MEOPAR/Rivers/RFraserCElse_y2015m07d05.nc')
print 'May 15'
find_points(river1.variables['rorunoff'][:,:,:])
print 'Jul 1'
find_points(river2.variables['rorunoff'][:,:,:])
print 'Jul 2'
find_points(river3.variables['rorunoff'][:,:,:])
print 'Jul 3'
find_points(river4.variables['rorunoff'][:,:,:])
print 'Jul 4'
find_points(river5.variables['rorunoff'][:,:,:])
print 'Jul 5'
find_points(river6.variables['rorunoff'][:,:,:])
May 15 411 324 49.0993995667 -123.083885193 0.72367 412 324 49.1033248901 -123.087242126 0.72367 414 334 49.1300430298 -123.041992188 7.24504 415 334 49.1339645386 -123.045349121 7.24504 416 334 49.1378898621 -123.048713684 7.24504 434 318 49.1783180237 -123.192321777 0.726198 Jul 1 411 324 49.0993995667 -123.083885193 0.674936 412 324 49.1033248901 -123.087242126 0.674936 414 334 49.1300430298 -123.041992188 8.04579 415 334 49.1339645386 -123.045349121 8.04579 416 334 49.1378898621 -123.048713684 8.04579 434 318 49.1783180237 -123.192321777 0.677294 Jul 2 411 324 49.0993995667 -123.083885193 0.674921 412 324 49.1033248901 -123.087242126 0.674921 414 334 49.1300430298 -123.041992188 7.98521 415 334 49.1339645386 -123.045349121 7.98521 416 334 49.1378898621 -123.048713684 7.98521 434 318 49.1783180237 -123.192321777 0.677279 Jul 3 411 324 49.0993995667 -123.083885193 0.676803 412 324 49.1033248901 -123.087242126 0.676803 414 334 49.1300430298 -123.041992188 7.92463 415 334 49.1339645386 -123.045349121 7.92463 416 334 49.1378898621 -123.048713684 7.92463 434 318 49.1783180237 -123.192321777 0.679167 Jul 4 411 324 49.0993995667 -123.083885193 0.681227 412 324 49.1033248901 -123.087242126 0.681227 414 334 49.1300430298 -123.041992188 7.86405 415 334 49.1339645386 -123.045349121 7.86405 416 334 49.1378898621 -123.048713684 7.86405 434 318 49.1783180237 -123.192321777 0.683607 Jul 5 411 324 49.0993995667 -123.083885193 0.674281 412 324 49.1033248901 -123.087242126 0.674281 414 334 49.1300430298 -123.041992188 7.80347 415 334 49.1339645386 -123.045349121 7.80347 416 334 49.1378898621 -123.048713684 7.80347 434 318 49.1783180237 -123.192321777 0.676636
ik = 425; jk = 302; d = 6
fig, ax = plt.subplots(1,1,figsize=(15,7.5))
imin = 390; imax = 435; jmin = 280; jmax = 398
cmap = plt.get_cmap('winter_r')
cmap.set_bad('burlywood')
mesh = ax.pcolormesh(depth[imin:imax,jmin:jmax], vmax = 10., cmap=cmap)
ax.set_xlim((0,110))
ax.set_xlabel('Grid Points')
ax.set_ylabel('Grid Points')
ax.text(40, 28, "Short Fraser River", fontsize=14)
cbar=fig.colorbar(mesh)
cbar.set_label('Depth (m)')
ax.plot(np.array((324,324,334,334,334,318))-jmin+0.5,np.array((411,412,414,415,416,434))-imin+0.5,'ko');