In [2]:

from __future__ import division, print_function

import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np

from salishsea_tools import (
    nc_tools,
    viz_tools,
    tidetools,
)

In [3]:

%matplotlib inline

In [4]:

grid = nc.Dataset('/data/dlatorne/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')

In [5]:

bathy = grid.variables['Bathymetry'][:]
lats = grid.variables['nav_lat'][:]
lons = grid.variables['nav_lon'][:]

Clover Point 48°24'09.8"N 123°20'55.6"W,¶

In [58]:

print (48+24/60.+09.8/60./60.)
dep = 65

48.4027222222

In [59]:

yind, xind = tidetools.find_closest_model_point(-123.348778, 48.402722,lons,lats,bathy, lon_tol = 0.082, lat_tol = 0.003) 
level = tidetools.find_model_level(dep, zlevels, fractional=False)

returns the y-index, yind (latitude) and x-index, xind (longitude) of the grid from specified lon/lat

In [62]:

print (yind, xind, level)
print(bathy[yind-1,xind])

289 210 24
70.0

y-index and x-index

In [63]:

yind = yind - 1
# Look at cell and variation
print (lats[yind,xind], lats[yind+1,xind], lats[yind-1,xind])
print (lons[yind, xind], lons[yind, xind+1], lons[yind, xind-1])

48.3971595764 48.4011230469 48.3931999207
-123.265525818 -123.260314941 -123.270736694

McCaulay 48°25'05.2"N 123°24'22.5"W 60m¶

In [91]:

lon = -(123+24/60.+22.5/60./60.)
lat = 48+25/60.+05.2/60./60.
dep = 60.

In [92]:

yind, xind = tidetools.find_closest_model_point(lon, lat,  lons,lats,bathy) 
level = tidetools.find_model_level(60, zlevels, fractional=False)
print (yind, xind, level)
print (bathy[yind-3,xind-2])

302 192 24
56.0

In [93]:

yind = yind - 3; xind = xind-2
# Look at cell and variation
print (lats[yind,xind], lats[yind+1,xind], lats[yind-1,xind])
print (lons[yind, xind], lons[yind, xind+1], lons[yind, xind-1])

48.4022827148 48.4062423706 48.3983192444
-123.406349182 -123.401138306 -123.411552429

Iona 49°12'20.6"N 123°15'42.8"W, 7.5 km offshore; depth of 90 m¶

In [89]:

lon = -(123+15/60.+42.8/60./60.)
lat = 49+12/60.+20.6/60./60.
dep = 90.
yind, xind = tidetools.find_closest_model_point(lon, lat,  lons,lats,bathy) 
level = tidetools.find_model_level(dep, zlevels, fractional=False)
print (yind, xind, level)
print (bathy[yind,xind-9])

444 311 26
103.0

Need to move 3 m up in depth and 9 grid points out from coast to get out of the topography

In [90]:

xind = xind-9
# Look at cell and variation
print (lats[yind,xind], lats[yind+1,xind], lats[yind-1,xind])
print (lons[yind, xind], lons[yind, xind+1], lons[yind, xind-1])

49.1873512268 49.1912727356 49.1834335327
-123.308982849 -123.303787231 -123.314178467

Nanaimo 49°13'50.0"N 123°58'04.5"W, 2.5 km offshore; 70m below sea level¶

In [47]:

lon = -(123+58/60.+04.5/60./60.)
lat = 49+13/60.+50.0/60./60.
dep = 70.
yind, xind = tidetools.find_closest_model_point(lon, lat,  lons,lats,bathy) 
level = tidetools.find_model_level(dep, zlevels, fractional=False)
print (yind, xind, level)
print (bathy[yind,xind+2])

500 211 25
83.8125

In [48]:

xind = xind+2
# Look at cell and variation
print (lats[yind,xind], lats[yind+1,xind], lats[yind-1,xind])
print (lons[yind, xind], lons[yind, xind+1], lons[yind, xind-1])

49.2388648987 49.2427902222 49.2349433899
-123.958374023 -123.953186035 -123.963562012

Campbell River 50°03'25.6"N 125°15'48.3"W, 540m offshore; minimum depth 35m¶

In [51]:

yind, xind = tidetools.find_closest_model_point(-125.263417, 50.057111,  lons,lats,bathy) 

returns the y-index, yind (latitude) and x-index, xind (longitude) of the grid from specified lon/lat

In [54]:

print (yind, xind)
print (bathy[yind,xind+1])

751 124
42.0

In [55]:

tidetools.find_model_level(35, zlevels, fractional=False)

Out[55]:

In [56]:

xind=xind+1
# Look at cell and variation
print (lats[yind,xind], lats[yind+1,xind], lats[yind-1,xind])
print (lons[yind, xind], lons[yind, xind+1], lons[yind, xind-1])

50.055103302 50.059009552 50.051197052
-125.250205994 -125.244987488 -125.2554245

Index of the model closest to the specified depth

In [18]:

v_vel = nc.Dataset('/ocean/dlatorne/MEOPAR/SalishSea/results/spin-up/17dec26dec/SalishSea_1d_20031217_20031226_grid_V.nc')
zlevels = v_vel.variables['depthv'][:]