The notebook will generate a two-dimensional bathymetry along the thalweg for simulation runs.

Two-dimensional simulations will allow us to quickly experiment with changes in resolution, mxing parameters etc

In [1]:
import netCDF4 as nc
import numpy as np
from salishsea_tools import tidetools, nc_tools
import matplotlib.pyplot as plt

%matplotlib inline

In [2]:
grid=nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
bathy,X,Y=tidetools.get_bathy_data(grid)

'/data/nsoontie/MEOPAR/tools/bathymetry/thalweg_working.txt',
delimiter=" ", unpack=False)
lines = lines.astype(int)

lines.shape

Out[2]:
(1539, 2)
In [3]:
plt.plot(-bathy[lines[:,0],lines[:,1]])

Out[3]:
[<matplotlib.lines.Line2D at 0x7f83967c5f10>]

To do:

1. Exclude northern SoG and Johnstone Strait
2. Close the right most boundary.
3. Smooth (make SoG more like a basin?)
4. Set a spatial scale (easiest to set a 500 m grid spacing between points)
5. Extend into y by 10 gridpoints

# 1. Exclude north¶

Exlcude the north by setting a cuttoff value.

In [4]:
ie=1100

thalweg_nonorth=-bathy[lines[0:ie,0],lines[0:ie,1]]
plt.plot(thalweg_nonorth)

Out[4]:
[<matplotlib.lines.Line2D at 0x7f8396646f10>]

Look for the index of the peak at the right:

In [5]:
imax=np.argmax(thalweg_nonorth[1000:ie])
thalweg_nonorth[1000:ie][imax]

Out[5]:
-169.0
In [6]:
ind=np.where(thalweg_nonorth==thalweg_nonorth[1000:ie][imax])
ind

Out[6]:
(array([ 388,  389,  418,  420,  439, 1040]),)
In [7]:
indp = ind[0][-1]
indp

Out[7]:
1040

# 2. Close right boundary¶

After the peak, extend to surface. Just use a line for now...

In [8]:
slope = (0-thalweg_nonorth[indp])/(len(thalweg_nonorth)-1-indp)
thalweg_right=thalweg_nonorth.copy()
thalweg_right[indp:] = (np.arange(indp,len(thalweg_nonorth))-(len(thalweg_nonorth)-1))*slope
plt.plot(thalweg_right)

Out[8]:
[<matplotlib.lines.Line2D at 0x7f8396590f50>]
In [9]:
thalweg_right[-1]

Out[9]:
0.0

# 3. Smooth¶

Smoothing from signal processing. Convolve with a window function. Try a common functions and check result

In [10]:
def smooth(x,window_len,window):
s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
if window == 'flat': #moving average
w=np.ones(window_len,'d')
else:
w=eval('np.'+window+'(window_len)')

y=np.convolve(w/w.sum(),s,mode='valid')
y= y[(window_len/2-1):-(window_len/2)]

return y

In [11]:
fig=plt.figure(figsize=(10,5))
w_len=30
windows=['flat','hanning','hamming', 'bartlett','blackman']
plt.plot(thalweg_right,label='original')
for w in windows:
y=smooth(thalweg_right,w_len,w)
print y.shape
plt.plot(y,label=w)

plt.legend(loc=0)

(1100,)
(1100,)
(1100,)
(1100,)
(1100,)

Out[11]:
<matplotlib.legend.Legend at 0x7f8395f2a910>

Many featurs remain present (like the SoG basin). The large peak around 600 is much smaller ...

Use blackman because it keeps the beak the highest and extends closes to surface at right.

In [12]:
w_len=30
y=smooth(thalweg_right,w_len,'blackman')
y[-1]=0;#keep right end closed
thalweg_smooth=y
plt.plot(thalweg_smooth,label='smoothed')

Out[12]:
[<matplotlib.lines.Line2D at 0x7f8395db6050>]
In [13]:
y[-1]

Out[13]:
0.0

# 4. Spatial scale¶

