Salish Sea NEMO Sub-domain Bathymetry¶

This notebook documents the bathymetry used for the initial Salish Sea NEMO runs on a sub-set of the whole region domain. This sub-domain was used for the runs known as JPP and WCSD_RUN_tide_M2_OW_ON_file_DAMP_ANALY.

The first part of the notebook explores and plots the bathymetry. The second part records the manual smoothing that was done to get the JPP M2 tidal forcing case to run.

In [1]:

%matplotlib inline
import netCDF4 as nc
import numpy as np

bathy_tools is a collection of Python functions that were developed to encapsulate the details of working with a bathymetry netCDF file.

In [2]:

from salishsea_tools import bathy_tools

Open the file and get some basic information about it. The r mode shown opens the file in read-only mode. Use r+ as the mode if you want to change the bathymetry.

In [3]:

bathy = nc.Dataset('../../NEMO-forcing/grid/SubDom_bathy_meter_NOBCchancomp.nc', 'r')
print bathy.file_format
bathy_tools.show_dimensions(bathy)
bathy_tools.show_variables(bathy)

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

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

[u'nav_lon', u'nav_lat', u'Bathymetry']

Assign short names to the variables and get some basic information about their values:

In [4]:

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

In [5]:

lat_stats = bathy_tools.min_mid_max(lats)
lon_stats = bathy_tools.min_mid_max(lons)
print lat_stats
print lon_stats
print np.max(depths)

(47.630783, 48.702590942382812, 49.751278)
(-125.20717, -123.60331726074219, -121.98944)
428.0

Plot a colour-mesh plot of the whole dataset. Use fig.savefig('filename.png') to store the plot as an image file in the pwd.

In [6]:

fig = bathy_tools.plot_colourmesh(
    bathy, 'Salish Sea Sub-domain Bathymetry', 
    axis_limits=(-125, -122.4, 48.08, 49.71))

Manual Smoothing¶

Manual smoothing of the sub-domain bathymetry was done iteratively via the following steps:

Run the JPP case and note the coordinates

at which the run aborted, typically due to high u-velocity, but occassionally due to negative surface salinity. 2. Susan examined the bathymetry in the region of the bad results point via:

bathy_tools.plot_colourmesh_zoom(bathy, centre=(fortran_i-1, fortran_j-1))

and decided on the bathymetry grid points to smooth and their new depth values. 3. Create a list of 3-tuples of the grid point to change and their new depths. 4. Print the present depths, grid of depth values in the region of the changes, and a zoomed colour-mesh plot of the region of the changes to verify that the changes will be in the correct locations. 5. Apply and verify the changes. 6. Verify that each set of previous changes are still present in the bathymetry. 7. Save and commit the updated bathymetry. 8. Repeat the cycle.

Below are the cummulative results of several passes through the above iteration cycle:

Centre of grid area to view, and coordinates of grid cell(s) to change (0-based indexing):

In [7]:

centre = (212, 93)
changes = [
    (211, 94, 14),
    (211, 92, 24),
    (212, 92, 23),
    (213, 92, 23),
]

Print the present depths, grid of depth values in the region of the changes, and a zoomed colour-mesh plot of the region of the changes to verify that the changes will be in the correct locations.

In [9]:

old_depths = (depths[jchg, ichg] for ichg, jchg, new_depth in changes)
print ' '.join(map(str, old_depths))

bathy_tools.show_region_depths(depths, centre)

bathy_tools.plot_colourmesh_zoom(bathy, centre)

14.0 24.0 23.0 23.0
[[-- -- -- 6.0 7.0 11.0 11.0 21.0 27.0 21.0]
 [-- 4.0 6.0 15.0 18.0 23.0 23.0 29.0 34.0 28.0]
 [-- 6.0 12.0 20.0 23.0 23.0 30.0 33.0 28.0 15.0]
 [-- 9.0 16.0 15.0 15.0 14.0 23.0 28.0 22.0 13.0]
 [4.0 14.0 14.0 14.0 14.0 15.0 21.0 19.0 17.0 8.0]
 [12.0 22.0 20.0 20.0 20.0 20.0 20.0 17.0 11.0 7.0]
 [24.0 25.0 24.0 26.0 24.0 23.0 23.0 23.0 12.0 4.0]
 [34.0 33.0 32.0 30.0 29.0 30.0 29.0 27.0 24.0 16.0]
 [46.0 50.0 49.0 49.0 41.0 35.0 34.0 32.0 36.0 34.0]
 [54.0 61.0 66.0 70.0 62.0 53.0 49.0 46.0 49.0 43.0]]

Apply changes:

In [365]:

for chg in changes:
    ichg, jchg, new_depth = chg
    depths[jchg, ichg] = new_depth

Confirm that Chatham Islands west fixes (94,211):13->14, (92,211):24->26, and (92,212&213):23->21&22 are present.

In [10]:

assert depths[94, 211] == 14
assert depths[92, 211] == 24
assert depths[92, 212] == 23
assert depths[92, 213] == 23
bathy_tools.show_region_depths(depths, centre)

