%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import matplotlib as mpl
import netCDF4 as NC
import numpy as np
import itertools
The following function will compute the grid points of the thalweg beginning at (ist,jst) and ending at ijend. The thalweg is the path that follows the maximum depth through the Strait. This function works by finding the maximum depth in a neighbourhood of (ist, jst) and then iterating through i or j (set by the direction flag) and terminates when the iteration reches ijend. At each iteration it records the location of the maximum depth in a neighbourhood of (i,j). The size of the neighbourhood is set by thres.
The resulting thalweg is a little bit choppy. This function must be applied several times to sections of the Strait.
def find_thalweg(bathy,ist,jst,direction,ijend,thres):
#Starting at (ist,jst) this function computes the gridpoints of thalweg by iterating through points until ijend.
#At each iteration it calculates the maximum depth in a neighborhood of (i,j)
#bathy: an array of the bathymetry
#ist: an int that defines the starting i grid point
#jst: an int that defines the starting j grid point
#direction: a string either 'i' or 'j' indicating the direction of the iterations
#thres: an int defining the size of the neighbourhood.
#returns thalweg: an array contatining the [j,i] thalweg gridpoints.
if direction is 'i':
#thres defines the search distance in the i-direction.
if ist<ijend:
thalweg=np.zeros((ijend-ist,2))
loop_start=ist+1; loop_end=ijend
forward=1; backward=0
else:
thalweg=np.zeros((ist-ijend,2))
loop_start=0; loop_end=ist-ijend-1;
forward=0; backward=1;
#find index of first maximum
bathmax=bathy[jst-thres:jst+thres,ist].max();
ind= np.where(bathy[:,ist]==bathmax)
ls=0
while ind[0][ls]< jst-thres:
ls=ls+1
j=ind[0][ls]
arr_ind=0*forward + (ist-ijend-1)*backward;
thalweg[arr_ind,0]=j
thalweg[arr_ind,1]=ist
#for loop for finding the thalweg
for k in range (loop_start,loop_end):
loop_ind = k*forward + (ist-k)*backward;
bathmax=bathy[j-thres:j+thres,loop_ind].max();
ind= np.where(bathy[:,loop_ind]==bathmax)
ls=0
while ind[0][ls]< j-thres:
ls=ls+1
j=ind[0][ls]
arr_ind=(k-ist)*forward + (loop_end-k-1)*backward;
thalweg[arr_ind,0]=j
thalweg[arr_ind,1]=loop_ind;
elif direction is 'j':
#thres defines the search distance in the j-direction.
if jst<ijend:
thalweg=np.zeros((ijend-jst,2))
loop_start=jst+1; loop_end=ijend
forward=1; backward=0
else:
thalweg=np.zeros((jst-ijend,2))
loop_start=0; loop_end=jst-ijend-1;
forward=0; backward=1;
#find index of first maximum
bathmax=bathy[jst,ist-thres:ist+thres].max();
ind= np.where(bathy[jst,:]==bathmax)
ls=0
while ind[0][ls]< ist-thres:
ls=ls+1
i=ind[0][ls]
arr_ind=0*forward + (jst-ijend-1)*backward;
thalweg[arr_ind,0]=jst
thalweg[arr_ind,1]=i
#for loop for finding the thalweg
for k in range (loop_start,loop_end):
loop_ind = k*forward + (jst-k)*backward;
bathmax=bathy[loop_ind,i-thres:i+thres].max();
ind= np.where(bathy[loop_ind,:]==bathmax)
ls=0
while ind[0][ls]< i-thres:
ls=ls+1
i=ind[0][ls]
arr_ind=(k-jst)*forward + (loop_end-k-1)*backward;
thalweg[arr_ind,0]=loop_ind
thalweg[arr_ind,1]=i;
else:
return 'Error'
return thalweg
def remove_section(thalweg, start, end, direction):
#This function is used to remove a jagged section of the thalweg and replace it with a straight line.
#thalweg: an array containing the thalweg grid coordinates
#start: an int containing the i or j coordintate where the cut begins
#end: an int contatining the i o j coordinate where the cut ends
#direction: indicates where start and end represent an i or j coordinate.
if direction is 'j':
flag1=0; flag2=1
elif direction is 'i':
flag1=1; flag2=0
else:
return 'Error'
#find index of cut points in other coordinate
index= np.where(thalweg[:,flag1]==start); index_st=index[0][0]
index= np.where(thalweg[:,flag1]==end); index_end=index[0][0]
#start and end points of
start2=thalweg[index_st,flag2]; end2=thalweg[index_end,flag2];
#draw a straight line to connect (start,start2) to (end,end2)
newinds = start2 + (end2-start2)/(end-start)*(thalweg[index_st:index_end,flag1]-start)
newinds=np.rint(newinds)
#replace the thalwegs with newinds
thalweg_copy=thalweg
thalweg_copy[index_st:index_end,flag2] = newinds
return thalweg_copy
def distance(thalweg,dx,dy):
#this function calculates the distance travelled along the thalweg
#thalweg is the array of grid coordinates.
#dx is the grid spacing in x (eventually we could send an array in but for now it is a single value).
#dy is the grid spacing in y (eventually we could send an array in but for now it is a single value).
leng=thalweg.shape
l =leng[0]
distance = np.zeros(l)
d=0
for k in range(1,l):
i1=thalweg[k-1,1]; i2=thalweg[k,1];
j1=thalweg[k-1,0]; j2=thalweg[k,0];
dseg =( (dx*(i2-i1))**2 + (dy*(j2-j1))**2)**(0.5)
d=d+dseg
distance[k]=d;
return distance
# test the function
grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea.nc','r')
bathy = grid.variables['Bathymetry'][:,:]
#masking
mu =bathy == 0
bathy_mask= np.ma.array(bathy,mask=mu)
#the thalweg is caclualted in sections:
#caculate first section of thalweg
step=4; ist=264; jst=415; ijend=700; direction = 'j';
thalweg5 = find_thalweg(bathy,ist,jst,direction,ijend,step)
#remove several jagged sections
jstart=600; jend=645;
thalweg5 = remove_section(thalweg5,jstart,jend,'j');
jstart=425; jend=447;
thalweg5 = remove_section(thalweg5,jstart,jend,'j');
#next section
jst=int(thalweg5[0,0]); ist= int(thalweg5[0,1]); ijend=346; direction='j';
thalweg4 = find_thalweg(bathy,ist,jst,direction,ijend,step)
#next section
jst=int(thalweg4[0,0]); ist= int(thalweg4[0,1]); ijend=235; direction='i';
thalweg3 = find_thalweg(bathy,ist,jst,direction,ijend,step)
#next section
jst=int(thalweg3[0,0]); ist= int(thalweg3[0,1]); ijend=246; direction='j';
thalweg2 = find_thalweg(bathy,ist,jst,direction,ijend,step)
#next section
jst=int(thalweg2[0,0]); ist= int(thalweg2[0,1]); ijend=0; direction='i';
thalweg1 = find_thalweg(bathy,ist,jst,direction,ijend,step)
#Combine all together
Tjs = [thalweg1[:,0],thalweg2[:,0],thalweg3[:,0],thalweg4[:,0],thalweg5[:,0]];
Tjs =list(itertools.chain.from_iterable(Tjs)); Tjs=np.array(Tjs);
Tis = [thalweg1[:,1],thalweg2[:,1],thalweg3[:,1],thalweg4[:,1],thalweg5[:,1]];
Tis =list(itertools.chain.from_iterable(Tis)); Tis=np.array(Tis);
thalweg=[Tjs, Tis]; thalweg=zip(*thalweg); thalweg=np.array(thalweg,'i');
plt.figure(figsize=(8,10))
cmin=0; cmax=400;
pylab.pcolormesh(bathy,vmin=cmin,vmax=cmax)
pylab.colorbar()
pylab.plot(thalweg[:,1]+1,thalweg[:,0]+1,marker='*',color='k')
for k in range (1,7):
pylab.plot(thalweg[k*100,1]+1, thalweg[k*100,0]+1,'o')
print thalweg[k*100,1] + 1, thalweg[k*100,0] + 1
#pylab.axis([240,300,300,400])
103 321 203 251 240 330 273 385 247 483 162 583
#look at the depths along the thalweg
ax=[0,500,-450,0]
depths = bathy[thalweg[:,0],thalweg[:,1]];
dx=400; dy=500
d = distance(thalweg,dx,dy); d=d/1000;
plt.figure(figsize=(14,6))
pylab.plot(d,-depths[:],color='b')
for k in range (1,7):
pylab.plot(d[k*100], -depths[k*100],'o')
pylab.title('Depth along thalweg')
pylab.axis(ax)
pylab.xlabel('km')
pylab.ylabel(' z (m)')
<matplotlib.text.Text at 0x3e34890>
plt.figure(figsize=(14,6))
pylab.plot(-depths[:],color='b')
for k in range (1,7):
pylab.plot(k*100, -depths[k*100],'o')
#plot a salinity field along the thalweg
f = NC.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/bottom/71d100d_bfri3e-3/output.abort.nc','r');
sal = f.variables['vosaline']
salPlot=sal[0,:,:,:];
smin=28; smax=34
dep = f.variables['deptht']
vs=np.arange(34,27.5,-0.5);
salThalweg=salPlot[:,thalweg[:,0],thalweg[:,1]];
#masking
mu =salThalweg == 0
salThalweg= np.ma.array(salThalweg,mask=mu)
XX,ZZ = np.meshgrid(d,-dep[:])
plt.figure(figsize=(19,6))
pylab.pcolormesh(XX,ZZ,salThalweg,vmin=smin,vmax=smax,cmap='rainbow')
pylab.colorbar()
CS=pylab.contour(XX,ZZ,salThalweg,vs, colors='black')
pylab.clabel(CS, fontsize=9, inline=1)
for k in range (1,7):
pylab.plot(d[k*100], -depths[k*100],'o')
#pylab.axis(ax)
pylab.title('Salinity Field along Thalweg')
pylab.axis(ax)
pylab.xlabel('km')
pylab.ylabel(' z (m)')
<matplotlib.text.Text at 0x3dde5d0>
#plot a salinity field along the thalweg
f = NC.Dataset('/ocean/dlatorne/MEOPAR/SalishSea/results/61d70d_nu200/SalishSea_3d_20021114_20021123_grid_T.nc','r');
sal = f.variables['vosaline']
salPlot=sal[0,:,:,:];
smin=28; smax=34
dep = f.variables['deptht']
vs=np.arange(34,27.5,-0.5);
salThalweg=salPlot[:,thalweg[:,0],thalweg[:,1]];
#masking
mu =salThalweg == 0
salThalweg= np.ma.array(salThalweg,mask=mu)
XX,ZZ = np.meshgrid(d,-dep[:])
plt.figure(figsize=(19,6))
pylab.pcolormesh(XX,ZZ,salThalweg,vmin=smin,vmax=smax,cmap='rainbow')
pylab.colorbar()
CS=pylab.contour(XX,ZZ,salThalweg,vs, colors='black')
pylab.clabel(CS, fontsize=9, inline=1)
for k in range (1,7):
pylab.plot(d[k*100], -depths[k*100],'o')
#pylab.axis(ax)
pylab.title('Salinity Field along Thalweg')
pylab.axis(ax)
pylab.xlabel('km')
pylab.ylabel(' z (m)')
<matplotlib.text.Text at 0x4222c90>
Thalweg index points are saved in: /data/nsoontie/MEOPAR/tools/analysis_tools/thalweg.txt. First column corresponds to j index, second column corresponds to i index (in python indexing convention).
This block of code is an example of how to load the data and plot:
#load thalweg points
lines = np.loadtxt('thalweg.txt', delimiter=" ", unpack=False)
lines = lines.astype(int)
#load data to plot
f = NC.Dataset('/data/nsoontie/MEOPAR/SalishSea/results/bottom/71d100d_bfri3e-3/output.abort.nc','r');
sal = f.variables['vosaline']
T_lon=f.variables['nav_lon']
T_lat=f.variables['nav_lat']
dep = f.variables['deptht']
#salinity, lon and lat along thalweg
salThal=sal[0,:,lines[:,0],lines[:,1]];
thalweg_lon = T_lon[lines[:,0],lines[:,1]]
thalweg_lat = T_lat[lines[:,0],lines[:,1]]
#array of thalweg indices
ds=np.arange(0,lines.shape[0],1);
#plotting parameters
smin=28; smax=34;
vs=np.arange(34,27.5,-0.5);
#masking
mu =salThal == 0; salThal= np.ma.array(salThal,mask=mu)
XX,ZZ = np.meshgrid(ds,-dep[:])
plt.figure(figsize=(19,12))
pylab.subplot(2,1,1)
pylab.pcolormesh(XX,ZZ,salThal,vmin=smin,vmax=smax,cmap='rainbow')
pylab.colorbar()
CS=pylab.contour(XX,ZZ,salThal,vs, colors='black')
pylab.clabel(CS, fontsize=9, inline=1)
pylab.title('Salinity along Thalweg')
<matplotlib.text.Text at 0x5ba3410>