Set hoizontal grid spacing to 500m.

In [14]:
dx=500
# x is variable for horiztonal space
x=dx*np.linspace(0,thalweg_smooth.shape[0]-1,thalweg_smooth.shape[0]);

plt.plot(x/1000,thalweg_smooth)

Out[14]:
[<matplotlib.lines.Line2D at 0x7f8395d9d590>]

# 5. Extend into y¶

In [15]:
Ny=10; dy=500
y = dy*np.linspace(0,Ny-1,Ny)

thalweg_y=np.tile(thalweg_smooth,Ny)
thalweg_y=thalweg_y.reshape((y.shape[0],thalweg_smooth.shape[0]))
thalweg_y.shape

Out[15]:
(10, 1100)
In [16]:
plt.plot(thalweg_y[0,:])

Out[16]:
[<matplotlib.lines.Line2D at 0x7f8395c33290>]
In [17]:
xx,yy=np.meshgrid(x,y)
xx.shape

Out[17]:
(10, 1100)
In [18]:
plt.pcolormesh(xx,yy,thalweg_y)
plt.colorbar()

Out[18]:
<matplotlib.colorbar.Colorbar instance at 0x7f8395afde18>

# Finalize and Create Netcdf file¶

In [19]:
thalweg_final=-thalweg_y #positive bathy values for saving

new_bathy = nc.Dataset('../bathy2D.nc', 'w')
new_bathy.createDimension('y', thalweg_final.shape[0])
new_bathy.createDimension('x', thalweg_final.shape[1])
nc_tools.show_dimensions(new_bathy)

<type 'netCDF4.Dimension'>: name = 'y', size = 10

<type 'netCDF4.Dimension'>: name = 'x', size = 1100


In [20]:
new_x = new_bathy.createVariable('x', float, ('y', 'x'), zlib=True)
new_x.setncattr('units', 'metres')
new_y = new_bathy.createVariable('y', float, ('y', 'x'), zlib=True)
new_y.setncattr('units', 'metres')
newdepths = new_bathy.createVariable(
'Bathymetry', float, ('y', 'x'),
zlib=True, least_significant_digit=1)
newdepths.setncattr('units', 'metres')

In [21]:
new_x[:]=xx;
new_y[:]=yy;
newdepths[:]=thalweg_final

In [22]:
new_bathy.title="""Thalweg Bathymetry for 2D model"""
new_bathy.institution= """
Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia"""
new_bathy.comment= """
Based on bathy_meter_SalishSea2.nc along thalweg with smoothing"""
new_bathy.reference= """
/data/nsoontie/MEOPAR/2Ddomain/notebooks/Generate_2D_bathy.ipynb"""

In [23]:
nc_tools.show_dataset_attrs(new_bathy)

file format: NETCDF4
title: Thalweg Bathymetry for 2D model
institution:
Dept of Earth, Ocean & Atmospheric Sciences, University of British Columbia
comment:
Based on bathy_meter_SalishSea2.nc along thalweg with smoothing
reference:
/data/nsoontie/MEOPAR/2Ddomain/Generate_2D_bathy.ipynb

In [24]:
new_bathy.close()

In [25]:
B=nc.Dataset('bathy2D.nc')
thal=B.variables['Bathymetry']
plt.pcolormesh(thal[:]); plt.colorbar()

Out[25]:
<matplotlib.colorbar.Colorbar instance at 0x7f83955354d0>
In [26]:
plt.plot(-thal[0,:])

Out[26]:
[<matplotlib.lines.Line2D at 0x7f839502d0d0>]
In [28]:
thal[0,-1]

Out[28]:
-0.0

Looks gooe enough to me.

# Next Steps¶

1. Create z grid
2. Generate initial conditions
3. Generate forcings
4. Look into coordinates defined in metres, not lon/lat (jphgr_mesh=2, set ppe1_m,ppe2_m). Also consider looking at GYRE configuration.
5. Run?

# Ideas and questions¶