[[-- -- -- 6.0 7.0 11.0 11.0 21.0 27.0 21.0]
 [-- 4.0 6.0 15.0 18.0 23.0 23.0 29.0 34.0 28.0]
 [-- 6.0 12.0 20.0 23.0 23.0 30.0 33.0 28.0 15.0]
 [-- 9.0 16.0 15.0 15.0 14.0 23.0 28.0 22.0 13.0]
 [4.0 14.0 14.0 14.0 14.0 15.0 21.0 19.0 17.0 8.0]
 [12.0 22.0 20.0 20.0 20.0 20.0 20.0 17.0 11.0 7.0]
 [24.0 25.0 24.0 26.0 24.0 23.0 23.0 23.0 12.0 4.0]
 [34.0 33.0 32.0 30.0 29.0 30.0 29.0 27.0 24.0 16.0]
 [46.0 50.0 49.0 49.0 41.0 35.0 34.0 32.0 36.0 34.0]
 [54.0 61.0 66.0 70.0 62.0 53.0 49.0 46.0 49.0 43.0]]

Confirm that Gordon Head fix (106,216):13->8 is present.

In [11]:

assert depths[106, 216] == 8
#print depths[111:101:-1, 211:221]

Confirm that Chatham Islands east fixes (93,219):12->9 and (93,220):22->20 are present.

In [12]:

assert depths[93, 220] == 20
assert depths[93, 219] == 9
#print depths[98:88:-1, 214:224]

Confirm that Constance Hole fixes (67&68, 215&216&217):120&121->119 are present.

In [13]:

assert depths[67, 215] == 119
assert depths[67, 216] == 119
assert depths[67, 217] == 119
assert depths[68, 215] == 119
assert depths[68, 216] == 119
assert depths[68, 217] == 119
#print depths[72:62:-1, 211:221]

Confirm that 6th Chatham Pass fixes (94,217&218):4->5, (95,217&218):4->6, (96,218):4->7 are present.

In [14]:

assert depths[94, 217] == 5
assert depths[94, 218] == 5
assert depths[95, 217] == 6
assert depths[95, 218] == 6
assert depths[96, 218] == 7
#print depths[100:90:-1, 212:222].data

Confirm that 5th Chatham Pass fix (98,217):5->8 and is present.

In [15]:

assert depths[98, 217] == 8
#print depths[103:93:-1, 212:222].data

Confirm that 4th Chatham Pass fixes (94,219):11->9 and (94,219):23->20 are present.

In [16]:

assert depths[94, 219] == 9
assert depths[94, 220] == 20
#print depths[99:89:-1, 214:224].data

Confirm that 3rd Chatham Pass fix (97,219):27->20 and (96,219):15->12 is present.

In [17]:

assert depths[97, 219] == 20
assert depths[96, 219] == 12
#print depths[101:91:-1, 214:224].data

Confirm that 2nd Chatham Pass fix (96,217)5->9 is present.

In [18]:

assert depths[96, 217] == 9
#print depths[101:91:-1, 212:222].data

Confirm that Poulier pass fix (230,239):4->12 is present.

In [19]:

assert depths[230, 239] == 12
#print depths[235:225:-1, 234:244].data

Confirm that 1st Chatham Pass fix (96,216):12->15 is present.

In [20]:

assert depths[96, 216] == 15
#print depths[101:91:-1, 211:221].data

Save the modified bathymetry:

In [5]:

bathy.close()

Algorithmic Smoothing¶

Susan developed a smoothing algorithm that successively adjusts the depths at adjacent grid point for pairs of grid points that have the greatest normalized difference in depth. The normalization factor is the average depth at the pair of grid points. The pair by pair adjustment of grid points is repeated until the maximum normalized depth difference has dropped below a threshold.

The bathy_tools.smooth() function used below is the final implementation of the algorithm. Its development can be seen in TowardSmoothing.ipynb notebook.

Bathymetry file to modify:

In [25]:

bathy = nc.Dataset('../../NEMO-forcing/grid/Subdom_bathy_meter_smooth.nc', 'r')
depths = bathy.variables['Bathymetry']

Smooth the bathymetry by successively adjusting depths that have the greatest normalized difference until all of those differences are below a threshhold. This takes several minutes...

In [14]:

depths = bathy_tools.smooth(depths, max_norm_depth_diff=0.8, smooth_factor=0.2)

Confirm that the threshold has been achieved.

In [16]:

diffs_lat = bathy_tools.calc_norm_depth_diffs(depths, delta_lat=1, delta_lon=0)
max_lat_diff_ij = bathy_tools.argmax(diffs_lat)
diffs_lon = bathy_tools.calc_norm_depth_diffs(depths, delta_lat=0, delta_lon=1)
max_lon_diff_ij = bathy_tools.argmax(diffs_lon)    

print diffs_lat[max_lat_diff_ij], diffs_lon[max_lon_diff_ij]

0.799839 0.799263

Compare the smoothed bathymetry to the original.

In [26]:

orig_bathy = nc.Dataset('../../NEMO-forcing/grid/Subdom_bathy_meter_orig.nc', 'r')
orig_depths = orig_bathy.variables['Bathymetry']

In [27]:

import matplotlib.pyplot as plt
plt.figure(figsize=(9, 9))
diffs = orig_depths[:] - depths[:]
plt.pcolormesh(-diffs[100:200,150:250], cmap=plt.cm.bwr)
plt.colorbar()

Out[27]:

<matplotlib.colorbar.Colorbar instance at 0x1058bf0e0>

In [28]:

np.min(diffs), np.max(diffs)

Out[28]:

(-98.399994, 98.399994)

In [24]:

bathy.close()

In [ ]: