Site Bathymetry

The final bathymetry for Samoa has been obtained from the combination of the following bathymetries with its own resolutions:

• 50m grid provided by NIWA covering Upolu and Savaii islands and obtained as explained here
  

• 250m grid provided by NIWA covering a larger area surrounding those islands
  
  The data was interpolated to the same 0.0005 degree resolution, which corresponds to 50 m. Finally, the resultant combined bathymetry was processed to smoothen the boundaries in the transitions between the different datasets.

# common
import warnings
warnings.filterwarnings('ignore')
import os
import os.path as op
import sys
import pickle as pk

# pip
import xarray as xr
import numpy as np

# DEV: bluemath
sys.path.insert(0, op.join(op.abspath(''), '..', '..', '..', '..'))

# bluemath modules
from bluemath.plotting.nearshore import Plot_bathymetry

Database and site parameters

# database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
site='Tongatapu'
p_site = op.join(p_data, site)

p_bathy = op.join(p_site, 'bathymetry')

# bathymetry file
p_bathyfile = op.join(p_bathy, 'tongatapu_smooth.nc')


# output files: SWAN grid depth and parameters
p_grid_depth = op.join(p_bathy, 'grid200_depth.pkl')
p_grid_param = op.join(p_bathy, 'grid200_params.pkl')

Load original bathymetry

bathy = xr.open_dataset(p_bathyfile)

# sort y dimension
bathy = bathy.sortby(bathy.y, ascending=True)
print(bathy)

# plot bathymetry
Plot_bathymetry(
    bathy, 
    min_z = -200, max_z=15, 
    figsize = [25, 12],
);

<xarray.Dataset>
Dimensions:  (x: 1196, y: 1648)
Coordinates:
  * y        (y) float64 -21.75 -21.75 -21.75 -21.75 ... -20.1 -20.1 -20.1
  * x        (x) float64 184.5 184.5 184.5 184.5 ... 185.7 185.7 185.7 185.7
Data variables:
    z        (y, x) float64 ...

Define region and resolution

This is the region that we will use to run the SWAN model both for the waves and the wind

# Cutting area
region= [-175.65+360, -174.71+360,  -21.58, -20.6] #Lon_min, Lon_max, Lat_min, Lat_max

bathy_cut = bathy.sel(
    x = slice(region[0], region[1]), 
    y = slice(region[2], region[3]),
)
print(bathy_cut)

<xarray.Dataset>
Dimensions:  (x: 792, y: 980)
Coordinates:
  * y        (y) float64 -21.58 -21.58 -21.58 -21.58 ... -20.6 -20.6 -20.6 -20.6
  * x        (x) float64 184.5 184.5 184.5 184.5 ... 185.3 185.3 185.3 185.3
Data variables:
    z        (y, x) float64 605.4 605.4 605.4 605.4 ... 27.55 27.36 27.18 27.26

Interpolate Grid

from scipy import interpolate

def Interpolate_bathymetry(bati_cut, resolution):
    
    '''
    bati_cut :   original bathymetry
    resolution:  resolution for the interpolation
    '''
    
    xn=np.arange(np.min(bati_cut.x), np.nanmax(bati_cut.x), resolution)
    yn=np.arange(np.min(bati_cut.y), np.nanmax(bati_cut.y), resolution)
    
    f = interpolate.interp2d(bati_cut.x, bati_cut.y, bati_cut.z, kind='linear')
    znew = f(xn,yn)
    bati_interp = xr.Dataset(
        {
            'z': (['y','x'], znew),
        },
        coords={
            'y': yn,
            'x': xn,
        },
    )
    
    return bati_interp

resolution=0.002
bathy_interp=Interpolate_bathymetry(bathy_cut, resolution)

Plotting

# plot bathymetry
Plot_bathymetry(
    bathy_interp, 
    min_z = -400, max_z=15, 
    figsize = [25, 12],
);

SWAN Pre-processing

Here we calculate the bathymetric parameters needed to run the SWAN cases

• [xpc, ypc] : Location of the lower left corner

• [xlenc, ylenc] : Area length in x and y

• [mxc, myc] : Number of grids in x and y

• [res_x, res_y] : Resolution in x and y

# Origin
xpc = np.nanmin(bathy_interp.x)
ypc = np.nanmin(bathy_interp.y)

# Area
xlenc = np.nanmax(bathy_interp.x) - np.nanmin(bathy_interp.x)
ylenc = np.nanmax(bathy_interp.y) - np.nanmin(bathy_interp.y)

# Grid cells
mxc = len(bathy_interp.x) - 1
myc = len(bathy_interp.y) - 1

# Resolution
res_x = bathy_interp.x[1].values - bathy_interp.x[0].values
res_y = bathy_interp.y[1].values - bathy_interp.y[0].values

# SWAN main grid parameters
grid = {
    'xpc': xpc,      # x origin
    'ypc': ypc,      # y origin
    'alpc': 0,       # x-axis direction 
    'xlenc': xlenc,  # grid length in x  
    'ylenc': ylenc,  # grid length in y
    'mxc': mxc,      # number mesh x, una menos pq si no SWAN peta
    'myc': myc,      # number mesh y, una menos pq si no SWAN peta
    'dxinp': res_x,  # size mesh x (resolution in x)
    'dyinp': res_y,  # size mesh y (resolution in y)
}

print('Bathymetry parameters:')
for k in grid.keys():
    print('{0:6}: {1}'.format(k, grid[k]))

Bathymetry parameters:
xpc   : 184.498816
ypc   : -21.57901999999979
alpc  : 0
xlenc : 0.7900000000037721
ylenc : 0.977999999999458
mxc   : 395
myc   : 489
dxinp : 0.0020000000000095497
dyinp : 0.0019999999999988916

---

Swan needs bathymetry to be flipped upside-down

depth = np.flipud(bathy_interp.z.values)

# store bathymetry data
pk.dump(grid, open(p_grid_param, "wb"))
pk.dump(depth, open(p_grid_depth, "wb"))