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
NEMO 3.6 is expecting obcs t start on a land/oean corner. This notebook will masked all of the edges.
import netCDF4 as nc
import numpy as np
from salishsea_tools import tidetools, nc_tools
import matplotlib.pyplot as plt
%matplotlib inline
grid=nc.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
bathy,X,Y=tidetools.get_bathy_data(grid)
lines = np.loadtxt(
'/data/nsoontie/MEOPAR/tools/bathymetry/thalweg_working.txt',
delimiter=" ", unpack=False)
lines = lines.astype(int)
lines.shape
(1539, 2)
plt.plot(-bathy[lines[:,0],lines[:,1]])
[<matplotlib.lines.Line2D at 0x7ff11bba4690>]
To do:
Exlcude the north by setting a cuttoff value.
ie=1100
thalweg_nonorth=-bathy[lines[0:ie,0],lines[0:ie,1]]
plt.plot(thalweg_nonorth)
[<matplotlib.lines.Line2D at 0x7ff11ba68c10>]
Look for the index of the peak at the right:
imax=np.argmax(thalweg_nonorth[1000:ie])
thalweg_nonorth[1000:ie][imax]
-169.0
ind=np.where(thalweg_nonorth==thalweg_nonorth[1000:ie][imax])
ind
(array([ 388, 389, 418, 420, 439, 1040]),)
indp = ind[0][-1]
indp
1040
After the peak, extend to surface. Just use a line for now...
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)
[<matplotlib.lines.Line2D at 0x7ff115bf4a50>]
thalweg_right[-1]
0.0
Smoothing from signal processing. Convolve with a window function. Try a common functions and check result
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
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,)
<matplotlib.legend.Legend at 0x7ff115b3a390>
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.
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')
[<matplotlib.lines.Line2D at 0x7ff115a40ed0>]
y[-1]
0.0
Set hoizontal grid spacing to 500m.
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)
[<matplotlib.lines.Line2D at 0x7ff115985590>]
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
(10, 1100)
plt.plot(thalweg_y[0,:])
[<matplotlib.lines.Line2D at 0x7ff1158c2750>]
xx,yy=np.meshgrid(x,y)
xx.shape
(10, 1100)
plt.pcolormesh(xx,yy,thalweg_y)
plt.colorbar()
<matplotlib.colorbar.Colorbar instance at 0x7ff11578ccb0>
/home/nsoontie/anaconda3/envs/py2/lib/python2.7/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison if self._edgecolors == str('face'):
thalweg_mask = np.copy(thalweg_y)
thalweg_mask[0:2,:]=0
thalweg_mask[-2:,:]=0
thalweg_mask[:,0]=0
thalweg_mask[:,-1:]=0
plt.pcolormesh(thalweg_mask,vmin=-50,vmax=0)
plt.colorbar()
plt.axis([0,10,0,10])
[0, 10, 0, 10]
thalweg_final=-thalweg_mask #positive bathy values for saving
new_bathy = nc.Dataset('../grid/bathy2D_36.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._netCDF4.Dimension'>: name = 'y', size = 10 <type 'netCDF4._netCDF4.Dimension'>: name = 'x', size = 1100
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')
new_x[:]=xx;
new_y[:]=yy;
newdepths[:]=thalweg_final
new_bathy.title="""Thalweg Bathymetry for 2D model NEMO 3.6"""
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"""
nc_tools.show_dataset_attrs(new_bathy)
file format: NETCDF4 title: Thalweg Bathymetry for 2D model NEMO 3.6 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/notebooks/Generate_2D_bathy.ipynb
new_bathy.close()
B=nc.Dataset('../grid/bathy2D_36.nc')
thal=B.variables['Bathymetry']
plt.pcolormesh(thal[:]); plt.colorbar()
plt.axis([0,10,0,10])
[0, 10, 0, 10]
plt.plot(-thal[2,:])
[<matplotlib.lines.Line2D at 0x7ff11540ef50>]
thal[2,-1]
-0.0
Looks good enough to me.