'''
This generates the NEMO Boundary. Creates indices for t, u and v points, plus rim gradient.
The variable names have been renamed to keep consistent with python standards and generalizing
the variable names eg. bdy_i is used instead of bdy_t
Ported from Matlab code by James Harle
@author: John Kazimierz Farey
@author: Srikanth Nagella
'''
# pylint: disable=E1103
# pylint: disable=no-name-in-module
#External Imports
import numpy as np
import logging
#Local Imports
from .utils.nemo_bdy_lib import sub2ind
class Boundary:
# Bearings for overlays
_NORTH = [1,-1,1,-1,2,None,1,-1]
_SOUTH = [1,-1,1,-1,None,-2,1,-1]
_EAST = [1,-1,1,-1,1,-1,2,None]
_WEST = [1,-1,1,-1,1,-1,None,-2]
def __init__(self, boundary_mask, settings, grid):
"""Generates the indices for NEMO Boundary and returns a Grid object with indices
Keyword arguments:
boundary_mask -- boundary mask
settings -- dictionary of setting values
grid -- type of the grid 't', 'u', 'v'
Attributes:
bdy_i -- index
bdy_r -- r index
"""
self.logger = logging.getLogger(__name__)
bdy_msk = boundary_mask
self.settings = settings
rw = self.settings['rimwidth']
rm = rw - 1
self.grid_type = grid.lower()
# Throw an error for wrong grid input type
if grid not in ['t', 'u', 'v', 'f']:
self.logger.error('Grid Type not correctly specified:'+grid)
raise ValueError("""%s is invalid grid grid_type;
must be 't', 'u', 'v' or 'f'""" %grid)
# Configure grid grid_type
if grid is not 't':
# We need to copy this before changing, because the original will be
# needed in calculating later grid boundary types
bdy_msk = boundary_mask.copy()
grid_ind = np.zeros(bdy_msk.shape, dtype=np.bool, order='C')
#NEMO works with a staggered 'C' grid. need to create a grid with staggered points
for fval in [-1, 0]: #-1 mask value, 0 Land, 1 Water/Ocean
if grid is 'u':
grid_ind[:, :-1] = np.logical_and(bdy_msk[:, :-1] == 1,
bdy_msk[:, 1:] == fval)
bdy_msk[grid_ind] = fval
elif grid is 'v':
grid_ind[:-1, :] = np.logical_and(bdy_msk[:-1, :] == 1,
bdy_msk[1:, :] == fval)
bdy_msk[grid_ind] = fval
elif grid is 'f':
grid_ind[:-1, :-1] = np.logical_and(bdy_msk[:-1,:-1] == 1,
bdy_msk[1:, 1:] == fval)
grid_ind[:-1, :] = np.logical_or(np.logical_and(bdy_msk[:-1, :] == 1,
bdy_msk[1:, :] == fval),
grid_ind[:-1, :] == 1)
grid_ind[:, :-1] = np.logical_or(np.logical_and(bdy_msk[:, :-1] == 1,
bdy_msk[:, 1:] == fval),
grid_ind[:, :-1] == 1)
bdy_msk[grid_ind] = fval
# Create padded array for overlays
msk = np.pad(bdy_msk,((1,1),(1,1)), 'constant', constant_values=(-1))
# create index arrays of I and J coords
igrid, jgrid = np.meshgrid(np.arange(bdy_msk.shape[1]), np.arange(bdy_msk.shape[0]))
SBi, SBj = self._find_bdy(igrid, jgrid, msk, self._SOUTH)
NBi, NBj = self._find_bdy(igrid, jgrid, msk, self._NORTH)
EBi, EBj = self._find_bdy(igrid, jgrid, msk, self._EAST)
WBi, WBj = self._find_bdy(igrid, jgrid, msk, self._WEST)
#create a 2D array index for the points that are on border
tij = np.column_stack((np.concatenate((SBi, NBi, WBi, EBi)),np.concatenate((SBj, NBj, WBj, EBj))))
bdy_i = np.tile(tij, (rw, 1, 1))
bdy_i = np.transpose(bdy_i, (1, 2, 0))
bdy_r = bdy_r = np.tile(np.arange(0,rw),(bdy_i.shape[0],1))
# Add points for relaxation zone over rim width
# In the relaxation zone with rim width. looking into the domain up to the rim width
# and select the points. S head North (0,+1) N head South(0,-1) W head East (+1,0)
# E head West (-1,0)
temp = np.column_stack((np.concatenate((SBi*0, NBi*0, WBi*0+1, EBi*0-1)),
np.concatenate((SBj*0+1, NBj*0-1, WBj*0, EBj*0))))
for i in range(rm):
bdy_i[:, :, i+1] = bdy_i[:, :, i] + temp
bdy_i = np.transpose(bdy_i, (1, 2, 0))
bdy_i = np.reshape(bdy_i,
(bdy_i.shape[0],bdy_i.shape[1]*bdy_i.shape[2]))
bdy_r = bdy_r.flatten(1)
## Remove duplicate and open sea points ##
bdy_i, bdy_r = self._remove_duplicate_points(bdy_i, bdy_r)
bdy_i, bdy_r, nonmask_index = self._remove_landpoints_open_ocean(bdy_msk, bdy_i, bdy_r)
### Fill in any gradients between relaxation zone and internal domain
### bdy_msk matches matlabs incarnation, r_msk is pythonic
r_msk = bdy_msk.copy()
r_msk[r_msk == 1] = rw
r_msk = np.float16(r_msk)
r_msk[r_msk < 1] = np.NaN
r_msk[bdy_i[:,1],bdy_i[:,0]] = np.float16(bdy_r)
r_msk_orig = r_msk.copy()
r_msk_ref = r_msk[1:-1, 1:-1]
self.logger.debug('Start r_msk bearings loop')
#Removes the shape gradients by smoothing it out
for i in range(rw-1):
# Check each bearing
for b in [self._SOUTH, self._NORTH, self._WEST, self._EAST]:
r_msk,r_msk_ref = self._fill(r_msk, r_msk_ref, b)
self.logger.debug('done loop')
# update bdy_i and bdy_r
new_ind = np.abs(r_msk - r_msk_orig) > 0
# The transposing gets around the Fortran v C ordering thing.
bdy_i_tmp = np.array([igrid.T[new_ind.T], jgrid.T[new_ind.T]])
bdy_r_tmp = r_msk.T[new_ind.T]
bdy_i = np.vstack((bdy_i_tmp.T, bdy_i))
uniqind = self._unique_rows(bdy_i)
bdy_i = bdy_i[uniqind, :]
bdy_r = np.hstack((bdy_r_tmp, bdy_r))
bdy_r = bdy_r[uniqind]
# sort by rimwidth
igrid = np.argsort(bdy_r, kind='mergesort')
bdy_r = bdy_r[igrid]
bdy_i = bdy_i[igrid, :]
self.bdy_i = bdy_i
self.bdy_r = bdy_r
self.logger.debug( 'Final bdy_i: %s', self.bdy_i.shape)
def _remove_duplicate_points(self, bdy_i, bdy_r):
""" Removes the duplicate points in the bdy_i and return the bdy_i and bdy_r
bdy_i -- bdy indexes
bdy_r -- bdy rim values
"""
bdy_i2 = np.transpose(bdy_i, (1, 0))
uniqind = self._unique_rows(bdy_i2)
bdy_i = bdy_i2[uniqind]
bdy_r = bdy_r[uniqind]
return bdy_i, bdy_r
def _remove_landpoints_open_ocean(self, mask, bdy_i, bdy_r):
""" Removes the land points and open ocean points """
unmask_index = mask[bdy_i[:,1],bdy_i[:,0]] != 0
bdy_i = bdy_i[unmask_index, :]
bdy_r = bdy_r[unmask_index]
return bdy_i, bdy_r, unmask_index
def _find_bdy(self, I, J, mask, brg):
"""Finds the border indexes by checking the change from ocean to land.
Returns the i and j index array where the shift happens.
Keyword arguments:
I -- I x direction indexes
J -- J y direction indexes
mask -- mask data
brg -- mask index range
"""
# subtract matrices to find boundaries, set to True
m1 = mask[brg[0]:brg[1], brg[2]:brg[3]]
m2 = mask[brg[4]:brg[5], brg[6]:brg[7]]
overlay = np.subtract(m1,m2)
# Create boolean array of bdy points in overlay
bool_arr = overlay==2
# index I or J to find bdies
bdy_I = I[bool_arr]
bdy_J = J[bool_arr]
return bdy_I, bdy_J
def _fill(self, mask, ref, brg):
""" """
tmp = mask[brg[4]:brg[5], brg[6]:brg[7]]
ind = (ref - tmp) > 1
ref[ind] = tmp[ind] + 1
mask[brg[0]:brg[1], brg[2]:brg[3]] = ref
return mask, ref
def _unique_rows(self, t):
""" This returns unique rows in the 2D array.
Returns indexes of unique rows in the input 2D array
t -- input 2D array
"""
sh = np.shape(t)
if (len(sh)> 2) or (sh[0] ==0) or (sh[1] == 0):
print('Warning: Shape of expected 2D array:', sh)
tlist = t.tolist()
sortt = []
indx = list(zip(*sorted([(val, i) for i,val in enumerate(tlist)])))[1]
indx = np.array(indx)
for i in indx:
sortt.append(tlist[i])
del tlist
for i,x in enumerate(sortt):
if x == sortt[i-1]:
indx[i] = -1
# all the rows are identical, set the first as the unique row
if sortt[0] == sortt[-1]:
indx[0] = 0
return indx[indx != -1]