Repository: British-Oceanographic-Data-Centre/COAsT (language: python)

Issue:
Climatology handling missing values
Added an argument to the create_climatology() routine to specify whether there are missing values in the dataset. If there are a different computation is used to create the climatology. Expanded the unit testing to check this is working.


Discussion / hints:

Fix patch:
diff --git a/coast/ALTIMETRY.py b/coast/ALTIMETRY.py
index b984a77a..f0e43335 100644
--- a/coast/ALTIMETRY.py
+++ b/coast/ALTIMETRY.py
@@ -193,8 +193,8 @@ def obs_operator(self, model, mod_var_name:str,
         self.dataset[new_var_name] = interpolated
     
     def crps(self, model_object, model_var_name, obs_var_name, 
-             nh_radius: float = 20, cdf_type:str='empirical', 
-             time_interp:str='linear', create_new_object = True):
+             nh_radius: float = 20, time_interp:str='linear', 
+             create_new_object = True):
         
         '''
         Comparison of observed variable to modelled using the Continuous
@@ -207,10 +207,7 @@ def crps(self, model_object, model_var_name, obs_var_name,
         model_object (model) : Model object (NEMO) containing model data
         model_var_name (str) : Name of model variable to compare.
         obs_var_name (str)   : Name of observed variable to compare.
-        nh_radius (float)    : Neighbourhood rad
-        cdf_type (str)       : Type of cumulative distribution to use for the
-                               model data ('empirical' or 'theoretical').
-                               Observations always use empirical.
+        nh_radius (float)    : Neighbourhood radius (km)
         time_interp (str)    : Type of time interpolation to use (s)
         create_new_obj (bool): If True, save output to new ALTIMETRY obj.
                                Otherwise, save to this obj.
@@ -234,7 +231,7 @@ def crps(self, model_object, model_var_name, obs_var_name,
                                obs_var.latitude.values, 
                                obs_var.values, 
                                obs_var.time.values, 
-                               nh_radius, cdf_type, time_interp )
+                               nh_radius, time_interp )
         if create_new_object:
             new_object = ALTIMETRY()
             new_dataset = self.dataset[['longitude','latitude','time']]
diff --git a/coast/CLIMATOLOGY.py b/coast/CLIMATOLOGY.py
new file mode 100644
index 00000000..3c1d7a28
--- /dev/null
+++ b/coast/CLIMATOLOGY.py
@@ -0,0 +1,79 @@
+from .COAsT import COAsT
+import numpy as np
+import xarray as xr
+import xarray.ufuncs as uf
+from dask.diagnostics import ProgressBar
+
+
+class CLIMATOLOGY(COAsT):
+    '''
+    A Python class containing methods for lazily creating climatologies of
+    NEMO data (or any xarray datasets) and writing to file. Also for resampling
+    methods.
+    '''
+    
+    def __init__(self):
+        return
+    
+    @staticmethod
+    def make_climatology(ds, output_frequency, monthly_weights = False, 
+                              time_var_name = 'time', time_dim_name = 't_dim',
+                              fn_out = None, missing_values = False):
+        '''
+        Calculates a climatology for all variables in a supplied dataset.
+        The resulting xarray dataset will NOT be loaded to RAM. Instead,
+        it is a set of dask operations. To load to RAM use, e.g. .compute().
+        However, if the original data was large, this may take a long time and
+        a lot of memory. Make sure you have the available RAM or chunking
+        and parallel processes are specified correctly.
+        
+        Otherwise, it is recommended that you access the climatology data
+        in an indexed way. I.E. compute only at specific parts of the data
+        are once.
+        
+        The resulting cliamtology dataset can be written to disk using
+        .to_netcdf(). Again, this may take a while for larger datasets.
+        
+        ds :: xarray dataset object from a COAsT object.
+        output_frequency :: any xarray groupby string. i.e:
+            'month'
+            'season'
+        time_var_name :: the string name of the time variable in dataset
+        time_dim_name :: the string name of the time dimension variable in dataset
+        fn_out :: string defining full output netcdf file path and name.
+        missing_values :: boolean where True indicates the data has missing values 
+            that should be ignored. Missing values must be represented by NaNs.
+        '''
+        
+        frequency_str = time_var_name + '.' + output_frequency
+        print('Calculating climatological mean')
+        
+        if missing_values:
+            ds_mean=xr.Dataset()
+            for varname, da in ds.data_vars.items():       
+                mask = xr.where(uf.isnan(da), 0, 1 )
+                data = da.groupby(frequency_str).sum(dim=time_dim_name) 
+                N = mask.groupby(frequency_str).sum(dim=time_dim_name)
+                ds_mean[varname] = data / N
+        else:
+            if monthly_weights:
+                month_length = ds[time_var_name].dt.days_in_month
+                grouped = month_length.groupby(frequency_str)
+            else:
+                ds['clim_mean_ones_tmp'] = (time_dim_name, np.ones(ds[time_var_name].shape[0]))
+                grouped = ds['clim_mean_ones_tmp'].groupby(frequency_str)
+
+            weights = grouped / grouped.sum()
+            ds_mean = (ds*weights).groupby(frequency_str).sum(dim=time_dim_name)
+
+            if not monthly_weights:
+                ds = ds.drop_vars('clim_mean_ones_tmp')
+
+        if fn_out is not None:
+            print('Saving to file. May take some time..')
+            with ProgressBar():
+                ds_mean.to_netcdf(fn_out)
+        
+        return ds_mean
+        
+        return
\ No newline at end of file
diff --git a/coast/CONTOUR.py b/coast/CONTOUR.py
index 9cb9daff..22a6466c 100755
--- a/coast/CONTOUR.py
+++ b/coast/CONTOUR.py
@@ -11,6 +11,7 @@
 from scipy.integrate import cumtrapz
 from sklearn.neighbors import BallTree
 from skimage import measure
+from .logging_util import get_slug, debug, info, warn, error
 
 # =============================================================================
 # The contour module is a place for code related to contours only
@@ -146,7 +147,7 @@ def __init__(self, nemo: COAsT, y_ind, x_ind, depth: int):
         '''
         Class defining a Contour type, which is a 3d dataset of points between a point A and 
         a point B defining an isobath contour. The dataset has a time, depth and contour dimension. 
-        The cotnour dimension defines the points along the contour.
+        The contour dimension defines the points along the contour.
         The supplied model Data is subsetted in its entirety along these dimensions and
         calculations can be performed on this dataset.
         
@@ -296,6 +297,10 @@ def calc_cross_contour_flow(self, nemo_u: COAsT, nemo_v: COAsT):
         (time, depth, position along contour) in m/s
         Contour_f.depth_integrated_normal_transport are the depth integrated 
         volume transports across the contour (time, position along contour) in Sv
+        
+        If the time dependent cell thicknesses (e3) on the u and v grids are 
+        present in the nemo_u and nemo_v datasets they will be used, if they 
+        are not then the initial cell thicknesses (e3_0) will be used.
 
         Parameters
         ----------
@@ -308,14 +313,40 @@ def calc_cross_contour_flow(self, nemo_u: COAsT, nemo_v: COAsT):
         -------
         None.
 
-        '''
-                
+        '''      
+        # compute transports flag; set to false if suitable e3 not found
+        compute_transports = True
+        
         # subset the u and v datasets 
         da_y_ind = xr.DataArray( self.y_ind, dims=['r_dim'] )
         da_x_ind = xr.DataArray( self.x_ind, dims=['r_dim'] )
         u_ds = nemo_u.dataset.isel(y_dim = da_y_ind, x_dim = da_x_ind)
         v_ds = nemo_v.dataset.isel(y_dim = da_y_ind, x_dim = da_x_ind)
-        
+                
+        # use time varying if e3 is present, if not default to e3_0
+        if 'e3' not in u_ds.data_vars:
+            if 'e3_0' not in u_ds.data_vars:
+                warn("e3 not found, transports will not be calculated")
+                compute_transports = False
+            else:
+                u_ds['e3'] = u_ds.e3_0.broadcast_like(u_ds.u_velocity) 
+        if 'e3' not in v_ds.data_vars:
+            if 'e3_0' not in v_ds.data_vars:
+                warn("e3 not found, transports will not be calculated")
+                compute_transports = False
+            else:
+                v_ds['e3'] = v_ds.e3_0.broadcast_like(v_ds.v_velocity)  
+           
+        # If time dimension is missing it can throw off the indexing so expand dims        
+        if 't_dim' not in u_ds.dims:
+            u_ds['u_velocity'] = u_ds.u_velocity.expand_dims('t_dim', axis=0)
+            if compute_transports:
+                u_ds['e3'] = u_ds.e3.expand_dims('t_dim', axis=0)
+        if 't_dim' not in v_ds.dims:
+            v_ds['v_velocity'] = v_ds.v_velocity.expand_dims('t_dim', axis=0)
+            if compute_transports:
+                v_ds['e3'] = v_ds.e3.expand_dims('t_dim', axis=0)
+                    
         dr_n = np.where(np.diff(self.y_ind)>0, np.arange(0,u_ds.r_dim.size-1), np.nan )
         dr_n = dr_n[~np.isnan(dr_n)].astype(int)
         dr_s = np.where(np.diff(self.y_ind)<0, np.arange(0,u_ds.r_dim.size-1), np.nan )
@@ -327,30 +358,31 @@ def calc_cross_contour_flow(self, nemo_u: COAsT, nemo_v: COAsT):
         
         # Note that subsetting the dataset first instead of subsetting each array seperately,
         # as we do here, is neater but significantly slower.
-        self.data_cross_flow['normal_velocities'] = xr.full_like(u_ds.vozocrtx, np.nan)        
-        self.data_cross_flow['normal_velocities'][:,:,dr_n] = u_ds.vozocrtx.data[:,:,dr_n+1]
-        self.data_cross_flow['normal_velocities'][:,:,dr_s] = -u_ds.vozocrtx.data[:,:,dr_s]
-        self.data_cross_flow['normal_velocities'][:,:,dr_e] = -v_ds.vomecrty.data[:,:,dr_e+1]
-        self.data_cross_flow['normal_velocities'][:,:,dr_w] = v_ds.vomecrty.data[:,:,dr_w]
+        self.data_cross_flow['normal_velocities'] = xr.full_like(u_ds.u_velocity, np.nan)        
+        self.data_cross_flow['normal_velocities'][:,:,dr_n] = u_ds.u_velocity.data[:,:,dr_n+1]
+        self.data_cross_flow['normal_velocities'][:,:,dr_s] = -u_ds.u_velocity.data[:,:,dr_s]
+        self.data_cross_flow['normal_velocities'][:,:,dr_e] = -v_ds.v_velocity.data[:,:,dr_e+1]
+        self.data_cross_flow['normal_velocities'][:,:,dr_w] = v_ds.v_velocity.data[:,:,dr_w]
         self.data_cross_flow['normal_velocities'].attrs = {'units':'m/s', \
                 'standard_name':'contour-normal velocities'}
-             
-        self.data_cross_flow['normal_transport'] = xr.full_like(u_ds.vozocrtx, np.nan)  
-        self.data_cross_flow['normal_transport'][:,:,dr_n] = ( u_ds.vozocrtx.data[:,:,dr_n+1] * 
-                                u_ds.e2.data[dr_n+1] * u_ds.e3_0.data[:,dr_n+1] )
-        self.data_cross_flow['normal_transport'][:,:,dr_s] = ( -u_ds.vozocrtx.data[:,:,dr_s] * 
-                                u_ds.e2.data[dr_s] * u_ds.e3_0.data[:,dr_s] )
-        self.data_cross_flow['normal_transport'][:,:,dr_e] = ( -v_ds.vomecrty.data[:,:,dr_e+1] *
-                                v_ds.e1.data[dr_e+1] * v_ds.e3_0.data[:,dr_e+1] )
-        self.data_cross_flow['normal_transport'][:,:,dr_w] = ( v_ds.vomecrty.data[:,:,dr_w] *
-                                v_ds.e1.data[dr_w] * v_ds.e3_0.data[:,dr_w] )
-        self.data_cross_flow['normal_transport'].attrs = {'units':'m^3/s', \
-                'standard_name':'contour-normal volume transport'}
-        
-        self.data_cross_flow['depth_integrated_normal_transport'] = (self.data_cross_flow
-                                .normal_transport.sum(dim='z_dim') / 1000000.)
-        self.data_cross_flow['depth_integrated_normal_transport'].attrs ={'units':'Sv', \
-                'standard_name':'contour-normal depth integrated volume transport'}
+        
+        if compute_transports:    
+            self.data_cross_flow['normal_transport'] = xr.full_like(u_ds.u_velocity, np.nan)  
+            self.data_cross_flow['normal_transport'][:,:,dr_n] = ( u_ds.u_velocity.data[:,:,dr_n+1] * 
+                                    u_ds.e2.data[dr_n+1] * u_ds.e3.data[:,:,dr_n+1] )
+            self.data_cross_flow['normal_transport'][:,:,dr_s] = ( -u_ds.u_velocity.data[:,:,dr_s] * 
+                                    u_ds.e2.data[dr_s] * u_ds.e3.data[:,:,dr_s] )
+            self.data_cross_flow['normal_transport'][:,:,dr_e] = ( -v_ds.v_velocity.data[:,:,dr_e+1] *
+                                    v_ds.e1.data[dr_e+1] * v_ds.e3.data[:,:,dr_e+1] )
+            self.data_cross_flow['normal_transport'][:,:,dr_w] = ( v_ds.v_velocity.data[:,:,dr_w] *
+                                    v_ds.e1.data[dr_w] * v_ds.e3.data[:,:,dr_w] )
+            self.data_cross_flow['normal_transport'].attrs = {'units':'m^3/s', \
+                    'standard_name':'contour-normal volume transport'}
+            
+            self.data_cross_flow['depth_integrated_normal_transport'] = (self.data_cross_flow
+                                    .normal_transport.sum(dim='z_dim') / 1000000.)
+            self.data_cross_flow['depth_integrated_normal_transport'].attrs ={'units':'Sv', \
+                    'standard_name':'contour-normal depth integrated volume transport'}
                                 
         self.__update_cross_flow_vars('depth_0',u_ds.depth_0,v_ds.depth_0,dr_n,dr_s,dr_e,dr_w,1)
         self.__update_cross_flow_vars('longitude',u_ds.longitude,v_ds.longitude,dr_n,dr_s,dr_e,dr_w,0)
@@ -359,8 +391,9 @@ def calc_cross_contour_flow(self, nemo_u: COAsT, nemo_v: COAsT):
         self.__update_cross_flow_vars('e1',u_ds.e1,v_ds.e1,dr_n,dr_s,dr_e,dr_w,0)
         self.data_cross_flow['e2'] = xr.full_like(self.data_contour.e2, np.nan)
         self.__update_cross_flow_vars('e2',u_ds.e2,v_ds.e2,dr_n,dr_s,dr_e,dr_w,0)
-        self.data_cross_flow['e3_0'] = xr.full_like(self.data_contour.e3_0, np.nan)
-        self.__update_cross_flow_vars('e3_0',u_ds.e3_0,v_ds.e3_0,dr_n,dr_s,dr_e,dr_w,1)
+        self.data_cross_flow['e3'] = xr.full_like(self.data_cross_flow.normal_velocities, np.nan)
+        if compute_transports:
+            self.__update_cross_flow_vars('e3',u_ds.e3,v_ds.e3,dr_n,dr_s,dr_e,dr_w,2)
         
         self.data_cross_flow['depth_0'].attrs = {'standard_name':'Depth at time zero \
                 on the contour-normal velocity grid points'}
@@ -368,6 +401,7 @@ def calc_cross_contour_flow(self, nemo_u: COAsT, nemo_v: COAsT):
                 the contour-normal velocity grid points'}
         self.data_cross_flow['longitude'].attrs = {'standard_name':'Longitude at \
                 the contour-normal velocity grid points'}
+        self.data_cross_flow = self.data_cross_flow.squeeze()
 
     
     def __update_cross_flow_vars(self, var, u_var,v_var, dr_n, dr_s, dr_e, dr_w, pos ):
diff --git a/coast/DISTRIBUTION.py b/coast/DISTRIBUTION.py
new file mode 100644
index 00000000..0edc513a
--- /dev/null
+++ b/coast/DISTRIBUTION.py
@@ -0,0 +1,192 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from . import stats_util
+from .logging_util import get_slug, debug, error
+import scipy.stats
+import xarray as xr
+
+class DISTRIBUTION:
+    '''
+    An object for storing Cumulative Distribution Function information.
+    Used primarily for calculating the Continuous Ranked Probability Score.
+    The object is initalisated by passing in a 'sample' vector of data for
+    which to determine a CDF. This sample is used to construct either a 
+    empirical or theoretical CDF, depending on the cdf_type argument.
+    '''
+    
+    def __init__(self, sample):
+        """Initialisation of CDF object.
+        Args:
+            sample (array): Data sample over which to estimate CDF
+            cdf_type (str): Either 'empirical' or 'theoretical'.
+            cdf_func (str): Function type if cdf_type='theoretical'. Presently
+                            only 'gaussian'.
+        Returns:
+            New CDF object.
+        """
+        debug(f"Creating a new {get_slug(self)}")
+        self.sample = sample
+        self.sample_size = len(sample)
+        self.mu       = np.nanmean(sample)
+        self.sigma    = np.nanstd(sample)
+        self.plot_xmin, self.plot_xmax = self.set_x_bounds()
+        debug(f"{get_slug(self)} initialised")
+        
+    def set_x_bounds(self):
+        ''' Calculate x bounds for CDF plotting '''
+        # Is input a single value (st. dev == 0)
+        single_value = True if self.sigma == 0 else False  # TODO This could just be: single_value = self.sigma == 0
+        if single_value: 
+            self.cdf_type = 'empirical'
+        
+        # Calculate x bounds as 5 std dev. either side of the mean
+        debug(f"Calculating x bounds for {get_slug(self)}")
+        if single_value:
+            xmin = self.mu-1
+            xmax = self.mu+1
+        else:
+            xmin = self.mu - 5*self.sigma
+            xmax = self.mu + 5*self.sigma
+        return xmin, xmax
+    
+    def build_discrete_cdf(self, x: np.ndarray = None, n_pts: int = 1000):
+        """Builds a discrete CDF for plotting and direct comparison.
+        Args:
+            x (array): x-values over which to calculate discrete CDF values.
+                       If none, these will be determined using mu and sigma.
+            n_pts (int): n_pts to use for x if x=None.
+        Returns:
+            x and y arrays for discrete CDF.
+        """
+        debug(f"Discrete CDF will be built for {get_slug(self)}")
+        # Build discrete X bounds
+        if x is None:
+            debug(f"Building discrete X bounds for {get_slug(self)}")
+            x = np.linspace(self.plot_xmin, self.plot_xmax, n_pts)
+        
+        if self.cdf_type == "empirical":
+            y = stats_util.empirical_distribution(x, self.sample)
+
+        else:
+            error(f"CDF type for {get_slug(self)} is , which is not acceptable, raising exception!")
+
+        return x, y
+        
+    @staticmethod
+    def normal_distribution(mu: float=0, sigma: float=1, 
+                            x: np.ndarray=None, n_pts: int=1000):
+        """Generates a discrete normal distribution.
+    
+        Keyword arguments:
+        x     -- Arbitrary array of x-values
+        mu    -- Distribution mean
+        sigma -- Distribution standard deviation
+        
+        return: Array of len(x) containing the normal values calculated from
+                the elements of x.
+        """
+        if x is None:
+            x = np.linspace( mu-5*sigma, mu+5*sigma, n_pts)
+        debug(f"Generating normal distribution for {get_slug(x)}")
+        term1 = sigma*np.sqrt( 2*np.pi )
+        term1 = 1/term1
+        exponent = -0.5*((x-mu)/sigma)**2
+        return term1*np.exp( exponent )
+    
+    @staticmethod
+    def cumulative_distribution(mu: float=0, sigma: float=1, 
+                                x: np.ndarray=None, cdf_func: str='gaussian'):
+        """Integrates under a discrete PDF to obtain an estimated CDF.
+    
+        Keyword arguments:
+        x   -- Arbitrary array of x-values
+        pdf -- PDF corresponding to values in x. E.g. as generated using
+               normal_distribution.
+        
+        return: Array of len(x) containing the discrete cumulative values 
+                estimated using the integral under the provided PDF.
+        """
+        debug(f"Estimating CDF using {get_slug(x)}")
+        if cdf_func=='gaussian': #If Gaussian, integrate under pdf
+            pdf = DISTRIBUTION.normal_distribution(mu=mu, sigma=sigma, x=x)
+            cdf = [np.trapz(pdf[:ii],x[:ii]) for ii in range(0,len(x))]
+        else: 
+            raise NotImplementedError
+        return np.array(cdf)
+    
+    @staticmethod
+    def empirical_distribution(x, sample):
+        """Estimates a CDF empirically.
+    
+        Keyword arguments:
+        x      -- Array of x-values over which to generate distribution
+        sample -- Sample to use to generate distribution
+        
+        return: Array of len(x) containing corresponding EDF values
+        """
+        debug(f"Estimating empirical distribution with {get_slug(x)}")
+        sample = np.array(sample)
+        sample = sample[~np.isnan(sample)]
+        sample = np.sort(sample)
+        edf = np.zeros(len(x))
+        n_sample = len(sample)
+        for ss in sample:
+            edf[x>ss] = edf[x>ss] + 1/n_sample
+        return xr.DataArray(edf)
+        
+    def get_common_x(self, other, n_pts=2000):
+        """Generates a common x vector for two CDF objects."""
+        debug(f"Generating common X vector for {get_slug(self)} and {get_slug(other)}")
+        xmin = min(self.plot_xmin, other.plot_xmin)
+        xmax = max(self.plot_xmax, other.plot_xmax)
+        common_x = np.linspace(xmin, xmax, n_pts)
+        return common_x
+    
+    def plot_cdf(self):
+        """ A quick plot showing the CDF contained in this object."""
+        debug(f"Generating quick plot for {get_slug(self)}")
+        ax = plt.subplot(111)
+        x,y = self.build_discrete_cdf()
+        ax.plot(x, y)
+        ax.grid()
+        return
+    
+    def integrate_cdf(self, other=None, plot=False):
+        """Returns the integral under CDF or between two CDFs. This is 
+        equivalent to the first order Wasserstein metric for two 
+        probability distributions"""
+        debug(f"Generating diff plot for {get_slug(self)} and {get_slug(other)}")
+
+        if other is None:
+            integral = scipy.stats.wasserstein_distance(self.sample, [0])
+            if plot:
+                x = self.get_common_x(other)
+                x, y1 = self.build_discrete_cdf(x)
+                dum, y2 = other.build_discrete_cdf(x)
+                fig = plt.figure()
+                ax = plt.subplot(111)
+                ax.plot(x, y1, c='k', linestyle='--')
+                ax.plot(x, y2, linestyle='--')
+                ax.fill_between(x, y1, y2, alpha=0.5)
+                ax.set_title('Area: ' + str(integral))
+                plt.legend(('1','2'))
+            
+        else:
+            integral = scipy.stats.wasserstein_distance(self.sample, other.sample)
+            
+            if plot:
+                x = self.get_common_x(other)
+                x, y1 = self.build_discrete_cdf(x)
+                dum, y2 = other.build_discrete_cdf(x)
+                fig = plt.figure()
+                ax = plt.subplot(111)
+                ax.plot(x, y1, c='k', linestyle='--')
+                ax.plot(x, y2, linestyle='--')
+                ax.fill_between(x, y1, y2, alpha=0.5)
+                ax.set_title('Area: ' + str(integral))
+                plt.legend(('1','2'))
+                
+        if plot:
+            return integral, fig, ax
+        else:
+            return integral
\ No newline at end of file
diff --git a/coast/MASK_MAKER.py b/coast/MASK_MAKER.py
new file mode 100644
index 00000000..01eedff0
--- /dev/null
+++ b/coast/MASK_MAKER.py
@@ -0,0 +1,142 @@
+from .COAsT import COAsT
+import xarray as xr
+import numpy as np
+import skimage.draw as draw
+from . import general_utils
+
+class MASK_MAKER(): 
+
+    def __init__(self):
+        
+        return
+    
+    @staticmethod
+    def fill_polygon_by_index(array_to_fill, vertices_r, vertices_c, 
+                              fill_value = 1, additive = False):
+        """
+        Draws and fills a polygon onto an existing numpy array based on array
+        indices. To create a new mask, give np.zeros(shape) as input. 
+        Polygon vertices are drawn in the order given.
+        
+        Parameters
+        ----------
+        array_to_fill (2D array): Array onto which to fill polygon
+        vertices_r (1D array): Row indices for polygon vertices
+        vertices_c (1D_array): Column indices for polygon vertices
+        fill_value (float, bool or int): Fill value for polygon (Default: 1)
+        additive (bool): If true, add fill value to existing array. Otherwise
+                         indices will be overwritten. (Default: False)
+
+        Returns
+        -------
+        Filled 2D array
+        """
+        array_to_fill = np.array(array_to_fill)
+        polygon_ind = draw.polygon(vertices_r, vertices_c, 
+                                       array_to_fill.shape)
+        if additive:
+            array_to_fill[polygon_ind[0], polygon_ind[1]] += fill_value
+        else:
+            array_to_fill[polygon_ind[0], polygon_ind[1]] = fill_value
+        return array_to_fill
+    
+    @staticmethod
+    def fill_polygon_by_lonlat(array_to_fill, longitude, latitude, 
+                               vertices_lon, vertices_lat, fill_value = 1,
+                               additive = False):
+        """
+        Draws and fills a polygon onto an existing numpy array based on 
+        vertices defined by longitude and latitude locations. This does NOT
+        draw a polygon on a sphere, but instead based on straight lines 
+        between points. This is OK for small regional areas, but not advisable
+        for large and global regions.
+        Polygon vertices are drawn in the order given.
+        
+        Parameters
+        ----------
+        array_to_fill (2D array): Array onto which to fill polygon
+        vertices_r (1D array): Row indices for polygon vertices
+        vertices_c (1D_array): Column indices for polygon vertices
+        fill_value (float, bool or int): Fill value for polygon (Default: 1)
+        additive (bool): If true, add fill value to existing array. Otherwise
+                         indices will be overwritten. (Default: False)
+
+        Returns
+        -------
+        Filled 2D array
+        """
+        array_to_fill = np.array(array_to_fill)
+        ind2D = general_utils.nearest_indices_2D(longitude, latitude, 
+                                                 vertices_lon, vertices_lat)
+        
+        polygon_ind = draw.polygon(ind2D[1], ind2D[0], array_to_fill.shape)
+        if additive:
+            array_to_fill[polygon_ind[0], polygon_ind[1]] += fill_value
+        else:
+            array_to_fill[polygon_ind[0], polygon_ind[1]] = fill_value
+        return array_to_fill
+    
+    @classmethod
+    def region_def_nws_north_sea(cls, longitude, latitude, bath):
+        '''
+        Regional definition for the North Sea (Northwest European Shelf)
+        Longitude, latitude and bath should be 2D arrays corresponding to model
+        coordinates and bathymetry. Bath should be positive with depth.
+        '''
+        vertices_lon = [-5.34,  -0.7,   7.5,    7.5,   9,     9, 
+                        6.3,    6.3,    5,      5,     4.126, 4.126, -1.071]
+        vertices_lat = [56.93,  54.09,  54.09,  56,    56,    57.859, 
+                        57.859, 58.121, 58.121, 58.59, 58.59, 60.5,  60.5]
+
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape),
+                                      longitude, latitude,
+                                      vertices_lon, vertices_lat)
+        mask = mask * (bath<200) * (bath>0) * (~xr.ufuncs.isnan(bath))
+        return mask
+
+    @classmethod
+    def region_def_nws_outer_shelf(cls, longitude, latitude, bath):
+        '''
+        Regional definition for the Outher Shelf (Northwest European Shelf)
+        Longitude, latitude and bath should be 2D arrays corresponding to model
+        coordinates and bathymetry. Bath should be positive with depth.
+        '''
+        vertices_lon = [-4,    -9.5,  -1,   3.171, 3.171, -3.76, -3.76, -12, 
+                        -12,   -12,   -4]
+        vertices_lat = [50.5,  52.71, 60.5, 60.45, 63.3,  63.3,  60.45, 60.45, 
+                        55.28, 48,    48]
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape),
+                                      longitude, latitude,
+                                      vertices_lon, vertices_lat)
+        mask = mask * (bath<200) * (bath>0) * (~xr.ufuncs.isnan(bath))
+        return mask
+
+    @classmethod
+    def region_def_nws_norwegian_trench(cls, longitude, latitude, bath):
+        '''
+        Regional definition for the Norwegian Trench (Northwest European Shelf)
+        Longitude, latitude and bath should be 2D arrays corresponding to model
+        coordinates and bathymetry. Bath should be positive with depth.
+        '''
+        vertices_lon = [10.65, 1.12,  1.12, 10.65]
+        vertices_lat = [61.83, 61.83, 48,   48]
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape),
+                                      longitude, latitude,
+                                      vertices_lon, vertices_lat)
+        mask = mask * (bath>200) * (bath>0) * (~xr.ufuncs.isnan(bath))
+        return mask
+
+    @classmethod    
+    def region_def_nws_english_channel(cls, longitude, latitude, bath):
+        '''
+        Regional definition for the English Channel (Northwest European Shelf)
+        Longitude, latitude and bath should be 2D arrays corresponding to model
+        coordinates and bathymetry. Bath should be positive with depth.
+        '''
+        vertices_lon = [7.57, 7.57,  -0.67,  -2,   -3.99, -3.99, -3.5, 12, 14]
+        vertices_lat = [56,   54.08, 54.08,  50.7, 50.7,  48.8,  48,   48, 56]
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape),
+                                      longitude, latitude,
+                                      vertices_lon, vertices_lat)
+        mask = mask * (bath<200) * (bath>0) * (~xr.ufuncs.isnan(bath))
+        return mask
\ No newline at end of file
diff --git a/coast/NEMO.py b/coast/NEMO.py
index 4c7b8a4f..6803b82e 100644
--- a/coast/NEMO.py
+++ b/coast/NEMO.py
@@ -56,13 +56,13 @@ def set_grid_vars(self):
         """ Define the variables to map from the domain file to the NEMO obj"""
         # Define grid specific variables to pull across
         if self.grid_ref == 'u-grid':
-            self.grid_vars = ['glamu', 'gphiu', 'e1u', 'e2u', 'e3u_0', 'depthu_0'] #What about e3vw
+            self.grid_vars = ['glamu', 'gphiu', 'e1u', 'e2u', 'e3u_0', 'depthu_0',] #What about e3vw
         elif self.grid_ref == 'v-grid':
             self.grid_vars = ['glamv', 'gphiv', 'e1v', 'e2v', 'e3v_0', 'depthv_0']
         elif self.grid_ref == 't-grid':
-            self.grid_vars = ['glamt', 'gphit', 'e1t', 'e2t', 'e3t_0', 'deptht_0', 'tmask']
+            self.grid_vars = ['glamt', 'gphit', 'e1t', 'e2t', 'e3t_0', 'deptht_0', 'tmask', 'bottom_level']
         elif self.grid_ref == 'w-grid':
-            self.grid_vars = ['glamt', 'gphit', 'e1t', 'e2t', 'e3w_0', 'depthw_0']
+            self.grid_vars = ['glamt', 'gphit', 'e1t', 'e2t', 'e3w_0', 'depthw_0', 'bottom_level']
         elif self.grid_ref == 'f-grid':
             self.grid_vars = ['glamf', 'gphif', 'e1f', 'e2f', 'e3f_0', 'depthf_0']
 
@@ -86,7 +86,11 @@ def set_variable_mapping(self):  # TODO Add a docstring
                             'toce' : 'temperature',
                             'so' : 'salinity',
                             'vosaline' : 'salinity',
-                            'voce' : 'salinity',
+                            'soce' : 'salinity',
+                            'voce' : 'v_velocity',
+                            'vomecrty' : 'v_velocity',
+                            'uoce' : 'u_velocity',
+                            'vozocrtx' : 'u_velocity',
                             'sossheig' : 'ssh',
                             'zos' : 'ssh' }
         # Variable names mapped from domain to NEMO object
@@ -108,7 +112,7 @@ def set_variable_mapping(self):  # TODO Add a docstring
                                    'depthf_0':'depth_0',
                                    'depthu_0':'depth_0', 'depthv_0':'depth_0',
                                    'depthw_0':'depth_0', 'deptht_0':'depth_0',
-                                   'ln_sco':'ln_sco'}
+                                   'ln_sco':'ln_sco', 'bottom_level':'bottom_level'}
 
     # TODO Add parameter type hints and a docstring
     def load_domain(self, fn_domain, chunks):  # TODO Do something with this unused parameter or remove it
@@ -635,6 +639,151 @@ def apply_doodson_x0_filter(self, var_str):
         if filtered is not None:
             self.dataset[new_var_str] = (old_dims, filtered)
         return
+       
+    
+    @staticmethod
+    def get_e3_from_ssh(nemo_t, e3t=True, e3u=False, e3v=False, e3f=False, 
+                        e3w=False, dom_fn:str=None):
+        '''
+        Where the model has been run with a nonlinear free surface
+        and z* variable volumne (ln_vvl_zstar=True) then the vertical scale factors
+        will vary in time (and space). This function will compute the vertical
+        scale factors e3t, e3u, e3v, e3f and e3w by using the sea surface height
+        field (ssh variable) and initial scale factors from the domain_cfg file.
+        The vertical scale factors will be computed at the same model time as the 
+        ssh and if the ssh field is averaged in time then the scale factors will
+        also be time averages. 
+        
+        A t-grid NEMO object containing the ssh variable must be passed in. Either 
+        the domain_cfg path must have been passed in as an argument when the NEMO 
+        object was created or it must be passed in here using the dom_fn argument.
+        
+        e.g. e3t,e3v,e3f = coast.NEMO.get_e3_from_ssh(nemo_t,true,false,true,true,false)
+        
+        Parameters
+        ----------
+        nemo_t : (COAsT.NEMO), NEMO object on the t-grid containing the ssh variable
+        e3t : (boolean), true if e3t is to be returned. Default True.
+        e3u : (boolean), true if e3u is to be returned. Default False.
+        e3v : (boolean), true if e3v is to be returned. Default False.
+        e3f : (boolean), true if e3f is to be returned. Default False.
+        e3w : (boolean), true if e3w is to be returned. Default False.
+        dom_fn : (str), Optional, path to domain_cfg file. 
+
+        Returns
+        -------
+        Tuple of xarray.DataArrays
+        (e3t, e3u, e3v, e3f, e3w)
+        Only those requested will be returned, but the ordering is always the same.
+
+        '''
+        e3_return = []
+        try:
+            ssh = nemo_t.dataset.ssh
+        except AttributeError:
+            print('The nemo_t dataset must contain the ssh variable.')
+            return
+        if 't_dim' not in ssh.dims:
+            ssh = ssh.expand_dims('t_dim', axis=0)
+
+        # Load domain_cfg
+        if dom_fn is None:
+            dom_fn = nemo_t.filename_domain
+        try:
+            ds_dom = xr.open_dataset(dom_fn).squeeze().rename(
+                {'z':'z_dim', 'x':'x_dim', 'y':'y_dim'})
+        except OSError:
+            print(f'Problem opening domain_cfg file: {dom_fn}')
+            return 
+            
+        e3t_0 = ds_dom.e3t_0
+    
+        # Water column thickness, i.e. depth of bottom w-level on horizontal t-grid
+        H = e3t_0.cumsum(dim='z_dim').isel(z_dim=ds_dom.bottom_level.astype("int")-1)
+        # Add correction to e3t_0 due to change in ssh
+        e3t_new = e3t_0 * ( 1 + ssh / H )
+        # preserve dimension ordering
+        e3t_new = e3t_new.transpose('t_dim','z_dim','y_dim','x_dim')
+        # mask out correction at layers below bottom level
+        e3t_new = e3t_new.where(e3t_new.z_dim<ds_dom.bottom_level,e3t_0.data)
+        # preserve any other t mask
+        e3t_new = e3t_new.where(~np.isnan(ssh))
+        if e3t:
+            e3_return.append(e3t_new.squeeze())
+        
+        if np.any([e3u,e3v,e3f]):
+            e1e2t = ds_dom.e1t * ds_dom.e2t
+        if np.any([e3u,e3v,e3w]):
+            e3t_dt = e3t_new - e3t_0
+            
+        # area averaged interpolation onto the u-grid to get e3u    
+        if np.any([e3u,e3f]):           
+            e1e2u = ds_dom.e1u * ds_dom.e2u  
+            # interpolate onto u-grid
+            e3u_temp = ( ( (0.5/e1e2u[:,:-1]) * ( (e1e2t[:,:-1] * e3t_dt[:,:,:,:-1]) 
+                    + (e1e2t[:,1:] * e3t_dt[:,:,:,1:]) ) )
+                    .transpose('t_dim','z_dim','y_dim','x_dim') )
+            # u mask
+            e3u_temp = e3u_temp.where(e3t_dt[:,:,:,1:]!=0, 0)
+            # mask out correction at layers below bottom level
+            e3u_temp = e3u_temp.where(e3u_temp.z_dim<ds_dom.bottom_level[:,:-1], 0)
+            # Add correction to e3u_0
+            e3u_temp = ( e3u_temp + ds_dom.e3u_0[:,:,:-1] )
+            e3u_new = xr.zeros_like(e3t_new)
+            e3u_new[:,:,:,:-1] = e3u_temp
+            e3u_new[:,:,:,-1] = ds_dom.e3u_0[:,:,-1]
+            e3u_new['longitude'] = ds_dom.glamu
+            e3u_new['latitude'] = ds_dom.gphiu
+            if e3u:
+                e3_return.append(e3u_new.squeeze())
+        
+        # area averaged interpolation onto the u-grid to get e3v 
+        if e3v:
+            e1e2v = ds_dom.e1v * ds_dom.e2v
+            e3v_temp = ( ( (0.5/e1e2v[:-1,:]) * ( (e1e2t[:-1,:] * e3t_dt[:,:,:-1,:]) 
+                    + (e1e2t[1:,:] * e3t_dt[:,:,1:,:]) ) )
+                    .transpose('t_dim','z_dim','y_dim','x_dim') )
+            e3v_temp = e3v_temp.where(e3t_dt[:,:,1:,:]!=0, 0)
+            e3v_temp = e3v_temp.where(e3v_temp.z_dim<ds_dom.bottom_level[:-1,:], 0)
+            e3v_temp = ( e3v_temp + ds_dom.e3v_0[:,:-1,:] )
+            e3v_new = xr.zeros_like(e3t_new)
+            e3v_new[:,:,:-1,:] = e3v_temp
+            e3v_new[:,:,-1,:] = ds_dom.e3v_0[:,-1,:]
+            e3v_new['longitude'] = ds_dom.glamv
+            e3v_new['latitude'] = ds_dom.gphiv
+            e3_return.append(e3v_new.squeeze())
+        
+        # area averaged interpolation onto the u-grid to get e3f
+        if e3f:
+            e1e2f = ds_dom.e1f * ds_dom.e2f
+            e3u_dt = e3u_new - ds_dom.e3u_0
+            e3f_temp = ( ( (0.5/e1e2f[:-1,:]) * ( (e1e2u[:-1,:] * e3u_dt[:,:,:-1,:]) 
+                    + (e1e2u[1:,:] * e3u_dt[:,:,1:,:] ) ) )
+                    .transpose('t_dim','z_dim','y_dim','x_dim') )
+            e3f_temp = e3f_temp.where(e3u_dt[:,:,1:,:]!=0, 0)
+            e3f_temp = e3f_temp.where(e3f_temp.z_dim<ds_dom.bottom_level[:-1,:], 0)        
+            e3f_temp = ( e3f_temp + ds_dom.e3f_0[:,:-1,:] )
+            e3f_new = xr.zeros_like(e3t_new)
+            e3f_new[:,:,:-1,:] = e3f_temp
+            e3f_new[:,:,-1,:] = ds_dom.e3f_0[:,-1,:]
+            e3f_new['longitude'] = ds_dom.glamf
+            e3f_new['latitude'] = ds_dom.gphif
+            e3_return.append(e3f_new.squeeze())
+           
+        # simple vertical interpolation for e3w. Special treatment of top and bottom levels   
+        if e3w:
+            # top levels correction same at e3t
+            e3w_new = (ds_dom.e3w_0 + e3t_dt).transpose('t_dim','z_dim','y_dim','x_dim')
+            # levels between top and bottom
+            e3w_new[dict(z_dim=slice(1,None))] = ( 0.5*e3t_dt[:,:-1,:,:] 
+                                    + 0.5*e3t_dt[:,1:,:,:] 
+                                    + ds_dom.e3w_0[1:,:,:] )
+            # bottom and below levels
+            e3w_new = e3w_new.where(e3w_new.z_dim<ds_dom.bottom_level, e3t_dt.shift(z_dim=1) + ds_dom.e3w_0)
+            e3_return.append(e3w_new.squeeze())
+            
+        return tuple(e3_return)
+            
     
     def harmonics_combine(self, constituents, components = ['x','y']):
         '''
diff --git a/coast/PROFILE.py b/coast/PROFILE.py
index ee0ca448..3fddfca2 100755
--- a/coast/PROFILE.py
+++ b/coast/PROFILE.py
@@ -2,6 +2,8 @@
 import xarray as xr
 from . import general_utils, plot_util, crps_util, COAsT
 import matplotlib.pyplot as plt
+import glob
+import datetime
 
 class PROFILE(COAsT):
     '''
@@ -24,14 +26,56 @@ def __init__(self):
         return
     
     def read_EN4(self,fn_en4, multiple = False, chunks = {}):
+        '''
+        Reads a single or multiple EN4 netCDF files into the COAsT profile
+        data structure.
+        '''
         if not multiple:
             self.dataset = xr.open_dataset(fn_en4, chunks = chunks)
         else:
-            self.dataset = xr.open_mfdataset(fn_en4, chunks = chunks)
-            
+            if type(fn_en4) is not list:
+                fn_en4 = [fn_en4]
+                
+            file_to_read = []
+            for file in fn_en4:
+                if '*' in file:
+                    wildcard_list = glob.glob(file)
+                    file_to_read = file_to_read + wildcard_list
+                else:
+                    file_to_read.append(file)
+                    
+            # Reorder files to read 
+            file_to_read = np.array(file_to_read)
+            dates = [ff[-9:-3] for ff in file_to_read]
+            dates = [datetime.datetime(int(dd[0:4]), int(dd[4:6]),1) for dd in dates]
+            sort_ind = np.argsort(dates)
+            file_to_read = file_to_read[sort_ind]
+                    
+            for ff in range(0,len(file_to_read)):
+                file = file_to_read[ff]
+                data_tmp = xr.open_dataset(file, chunks = chunks)
+                if ff==0:
+                    self.dataset = data_tmp
+                else:
+                    self.dataset=xr.concat((self.dataset, data_tmp),dim='N_PROF')
+                
+        
         rename_vars = {'LATITUDE':'latitude', 'LONGITUDE' : 'longitude',
-                       'DEPH_CORRECTED' : 'depth'}
-        self.dataset = self.dataset.rename(rename_vars)
+                       'DEPH_CORRECTED' : 'depth', 'JULD':'time',
+                       'POTM_CORRECTED':'potential_temperature',
+                       'TEMP':'temperature', 'PSAL_CORRECTED':'practical_salinity',
+                       'POTM_CORRECTED_QC':'qc_potential_temperature',
+                       'PSAL_CORRECTED_QC':'qc_practical_salinity',
+                       'DEPH_CORRECTED_QC':'qc_depth',
+                       'JULD_QC':'qc_time'}
+        rename_dims = {'N_PROF':'profile', 'N_PARAM':'parameter', 
+                       'N_LEVELS':'z_dim',}
+        vars_to_keep = list(rename_vars.keys())
+        coords = ['LATITUDE','LONGITUDE','JULD']
+        self.dataset = self.dataset.set_coords(coords)
+        self.dataset = self.dataset.rename(rename_dims)
+        self.dataset = self.dataset[vars_to_keep].rename_vars(rename_vars)
+        return
         
 ##############################################################################
 ###                ~            Manipulate           ~                     ###
@@ -39,7 +83,6 @@ def read_EN4(self,fn_en4, multiple = False, chunks = {}):
     
     def subset_indices_lonlat_box(self, lonbounds, latbounds):
         """Generates array indices for data which lies in a given lon/lat box.
-
         Keyword arguments:
         lon       -- Longitudes, 1D or 2D.
         lat       -- Latitudes, 1D or 2D
@@ -63,12 +106,12 @@ def plot_profile(self, var:str, profile_indices=None ):
         fig = plt.figure(figsize=(7,10))
         
         if profile_indices is None:
-            profile_indices=np.arange(0,self.dataset.dims['N_PROF'])
+            profile_indices=np.arange(0,self.dataset.dims['profile'])
             pass
 
         for ii in profile_indices:
-            prof_var = self.dataset[var].isel(N_PROF=ii)
-            prof_depth = self.dataset.depth.isel(N_PROF=ii)
+            prof_var = self.dataset[var].isel(profile=ii)
+            prof_depth = self.dataset.depth.isel(profile=ii)
             ax = plt.plot(prof_var, prof_depth)
             
         plt.gca().invert_yaxis()
@@ -77,21 +120,29 @@ def plot_profile(self, var:str, profile_indices=None ):
         plt.grid()
         return fig, ax
     
-    def plot_map(self, profile_indices=None):
+    def plot_map(self, profile_indices=None, var_str=None, depth_index=None):
         
         if profile_indices is None:
-            profile_indices=np.arange(0,self.dataset.dims['N_PROF'])
+            profile_indices=np.arange(0,self.dataset.dims['profile'])
+        
+        profiles=self.dataset.isel(profile=profile_indices)
         
-        profiles = self.dataset[['longitude','latitude']]
-        profiles=profiles.isel(N_PROF=profile_indices)
-        fig, ax = plot_util.geo_scatter(profiles.longitude.values,
-                                        profiles.latitude.values)
+        if var_str is None:
+            fig, ax = plot_util.geo_scatter(profiles.longitude.values,
+                                            profiles.latitude.values, s=5)
+        else:
+            print(profiles)
+            c = profiles[var_str].isel(level=depth_index)
+            fig, ax = plot_util.geo_scatter(profiles.longitude.values,
+                                            profiles.latitude.values,
+                                            c = c, s=5)
         
         return fig, ax
     
-    def plot_ts_diagram(self, profile_index, var_t='POTM_CORRECTED', var_s='PSAL_CORRECTED'):
+    def plot_ts_diagram(self, profile_index, var_t='potential_temperature', 
+                        var_s='practical_salinity'):
         
-        profile = self.dataset.isel(N_PROF=profile_index)
+        profile = self.dataset.isel(profile=profile_index)
         temperature = profile[var_t].values
         salinity = profile[var_s].values
         depth = profile.depth.values
diff --git a/coast/TIDEGAUGE.py b/coast/TIDEGAUGE.py
index 590dc7af..46efd8de 100755
--- a/coast/TIDEGAUGE.py
+++ b/coast/TIDEGAUGE.py
@@ -917,7 +917,7 @@ def plot_on_map_multiple(cls,tidegauge_list, color_var_str = None):
             fig, ax =  plot_util.geo_scatter(X, Y, title=title)
         else:
             fig, ax =  plot_util.geo_scatter(X, Y, title=title, c = C)
-            
+
         ax.set_xlim((min(X)-10, max(X)+10))
         ax.set_ylim((min(Y)-10, max(Y)+10))
         return fig, ax
@@ -1029,8 +1029,8 @@ def obs_operator(self, model, mod_var_name:str, time_interp = 'nearest',
         return
 
     def crps(self, model_object, model_var_name, obs_var_name:str='sea_level',
-         nh_radius: float = 20, cdf_type:str='empirical',
-         time_interp:str='linear', create_new_obj = True):
+             nh_radius: float = 20, time_interp:str='linear', 
+             create_new_obj = True):
         '''
         Comparison of observed variable to modelled using the Continuous
         Ranked Probability Score. This is done using this TIDEGAUGE object.
@@ -1065,11 +1065,11 @@ def crps(self, model_object, model_var_name, obs_var_name:str='sea_level',
 
         crps_list, n_model_pts, contains_land = crps_util.crps_sonf_fixed(
                                mod_var,
-                               self.dataset.longitude,
-                               self.dataset.latitude,
+                               self.dataset.longitude.values,
+                               self.dataset.latitude.values,
                                obs_var.values,
                                obs_var.time.values,
-                               nh_radius, cdf_type, time_interp )
+                               nh_radius, time_interp )
         if create_new_obj:
             new_object = TIDEGAUGE()
             new_dataset = self.dataset[['longitude','latitude','time']]
@@ -1252,8 +1252,11 @@ def find_high_and_low_water(self, var_str, method='comp',
         'comp' :: Find maxima by comparison with neighbouring values.
                   Uses scipy.signal.find_peaks. **kwargs passed to this routine
                   will be passed to scipy.signal.find_peaks.
-        DB NOTE: Currently only the 'comp' method is implemented. Future
-                 methods include linear interpolation and cublic splines.
+        'cubic':: Find the maxima using the roots of cubic spline.
+                  Uses scipy.interpolate.InterpolatedUnivariateSpline
+                  and scipy.signal.argrelmax. **kwargs are not activated.
+        NOTE: Currently only the 'comp' and 'cubic' methods implemented. Future
+                  methods include linear interpolation or refinements.
         '''
 
         x = self.dataset.time
diff --git a/coast/TRANSECT.py b/coast/TRANSECT.py
index b208f4c3..f1e10ee8 100644
--- a/coast/TRANSECT.py
+++ b/coast/TRANSECT.py
@@ -9,6 +9,7 @@
 from .logging_util import get_slug, debug, error
 import warnings
 import traceback
+from .logging_util import get_slug, debug, info, warn, error
 
 # =============================================================================
 # The TRANSECT module is a place for code related to transects only
@@ -286,7 +287,9 @@ def __init__(self, nemo_f: COAsT, point_A: tuple=None, point_B: tuple=None,
         super().__init__(nemo_f, point_A, point_B, y_indices, x_indices)
         
         
-    def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
+    def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT): 
+        #TODO the code here has become a little messy, could do with being 
+        # rewritten as in similar function in CONTOUR
         """
     
         Computes the flow through the transect at each segment and creates a new 
@@ -299,6 +302,10 @@ def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
         The latitude, longitude and the horizontal and vertical scale factors
         on the normal velocity points are also stored in the dataset.
         
+        If the time dependent cell thicknesses (e3) on the u and v grids are 
+        present in the nemo_u and nemo_v datasets they will be used, if they 
+        are not then the initial cell thicknesses (e3_0) will be used.
+        
         parameters
         ----------
         nemo_u : Nemo object on the u-grid containing the i-component velocities
@@ -306,6 +313,9 @@ def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
         
         """
         debug(f"Computing flow across the transect for {get_slug(self)}")
+        
+        # compute transports flag; set to false if suitable e3 not found
+        compute_transports = True
                 
         # subset the u and v datasets 
         da_y_ind = xr.DataArray( self.y_ind, dims=['r_dim'] )
@@ -313,12 +323,32 @@ def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
         u_ds = nemo_u.dataset.isel(y_dim = da_y_ind, x_dim = da_x_ind)
         v_ds = nemo_v.dataset.isel(y_dim = da_y_ind, x_dim = da_x_ind)
         
+        # use time varying if e3 is present, if not default to e3_0
+        if 'e3' not in u_ds.data_vars:
+            if 'e3_0' not in u_ds.data_vars:
+                warn("e3 not found, transports will not be calculated")
+                compute_transports = False
+            else:
+                u_ds['e3'] = u_ds.e3_0.broadcast_like(u_ds.u_velocity) 
+        if 'e3' not in v_ds.data_vars:
+            if 'e3_0' not in v_ds.data_vars:
+                warn("e3 not found, transports will not be calculated")
+                compute_transports = False
+            else:
+                v_ds['e3'] = v_ds.e3_0.broadcast_like(v_ds.v_velocity) 
+        
         # If there is no time dimension, add one. This is so
         # indexing can assume a time dimension exists
+        droptime=False
         if 't_dim' not in u_ds.dims:
-            u_ds = u_ds.expand_dims(dim={'t_dim':1},axis=0)
+            u_ds["u_velocity"] = u_ds.u_velocity.expand_dims(dim={'t_dim':1},axis=0)
+            if compute_transports:
+                u_ds["e3"] = u_ds.e3.expand_dims(dim={'t_dim':1},axis=0)
+            droptime=True
         if 't_dim' not in v_ds.dims:
-            v_ds = v_ds.expand_dims(dim={'t_dim':1},axis=0)
+            v_ds["v_velocity"] = v_ds.v_velocity.expand_dims(dim={'t_dim':1},axis=0)
+            if compute_transports:
+                v_ds["e3"] = v_ds.e3.expand_dims(dim={'t_dim':1},axis=0)
         
         velocity = np.ma.zeros( (u_ds.t_dim.size, u_ds.z_dim.size, u_ds.r_dim.size-1) )
         vol_transport = np.ma.zeros( (u_ds.t_dim.size, u_ds.z_dim.size, u_ds.r_dim.size-1) )
@@ -327,7 +357,7 @@ def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
         longitude = np.ma.zeros( (u_ds.r_dim.size-1) )
         e1 = np.ma.zeros( (u_ds.r_dim.size-1) )
         e2 = np.ma.zeros( (u_ds.r_dim.size-1) )
-        e3_0 = np.ma.zeros( (u_ds.z_dim.size, u_ds.r_dim.size-1) )
+        e3 = np.ma.zeros( (u_ds.t_dim.size, u_ds.z_dim.size, u_ds.r_dim.size-1) )
         
         # Find the indices where the derivative of the transect in the north, south, east and west
         # directions are positive.
@@ -339,62 +369,79 @@ def calc_flow_across_transect(self, nemo_u: COAsT, nemo_v: COAsT):
         dr_w = dr_w[~np.isnan(dr_w)].astype(int)
         
         # u flux (+ in)
-        velocity[:,:,dr_n] = u_ds.vozocrtx.to_masked_array()[:,:,dr_n+1]
-        vol_transport[:,:,dr_n] = ( velocity[:,:,dr_n] * u_ds.e2.to_masked_array()[dr_n+1] *
-                                          u_ds.e3_0.to_masked_array()[:,dr_n+1] )
+        velocity[:,:,dr_n] = u_ds.u_velocity.to_masked_array()[:,:,dr_n+1]
+        if compute_transports:
+            vol_transport[:,:,dr_n] = ( velocity[:,:,dr_n] * u_ds.e2.to_masked_array()[dr_n+1] *
+                                          u_ds.e3.to_masked_array()[:,:,dr_n+1] )
+            e3[:,:,dr_n] = u_ds.e3.values[:,:,dr_n+1]
         depth_0[:,dr_n] = u_ds.depth_0.to_masked_array()[:,dr_n+1]
         latitude[dr_n] = u_ds.latitude.values[dr_n+1]
         longitude[dr_n] = u_ds.longitude.values[dr_n+1]
         e1[dr_n] = u_ds.e1.values[dr_n+1]
         e2[dr_n] = u_ds.e2.values[dr_n+1]
-        e3_0[:,dr_n] = u_ds.e3_0.values[:,dr_n+1]
-        
+               
         # v flux (- in) 
-        velocity[:,:,dr_e] = - v_ds.vomecrty.to_masked_array()[:,:,dr_e+1]
-        vol_transport[:,:,dr_e] = ( velocity[:,:,dr_e] * v_ds.e1.to_masked_array()[dr_e+1] *
-                                  v_ds.e3_0.to_masked_array()[:,dr_e+1] )
+        velocity[:,:,dr_e] = - v_ds.v_velocity.to_masked_array()[:,:,dr_e+1]
+        if compute_transports:
+            vol_transport[:,:,dr_e] = ( velocity[:,:,dr_e] * v_ds.e1.to_masked_array()[dr_e+1] *
+                                  v_ds.e3.to_masked_array()[:,:,dr_e+1] )
+            e3[:,:,dr_e] = v_ds.e3.values[:,:,dr_e+1]
         depth_0[:,dr_e] = v_ds.depth_0.to_masked_array()[:,dr_e+1]
         latitude[dr_e] = v_ds.latitude.values[dr_e+1]
         longitude[dr_e] = v_ds.longitude.values[dr_e+1]
         e1[dr_e] = v_ds.e1.values[dr_e+1]
         e2[dr_e] = v_ds.e2.values[dr_e+1]
-        e3_0[:,dr_e] = v_ds.e3_0.values[:,dr_e+1]
-        
+               
         # v flux (+ in)
-        velocity[:,:,dr_w] = v_ds.vomecrty.to_masked_array()[:,:,dr_w]
-        vol_transport[:,:,dr_w] = ( velocity[:,:,dr_w] * v_ds.e1.to_masked_array()[dr_w] *
-                                  v_ds.e3_0.to_masked_array()[:,dr_w] )
+        velocity[:,:,dr_w] = v_ds.v_velocity.to_masked_array()[:,:,dr_w]
+        if compute_transports:
+            vol_transport[:,:,dr_w] = ( velocity[:,:,dr_w] * v_ds.e1.to_masked_array()[dr_w] *
+                                  v_ds.e3.to_masked_array()[:,:,dr_w] )
+            e3[:,:,dr_w] = v_ds.e3.values[:,:,dr_w]
         depth_0[:,dr_w] = v_ds.depth_0.to_masked_array()[:,dr_w]
         latitude[dr_w] = v_ds.latitude.values[dr_w]
         longitude[dr_w] = v_ds.longitude.values[dr_w]
         e1[dr_w] = v_ds.e1.values[dr_w]
         e2[dr_w] = v_ds.e2.values[dr_w]
-        e3_0[:,dr_w] = v_ds.e3_0.values[:,dr_w]
-           
-        # Add DataArrays to dataset           
-        self.data_cross_tran_flow['normal_velocities'] = xr.DataArray( np.squeeze(velocity), 
-                    coords={'time': (('t_dim'), u_ds.time.values),'depth_0': (('z_dim','r_dim'), depth_0)
+                   
+        # Add DataArrays to dataset 
+        if droptime:
+            self.data_cross_tran_flow['normal_velocities'] = xr.DataArray( velocity.squeeze(), 
+                    coords={'depth_0': (('z_dim','r_dim'), depth_0) 
                             ,'latitude': (('r_dim'), latitude), 'longitude': (('r_dim'), longitude) },
-                    dims=['t_dim', 'z_dim', 'r_dim'] )            
-        self.data_cross_tran_flow['normal_transports'] \
-                    = xr.DataArray( np.squeeze(np.sum(vol_transport, axis=1)) / 1000000.,
-                    coords={'time': (('t_dim'), u_ds.time.values)
-                            ,'latitude': (('r_dim'), latitude), 'longitude': (('r_dim'), longitude)},
-                    dims=['t_dim', 'r_dim'] )          
+                    dims=['z_dim', 'r_dim'] ) 
+            if compute_transports:
+                self.data_cross_tran_flow['normal_transports'] \
+                    = xr.DataArray( np.sum(vol_transport.squeeze(), axis=0) / 1000000., \
+                    coords={'latitude': (('r_dim'), latitude), 'longitude': (('r_dim'), longitude)},\
+                    dims=['r_dim'] ).squeeze() 
+        else:
+            self.data_cross_tran_flow['normal_velocities'] = xr.DataArray( velocity, 
+                        coords={'time': (('t_dim'), u_ds.time.values),'depth_0': (('z_dim','r_dim'), depth_0)
+                                ,'latitude': (('r_dim'), latitude), 'longitude': (('r_dim'), longitude) },
+                        dims=['t_dim', 'z_dim', 'r_dim'] )   
+            if compute_transports:
+                self.data_cross_tran_flow['normal_transports'] \
+                        = xr.DataArray( np.sum(vol_transport, axis=1) / 1000000.,
+                        coords={'time': (('t_dim'), u_ds.time.values)
+                                ,'latitude': (('r_dim'), latitude), 'longitude': (('r_dim'), longitude)},
+                        dims=['t_dim', 'r_dim'] ).squeeze()          
         self.data_cross_tran_flow['e1'] = xr.DataArray( e1, dims=['r_dim'] ) 
         self.data_cross_tran_flow['e2'] = xr.DataArray( e2, dims=['r_dim'] ) 
-        self.data_cross_tran_flow['e3_0'] = xr.DataArray( e3_0, dims=['z_dim','r_dim'] ) 
+        if compute_transports:
+            self.data_cross_tran_flow['e3'] = xr.DataArray( e3, dims=['t_dim','z_dim','r_dim'] ) 
         # DataArray attributes   
         self.data_cross_tran_flow.normal_velocities.attrs['units'] = 'm/s'
         self.data_cross_tran_flow.normal_velocities.attrs['standard_name'] = 'velocity across the transect'
-        self.data_cross_tran_flow.normal_velocities.attrs['long_name'] = 'velocity across the transect defined on the normal velocity grid points'    
-        self.data_cross_tran_flow.normal_transports.attrs['units'] = 'Sv'
-        self.data_cross_tran_flow.normal_transports.attrs['standard_name'] = 'depth integrated volume transport across transect'
-        self.data_cross_tran_flow.normal_transports.attrs['long_name'] = 'depth integrated volume transport across the transect defined on the normal velocity grid points' 
+        self.data_cross_tran_flow.normal_velocities.attrs['long_name'] = 'velocity across the transect defined on the normal velocity grid points'  
+        if compute_transports:
+            self.data_cross_tran_flow.normal_transports.attrs['units'] = 'Sv'
+            self.data_cross_tran_flow.normal_transports.attrs['standard_name'] = 'depth integrated volume transport across transect'
+            self.data_cross_tran_flow.normal_transports.attrs['long_name'] = 'depth integrated volume transport across the transect defined on the normal velocity grid points' 
         self.data_cross_tran_flow.depth_0.attrs['units'] = 'm'
         self.data_cross_tran_flow.depth_0.attrs['standard_name'] = 'depth'
         self.data_cross_tran_flow.depth_0.attrs['long_name'] = 'Initial depth at time zero defined at the normal velocity grid points'
-                        
+        self.data_cross_tran_flow = self.data_cross_tran_flow.squeeze()
   
     def __pressure_grad_fpoint(self, ds_T, ds_T_j1, ds_T_i1, ds_T_j1i1, r_ind, velocity_component):
         """
diff --git a/coast/__init__.py b/coast/__init__.py
index e605b2aa..ca5dc147 100644
--- a/coast/__init__.py
+++ b/coast/__init__.py
@@ -4,10 +4,12 @@
 from .TRANSECT import Transect, Transect_f, Transect_t
 from .ALTIMETRY import ALTIMETRY
 from .OBSERVATION import OBSERVATION
-from .CDF import CDF
+from .DISTRIBUTION import DISTRIBUTION
 from .INTERNALTIDE import INTERNALTIDE
 from .TIDEGAUGE import TIDEGAUGE
 from .PROFILE import PROFILE
+from .CLIMATOLOGY import CLIMATOLOGY
+from .MASK_MAKER import MASK_MAKER
 from . import logging_util
 from . import general_utils
 from . import plot_util
diff --git a/coast/crps_util.py b/coast/crps_util.py
index 1d8c12d0..b180b8ae 100755
--- a/coast/crps_util.py
+++ b/coast/crps_util.py
@@ -1,7 +1,6 @@
 '''
 Python definitions used to aid in the calculation of Continuous Ranked
 Probability Score.
-
 *Methods Overview*
     -> crps_sonf_fixed(): Single obs neighbourhood forecast CRPS for fixed obs
     -> crps_song_moving(): Same as above for moving obs
@@ -9,16 +8,131 @@
 
 import numpy as np
 import xarray as xr
-from .CDF import CDF
+from . import general_utils
+
+def crps_empirical(sample, obs):
+        """Calculates CRPS for a single observations against a sample of values.
+        This sample of values may be an ensemble of model forecasts or a model
+        neighbourhood. This is a comparison of a Heaviside function defined by
+        the observation value and an Empirical Distribution Function (EDF)
+        defined by the sample of values. This sample is sorted to create the 
+        EDF.
+        
+        The calculation method is that outlined by Hersbach et al. (2000). 
+        Each member of a supplied sample is weighted equally.
+        
+        Args:
+            sample (array): Array of points (ensemble or neighbourhood)
+            xa (float): A single 'observation' value which to compare against
+                        sample CDF.
+        Returns:
+            A single CRPS value.
+        """
+
+        def calc(alpha, beta, p):
+            return alpha * p**2 + beta*(1 - p)**2
+        
+        xa = float(obs)
+        crps_integral = 0
+        sample = np.array(sample)
+        
+        if all(np.isnan(sample)) or np.isnan(obs):
+            return np.nan
+        
+        sample = sample[~np.isnan(sample)]
+        sample = np.sort(sample)
+        sample_size = len(sample)
+        
+        alpha = np.zeros(sample_size-1)
+        beta= np.zeros(sample_size-1)
+        # sample[1:] = upper bounds, and vice versa
+        
+        tmp = sample[1:] - sample[:-1]
+        tmp_logic = sample[1:]<xa
+        alpha[tmp_logic] = tmp[tmp_logic]
+                               
+        tmp_logic = sample[:-1]>xa
+        beta[tmp_logic] = tmp[tmp_logic]
+        
+        tmp_logic = ( sample[1:]>xa )*( sample[:-1]<xa )
+        tmp = xa - sample[:-1]
+        alpha[tmp_logic] = tmp[tmp_logic]
+        tmp = sample[1:] - xa
+        beta[tmp_logic] = tmp[tmp_logic]
+        
+        p = np.arange(1,sample_size)/sample_size
+        c = alpha*p**2 + beta*(1-p)**2
+        crps_integral = np.sum(c)
+        
+        # Intervals 0 and N, where p = 0 and 1 respectively
+        if xa < sample[0]:
+            p=0
+            alpha = 0
+            beta = sample[0] - xa
+            crps = calc(alpha, beta, p)
+            crps_integral += crps
+        elif xa > sample[-1]:
+            p=1
+            alpha = xa - sample[-1]
+            beta = 0
+            crps = calc(alpha, beta, p)
+            crps_integral += crps
+            
+        crps = crps_integral
+            
+        return crps
+    
+def crps_empirical_loop(sample, obs):
+        """Like crps_empirical, however a loop is used instead of numpy 
+        boolean indexing. For large samples, will be slower but consume less
+        memory.
+        """
+        
+        def calc(alpha, beta, p):
+            return alpha * p**2 + beta*(1 - p)**2
+        crps_integral = 0
+        sample = np.array(sample)
+        sample = sample[~np.isnan(sample)]
+        sample = np.sort(sample)
+        sample_size = len(sample)
+        
+        # All intervals within range of the sample distribution
+        for ii in range(0, sample_size-1):
+            p = (ii+1)/sample_size
+            if obs > sample[ii+1]:
+                alpha = sample[ii+1] - sample[ii]
+                beta = 0
+            elif obs < sample[ii]:
+                alpha = 0
+                beta = sample[ii+1] - sample[ii]
+            else:
+                alpha = obs - sample[ii]
+                beta = sample[ii+1] - obs
+            crps = calc(alpha, beta, p)
+            crps_integral += crps
+        # Intervals 0 and N, where p = 0 and 1 respectively
+        if obs < sample[0]:
+            p=0
+            alpha = 0
+            beta = sample[0] - obs
+            crps = calc(alpha, beta, p)
+            crps_integral += crps
+        elif obs > sample[-1]:
+            p=1
+            alpha = obs - sample[-1]
+            beta = 0
+            crps = calc(alpha, beta, p)
+            crps_integral += crps
+            
+        return crps_integral
 
 def crps_sonf_fixed( mod_array, obs_lon, obs_lat, obs_var, obs_time, 
-                      nh_radius: float, cdf_type:str, time_interp:str
+                      nh_radius: float, time_interp:str,
     ):
     '''
     Handles the calculation of single-observation neighbourhood forecast CRPS
     for a time series at a fixed observation location. Differs from 
     crps_sonf_moving in that it only need calculate a model neighbourhood once.
-
     Parameters
     ----------
     mod_array   : (xarray DataArray) DataArray from a Model Dataset
@@ -27,10 +141,7 @@ def crps_sonf_fixed( mod_array, obs_lon, obs_lat, obs_var, obs_time,
     obs_var     : (array) of floatArray of variable values, e.g time series
     obs_time    : (array) of datetimeArray of times, corresponding to obs_var
     nh_radius   : (float) Neighbourhood radius in km
-    cdf_type    : (str) Type of CDF to use for model data. Either 'empirical' 
-                   or 'theoretical'.
     time_interp : (str) Type of time interpolation to use
-
     Returns
     -------
     crps_list     : Array of CRPS values
@@ -38,7 +149,6 @@ def crps_sonf_fixed( mod_array, obs_lon, obs_lat, obs_var, obs_time,
                     each CRPS value
     contains_land : Array of bools indicating where a model neighbourhood 
                     contained land.
-
     '''
 
     # Define output arrays
@@ -51,10 +161,10 @@ def crps_sonf_fixed( mod_array, obs_lon, obs_lat, obs_var, obs_time,
     neighbourhood_indices = np.arange(0,n_neighbourhoods)
     
     # Get model neighbourhood subset using specified method
-    subset_ind = subset_indices_by_distance(
-                     mod_array.longitude, mod_array.latitude, 
+    subset_ind = general_utils.subset_indices_by_distance(
+                     mod_array.longitude.values, mod_array.latitude.values, 
                      obs_lon, obs_lat, nh_radius)
-    mod_subset = mod_array.isel(y_dim = subset_ind[0],
+    mod_subset = mod_array.isel(  y_dim = subset_ind[0],
                                   x_dim = subset_ind[1])
     mod_subset = mod_subset.swap_dims({'t_dim':'time'})
     
@@ -76,35 +186,32 @@ def crps_sonf_fixed( mod_array, obs_lon, obs_lat, obs_var, obs_time,
                 if all(np.isnan(mod_subset_time)):
                     pass
                 else:
-                    # Create model and observation CDF objects
-                    mod_cdf = CDF(mod_subset_time, cdf_type = cdf_type)
-                
                     # Calculate CRPS and put into output array
-                    crps_list[ii] = mod_cdf.crps_fast(obs_var[ii])
+                    crps_list[ii] = crps_empirical(mod_subset_time.values, 
+                                                   obs_var[ii])
                     n_model_pts[ii] = int(mod_subset.shape[0])
 
     return crps_list, n_model_pts, contains_land
 
 def crps_sonf_moving( mod_array, obs_lon, obs_lat, obs_var, obs_time, 
-                      nh_radius: float, cdf_type:str, time_interp:str
+                      nh_radius: float, time_interp:str,
+                      obs_batch=10
     ):
     '''
     Handles the calculation of single-observation neighbourhood forecast CRPS
     for a moving observation instrument. Differs from crps_sonf_fixed in that 
-    latitude and longitude are arrays of locations.
-
+    latitude and longitude are arrays of locations. Mod_array must contain
+    dimensions x_dim, y_dim and t_dim and coordinates longitude, latitude,
+    time.
     Parameters
     ----------
     mod_array   : (xarray DataArray) DataArray from a Model Dataset
-    obs_lon     : (array) Longitudes of fixed observation point
-    obs_lat     : (array) Latitudes of fixed observation point
-    obs_var     : (array) of floatArray of variable values, e.g time series
-    obs_time    : (array) of datetimeArray of times, corresponding to obs_var
+    obs_lon     : (1Darray) Longitudes of fixed observation point
+    obs_lat     : (1Darray) Latitudes of fixed observation point
+    obs_var     : (1Darray) of floatArray of variable values, e.g time series
+    obs_time    : (1Darray) of datetimeArray of times, corresponding to obs_var
     nh_radius   : (float) Neighbourhood radius in km
-    cdf_type    : (str) Type of CDF to use for model data. Either 'empirical' 
-                   or 'theoretical'.
     time_interp : (str) Type of time interpolation to use
-
     Returns
     -------
     crps_list     : Array of CRPS values
@@ -119,20 +226,17 @@ def crps_sonf_moving( mod_array, obs_lon, obs_lat, obs_var, obs_time,
     crps_list     = np.zeros( n_neighbourhoods )*np.nan
     n_model_pts   = np.zeros( n_neighbourhoods )*np.nan
     contains_land = np.zeros( n_neighbourhoods , dtype=bool)
-
     # Loop over neighbourhoods
     neighbourhood_indices = np.arange(0,n_neighbourhoods)
     for ii in neighbourhood_indices:
-        
         # Neighbourhood centre
         cntr_lon = obs_lon[ii]
         cntr_lat = obs_lat[ii]
     
         # Get model neighbourhood subset using specified method
-        subset_ind = subset_indices_by_distance(
+        subset_ind = general_utils.subset_indices_by_distance(
                          mod_array.longitude, mod_array.latitude, 
                          cntr_lon, cntr_lat, nh_radius)
-        
         # Check that the model neighbourhood contains points
         if subset_ind[0].shape[0] == 0 or subset_ind[1].shape[0] == 0:
             crps_list[ii] = np.nan
@@ -144,7 +248,7 @@ def crps_sonf_moving( mod_array, obs_lon, obs_lat, obs_var, obs_time,
             mod_subset = mod_subset.interp(
                              time = obs_time[ii], method = time_interp,
                              kwargs={'fill_value':'extrapolate'})
-            
+
             #Check if neighbourhood contains a land value (TODO:mask)
             if any(np.isnan(mod_subset)):
                 contains_land[ii] = True
@@ -152,65 +256,9 @@ def crps_sonf_moving( mod_array, obs_lon, obs_lat, obs_var, obs_time,
             if all(np.isnan(mod_subset)):
                 pass
             else:
-                # Create model and observation CDF objects
-                mod_cdf = CDF(mod_subset, cdf_type = cdf_type)
-            
                 # Calculate CRPS and put into output array
-                crps_list[ii] = mod_cdf.crps_fast(obs_var[ii])
+                crps_list[ii] = crps_empirical(mod_subset, obs_var[ii])
                 n_model_pts[ii] = int(mod_subset.shape[0])
+        crps_mean = np.nanmean(crps_list)
                 
-    return crps_list, n_model_pts, contains_land
-
-def subset_indices_by_distance(
-        longitude, latitude, centre_lon: float, centre_lat: float, 
-        radius: float
-    ):
-    """
-    This method returns a `tuple` of indices within the `radius` of the lon/lat point given by the user.
-
-    Distance is calculated as haversine - see `self.calculate_haversine_distance`
-
-    :param centre_lon: The longitude of the users central point
-    :param centre_lat: The latitude of the users central point
-    :param radius: The haversine distance (in km) from the central point
-    :return: All indices in a `tuple` with the haversine distance of the central point
-    """
-
-    # Calculate the distances between every model point and the specified
-    # centre. Calls another routine dist_haversine.
-
-    dist = calculate_haversine_distance(centre_lon, centre_lat, 
-                                        longitude, latitude)
-    indices_bool = dist < radius
-    indices = np.where(indices_bool.compute())
-
-    return xr.DataArray(indices[0]), xr.DataArray(indices[1])
-
-def calculate_haversine_distance(lon1, lat1, lon2, lat2, r = 6371.00717):
-    '''
-    # Estimation of geographical distance using the Haversine function.
-    # Input can be single values or 1D arrays of locations. This
-    # does NOT create a distance matrix but outputs another 1D array.
-    # This works for either location vectors of equal length OR a single loc
-    # and an arbitrary length location vector.
-    #
-    # lon1, lat1 :: Location(s) 1.
-    # lon2, lat2 :: Location(s) 2.
-    '''
-
-    # Convert to radians for calculations
-    lon1 = xr.ufuncs.deg2rad(lon1)
-    lat1 = xr.ufuncs.deg2rad(lat1)
-    lon2 = xr.ufuncs.deg2rad(lon2)
-    lat2 = xr.ufuncs.deg2rad(lat2)
-
-    # Latitude and longitude differences
-    dlat = (lat2 - lat1) / 2
-    dlon = (lon2 - lon1) / 2
-
-    # Haversine function.
-    distance = xr.ufuncs.sin(dlat) ** 2 + xr.ufuncs.cos(lat1) *  \
-        xr.ufuncs.cos(lat2) * xr.ufuncs.sin(dlon) ** 2
-    distance = 2 * r * xr.ufuncs.arcsin(xr.ufuncs.sqrt(distance))
-
-    return distance
+    return crps_list, n_model_pts, contains_land
\ No newline at end of file
diff --git a/coast/eof.py b/coast/eof.py
index cdd80c4e..c21bbcda 100755
--- a/coast/eof.py
+++ b/coast/eof.py
@@ -195,7 +195,7 @@ def _compute(signal, full_matrices, active_ind, number_points):
     '''
     P, D, Q = linalg.svd(signal, full_matrices=full_matrices)   
     mode_count = P.shape[-1]     
-    EOFs = np.zeros( (number_points, mode_count), dtype=P.dtype)
+    EOFs = np.full( (number_points, mode_count), np.nan, dtype=P.dtype)
     EOFs[active_ind,:] = P 
     
     # Calculate variance explained
diff --git a/coast/general_utils.py b/coast/general_utils.py
index ad9001a0..4716b937 100755
--- a/coast/general_utils.py
+++ b/coast/general_utils.py
@@ -35,9 +35,6 @@ def subset_indices_by_distance_BT(longitude, latitude, centre_lon, centre_lat,
     If longitude is 2D:
         Returns arrays of x and y indices per central location.
         ind_y corresponds to row indices of the original input arrays.
-        
-    DB: This is really slow at the moment. But the bottleneck seems to be right
-    at the start somewhere.. Best stick to subset_indices_by_distance for now.
     """
     # change inputs to numpy
     longitude = np.array(longitude)
@@ -109,7 +106,6 @@ def subset_indices_by_distance(
 
     # Calculate the distances between every model point and the specified
     # centre. Calls another routine dist_haversine.
-
     dist = calculate_haversine_distance(centre_lon, centre_lat, 
                                         longitude, latitude)
     indices_bool = dist < radius
diff --git a/coast/plot_util.py b/coast/plot_util.py
index ea6ee575..98812d23 100755
--- a/coast/plot_util.py
+++ b/coast/plot_util.py
@@ -10,6 +10,144 @@
 from .logging_util import get_slug, debug, info, warn, error
 import numpy as np
 
+def r2_lin(x, y, fit):
+    '''For calculating r-squared of a linear fit. Fit should be a python polyfit
+    object'''
+
+    fity = fit(x)
+    diff = (y - fity)**2
+    ybar = np.nanmean(y)
+    ymybar = (y - ybar)**2
+
+    SStot = np.nansum(ymybar)
+    SSres = np.nansum(diff)
+
+    R2 = 1 - SSres/SStot
+
+    return R2
+
+def scatter_with_fit(x,y,s=10,c='k',yex=True, dofit=True):
+    ''' Does a scatter plot with a linear fit. Will also draw y=x for
+    comparison.
+    
+    Parameters
+    ----------
+    x     : (array) Values for the x-axis
+    y     : (array) Values for the y-axis
+    s     : (float or array) Marker size(s)
+    c     : (float or array) Marker colour(s)
+    yex   : (bool) True to plot y=x
+    dofit : (bool) True to calculate and plot linear fit
+    
+    Returns
+    -------
+    Figure and axis objects for further customisation
+    
+    Example Useage
+    -------
+    x = np.arange(0,50)
+    y = np.arange(0,50)/1.5
+    f,a = scatter_with_fit(x,y)
+    a.set_title('Example scatter with fit')
+    a.set_xlabel('Example x axis')
+    a.set_ylabel('Example y axis')
+    '''
+    
+    fig = plt.figure()
+    ax = fig.add_subplot(1,1,1)
+    
+    x = np.ma.masked_invalid(x)
+    y = np.ma.masked_invalid(y)
+    combined_mask = np.ma.mask_or(x.mask,y.mask)
+    x.mask = combined_mask
+    y.mask = combined_mask
+    
+    xmax = np.ma.max(x)
+    xmin = np.ma.min(x)
+    ymax = np.ma.max(y)
+    ymin = np.ma.min(y)
+    axmax0 = np.max([xmax, ymax])
+    axmin0 = np.min([xmin, ymin])
+    axmin = axmin0 - 0.1*np.abs(axmax0-axmin0)
+    axmax = axmax0 + 0.1*np.abs(axmax0-axmin0)
+    
+    if yex:
+        lineX = [axmin,axmax]
+        fityx = np.poly1d([1,0])
+        ax.plot(lineX, fityx(lineX),c=[0.5,0.5,0.5],linewidth=1)
+    
+    sca = ax.scatter(x,y, c=c, s=s)
+    
+    if dofit:
+        lineX = [axmin,axmax]
+         #Calculate data fit and cast to poly1d object
+        fit_tmp = np.ma.polyfit(x, y, 1)
+        fit = np.poly1d( fit_tmp )
+        ax.plot(lineX, fit(lineX),c=[1,128/255,0],linewidth=1.5)
+        r2 = r2_lin(x, y, fit)
+    
+    ax.set_xlim(axmin, axmax)
+    ax.set_ylim(axmin, axmax)
+    ax.set_aspect('equal', adjustable='box')
+    ax.grid()
+    
+    if dofit:
+        ax.text(0.4,0.125,'{} {:03.2f} {} {:03.2f}'.format('y =',
+                            fit_tmp[0],'x +',fit_tmp[1]),
+                            transform=ax.transAxes)
+        ax.text(0.4,0.05,'{} {:03.2f} '.format('$R^2$ =',
+                    r2),transform=ax.transAxes)
+    
+    return fig, ax
+
+def create_geo_axes(lonbounds, latbounds):
+    '''
+    A routine for creating an axis for any geographical plot. Within the
+    specified longitude and latitude bounds, a map will be drawn up using
+    cartopy. Any type of matplotlib plot can then be added to this figure.
+    For example:
+        
+    Example Useage
+    #############
+    
+        f,a = create_geo_axes(lonbounds, latbounds)
+        sca = a.scatter(stats.longitude, stats.latitude, c=stats.corr, 
+                        vmin=.75, vmax=1,
+                        edgecolors='k', linewidths=.5, zorder=100)
+        f.colorbar(sca)
+        a.set_title('SSH correlations \n Monthly PSMSL tide gauge vs CO9_AMM15p0', 
+                    fontsize=9)
+        
+    * Note: For scatter plots, it is useful to set zorder = 100 (or similar
+            positive number)
+    '''
+
+    import cartopy.crs as ccrs  # mapping plots
+    from cartopy.feature import NaturalEarthFeature
+        
+    # If no figure or ax is provided, create a new one
+    fig = plt.figure(1)
+    fig.clf()
+    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
+        
+    coast = NaturalEarthFeature(category='physical', facecolor=[0.9,0.9,0.9], name='coastline',
+                            scale='50m')
+    ax.add_feature(coast, edgecolor='gray')
+    #ax.coastlines(facecolor=[0.8,0.8,0.8])
+    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
+                      linewidth=0.5, color='gray', linestyle='-')
+    gl.top_labels = False
+    gl.bottom_labels = True
+    gl.right_labels = False
+    gl.left_labels = True
+    
+    ax.set_xlim(lonbounds[0], lonbounds[1])
+    ax.set_ylim(latbounds[0], latbounds[1])
+    ax.set_aspect('auto')
+
+    plt.show()
+    return fig, ax
+
 def ts_diagram(temperature, salinity, depth):
     
     fig = plt.figure(figsize = (10,7))
@@ -86,4 +224,32 @@ def geo_scatter(longitude, latitude, c=None, s = None,
     ax.set_aspect('auto')
 
     plt.show()
-    return fig, ax
\ No newline at end of file
+    return fig, ax
+
+def determine_colorbar_extension(color_data, vmin, vmax):
+    ''' Can be used to automatically determine settings for colorbar 
+    extension arrows. Color_data is the data used for the colormap, vmin
+    and vmax are the colorbar limits. Will output a string: "both", "max",
+    "min" or "neither", which can be inserted straight into a call to
+    matplotlib.pyplot.colorbar().
+    '''
+    extend_max = np.nanmax(color_data) > vmax
+    extend_min = np.nanmin(color_data) < vmin
+
+    if extend_max and extend_min: return "both"
+    elif extend_max and not extend_min: return 'max'
+    elif not extend_max and extend_min: return 'min'
+    else: return 'neither'
+
+def determine_clim_by_standard_deviation(color_data, n_std_dev=2.5):
+    ''' Automatically determine color limits based on number of standard
+    deviations from the mean of the color data (color_data). Useful if there
+    are outliers in the data causing difficulties in distinguishing most of 
+    the data. Outputs vmin and vmax which can be passed to plotting routine
+    or plt.clim().
+    '''
+    color_data_mean = np.nanmean(color_data)
+    color_data_std = np.nanstd(color_data)
+    vmin = color_data_mean - n_std_dev*color_data_std
+    vmax = color_data_mean + n_std_dev*color_data_std
+    return vmin, vmax 
\ No newline at end of file
diff --git a/coast/stats_util.py b/coast/stats_util.py
index 4971aaca..618661ae 100755
--- a/coast/stats_util.py
+++ b/coast/stats_util.py
@@ -12,61 +12,169 @@
 from .logging_util import get_slug, debug, info, warn, error
 import scipy
 
+def quadratic_spline_roots(spl):
+    """
+    A custom function for the roots of a quadratic spline. Cleverness found at
+    https://stackoverflow.com/questions/50371298/find-maximum-minimum-of-a-1d-interpolated-function
+    Used in find_maxima().
+
+    Example usage:
+    see example_scripts/tidegauge_tutorial.py
+    """
+    roots = []
+    knots = spl.get_knots()
+    for a, b in zip(knots[:-1], knots[1:]):
+        u, v, w = spl(a), spl((a+b)/2), spl(b)
+        t = np.roots([u+w-2*v, w-u, 2*v])
+        t = t[np.isreal(t) & (np.abs(t) <= 1)]
+        roots.extend(t*(b-a)/2 + (b+a)/2)
+    return np.sort(roots)
+
 def find_maxima(x, y, method='comp', **kwargs):
     '''
     Finds maxima of a time series y. Returns maximum values of y (e.g heights)
-    and corresponding values of x (e.g. times). 
+    and corresponding values of x (e.g. times).
     **kwargs are dependent on method.
-    
+
         Methods:
         'comp' :: Find maxima by comparison with neighbouring values.
                   Uses scipy.signal.find_peaks. **kwargs passed to this routine
                   will be passed to scipy.signal.find_peaks.
-        DB NOTE: Currently only the 'comp' method is implemented. Future
-                 methods include linear interpolation and cublic splines.
+        'cubic' :: Find the maxima and minima by fitting a cubic spline and
+                    finding the roots of its derivative.
+                    Expect input as xr.DataArrays
+        DB NOTE: Currently only the 'comp' and 'cubic' method are implemented.
+                Future methods include linear interpolation.
+
+        JP NOTE: Cubic method:
+            i) has ineligent fix for NaNs,
+
+    Example usage:
+    see example_scripts/tidegauge_tutorial.py
     '''
+    
+    if method == 'cubic':
+        if (type(x) != xr.DataArray) or (type(y) != xr.DataArray):
+            msg = 'With method {} require input to be type: xr.DataArray' \
+                + ' not {} and {}\n Reset method as comp'
+            print(msg.format( method, type(x), type(y) ))
+            method = 'comp'
+
     if method == 'comp':
         peaks, props = scipy.signal.find_peaks(y, **kwargs)
         return x[peaks], y[peaks]
 
+    if method == 'cubic':
+        """
+        Cleverness found at
+        https://stackoverflow.com/questions/50371298/find-maximum-minimum-of-a-1d-interpolated-function
+
+        Find the extrema on a cubic spline fitted to 1d array. E.g:
+        # Some data
+        x_axis = xr.DataArray([ 2.14414414,  2.15270826,  2.16127238,  2.1698365 ,  2.17840062, 2.18696474,  2.19552886,  2.20409298,  2.2126571 ,  2.22122122])
+        y_axis = xr.DataArray([ 0.67958442,  0.89628424,  0.78904004,  3.93404167,  6.46422317, 6.40459954,  3.80216674,  0.69641825,  0.89675386,  0.64274198])
+        # Fit cubic spline
+        f = scipy.interpolate.interp1d(x_axis, y_axis, kind = 'cubic')
+        x_new = np.linspace(x_axis[0], x_axis[-1],100)
+        cr_pts, cr_vals = stats_utils.find_maxima(x_axis, y_axis, method='cubic')
+        fig = plt.subplots()
+        plt.plot(x_axis, y_axis, 'r+') # The fitted spline
+        plt.plot(x_new, f(x_new)) # The fitted spline
+        plt.plot(cr_pts, cr_vals, 'o') # The extrema
+
+        """
+        # Remove NaNs
+        I = np.isnan(y)
+        if sum(I) > 0:
+            print('find_maxima(): There were NaNs in timeseries')
+            x = x[np.logical_not(I)]
+            y = y[np.logical_not(I)]
+
+        # Sort over time. Monotonic increasing
+        y = y.sortby(x)
+        x = x.sortby(x)
+
+        # Convert x to float64 (assuming y is/similar to np.float64)
+        if type(x.values[0]) == np.datetime64: # convert to decimal sec since 1970
+            x_float = ((x.values - np.datetime64('1970-01-01T00:00:00'))
+                           / np.timedelta64(1, 's')).astype('float64')
+            #x_float = x.values.astype('float64')
+            y_float = y.values.astype('float64')
+            flag_dt64 = True
+        else:
+            x_float = x.values.astype('float64')
+            y_float = y.values.astype('float64')
+            flag_dt64 = False
+
+        if type(y.values[0]) != np.float64:
+            print('find_maxima(): type(y)=',type(y))
+            print('I was expecting a np.float64')
+
+        # Do the interpolation
+        f = scipy.interpolate.InterpolatedUnivariateSpline(x_float, y, k=3)
+        # Find the extrema as roots
+        extr_x_vals = quadratic_spline_roots(f.derivative()) # x values of extrema
+        # Find the maxima roots
+        extr_x_vals = np.hstack( [x_float[0], extr_x_vals, x_float[-1]]) # add buffer points to ensure extrema are within
+        ind = scipy.signal.argrelmax( f(extr_x_vals) )[0] # index that gives max(f) over extrema x locations
+        max_vals = f(extr_x_vals[ind])
+
+        # Convert back to datetime64 if appropriate
+        y_out = max_vals
+        if flag_dt64:
+            N = len(extr_x_vals[ind])
+            x_out = [np.datetime64('1970-01-01T00:00:00')
+             + np.timedelta64(int(extr_x_vals[ind[i]]), 's') for i in range(N)]
+        else:
+            x_out = extr_x_vals[ind]
+
+        # restore xarray structure
+        new_x = xr.DataArray( x_out, coords=[x_out], dims=x.dims)
+        new_x.name = x.name
+        new_y = xr.DataArray( y_out, coords=[x_out], dims=y.dims)
+        new_y.name = y.name
+
+        return new_x, new_y
+
+
 def doodson_x0_filter(elevation, ax=0):
-    ''' 
-    The Doodson X0 filter is a simple filter designed to damp out the main 
-    tidal frequencies. It takes hourly values, 19 values either side of the 
+    '''
+    The Doodson X0 filter is a simple filter designed to damp out the main
+    tidal frequencies. It takes hourly values, 19 values either side of the
     central one and applies a weighted average using:
               (1010010110201102112 0 2112011020110100101)/30.
     ( http://www.ntslf.org/files/acclaimdata/gloup/doodson_X0.html )
-    
+
     In "Data Analaysis and Methods in Oceanography":
-    
+
     "The cosine-Lanczos filter, the transform filter, and the
     Butterworth filter are often preferred to the Godin filter,
     to earlier Doodson filter, because of their superior ability
     to remove tidal period variability from oceanic signals."
-    
+
     This routine can be used for any dimension input array.
-    
+
     Parameters
     ----------
         elevation (ndarray) : Array of hourly elevation values.
-        axis (int) : Time axis of input array. This axis must have >= 39 
+        axis (int) : Time axis of input array. This axis must have >= 39
         elements
-       
+
     Returns
     -------
         Filtered array of same rank as elevation.
-    ''' 
+    '''
     if elevation.shape[ax] < 39:
         print('Doodson_XO: Ensure time axis has >=39 elements. Returning.')
         return
     # Define DOODSON XO weights
-    kern = np.array([1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 2, 0, 1, 1, 0, 2, 1, 1, 2, 
+    kern = np.array([1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 2, 0, 1, 1, 0, 2, 1, 1, 2,
                      0,
                      2, 1, 1, 2, 0, 1, 1, 0, 2, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1])
     kern = kern/30
 
     # Convolve input array with weights along the specified axis.
-    filtered = np.apply_along_axis(lambda m: np.convolve(m, kern, mode=1), 
+    filtered = np.apply_along_axis(lambda m: np.convolve(m, kern, mode=1),
                                    axis=ax, arr=elevation)
 
     # Pad out boundary areas with NaNs for given (arbitrary) axis.
@@ -76,63 +184,3 @@ def doodson_x0_filter(elevation, ax=0):
     filtered[-19:] = np.nan
     filtered = filtered.swapaxes(0,ax)
     return filtered
-
-
-def normal_distribution(mu: float=0, sigma: float=1, 
-                        x: np.ndarray=None, n_pts: int=1000):
-    """Generates a discrete normal distribution.
-
-    Keyword arguments:
-    x     -- Arbitrary array of x-values
-    mu    -- Distribution mean
-    sigma -- Distribution standard deviation
-    
-    return: Array of len(x) containing the normal values calculated from
-            the elements of x.
-    """
-    if x is None:
-        x = np.linspace( mu-5*sigma, mu+5*sigma, n_pts)
-    debug(f"Generating normal distribution for {get_slug(x)}")
-    term1 = sigma*np.sqrt( 2*np.pi )
-    term1 = 1/term1
-    exponent = -0.5*((x-mu)/sigma)**2
-    return term1*np.exp( exponent )
-
-def cumulative_distribution(mu: float=0, sigma: float=1, 
-                            x: np.ndarray=None, cdf_func: str='gaussian'):
-    """Integrates under a discrete PDF to obtain an estimated CDF.
-
-    Keyword arguments:
-    x   -- Arbitrary array of x-values
-    pdf -- PDF corresponding to values in x. E.g. as generated using
-           normal_distribution.
-    
-    return: Array of len(x) containing the discrete cumulative values 
-            estimated using the integral under the provided PDF.
-    """
-    debug(f"Estimating CDF using {get_slug(x)}")
-    if cdf_func=='gaussian': #If Gaussian, integrate under pdf
-        pdf = normal_distribution(mu=mu, sigma=sigma, x=x)
-        cdf = [np.trapz(pdf[:ii],x[:ii]) for ii in range(0,len(x))]
-    else: 
-        raise NotImplementedError
-    return np.array(cdf)
-
-def empirical_distribution(x, sample):
-    """Estimates a CDF empirically.
-
-    Keyword arguments:
-    x      -- Array of x-values over which to generate distribution
-    sample -- Sample to use to generate distribution
-    
-    return: Array of len(x) containing corresponding EDF values
-    """
-    debug(f"Estimating empirical distribution with {get_slug(x)}")
-    sample = np.array(sample)
-    sample = sample[~np.isnan(sample)]
-    sample = np.sort(sample)
-    edf = np.zeros(len(x))
-    n_sample = len(sample)
-    for ss in sample:
-        edf[x>ss] = edf[x>ss] + 1/n_sample
-    return xr.DataArray(edf)
\ No newline at end of file
diff --git a/example_scripts/tidegauge_tutorial.py b/example_scripts/tidegauge_tutorial.py
index 2f8b503d..1265b383 100755
--- a/example_scripts/tidegauge_tutorial.py
+++ b/example_scripts/tidegauge_tutorial.py
@@ -8,6 +8,7 @@
 import coast
 import datetime
 import numpy as np
+import matplotlib.pyplot as plt
 
 # And by defining some file paths
 fn_nemo_dat  = './example_files/COAsT_example_NEMO_data.nc'
@@ -118,26 +119,6 @@
 fig, ax = TIDEGAUGE.plot_on_map_multiple(tidegauge_list, color_var_str='rmse')
 
 
-#%%  Additionally, alternative data streams can be read in and similarly
-# processed. For example the BODC processed data from the UK Tidegauge network.
-# Data name: UK Tide Gauge Network, processed data.
-# Source: https://www.bodc.ac.uk/
-
-# Load and plot BODC processed data
-fn_bodc = 'example_files/LIV2010.txt'
-
-# Set the start and end dates
-date_start = np.datetime64('2020-10-12 23:59')
-date_end = np.datetime64('2020-10-14 00:01')
-
-# Initiate a TIDEGAUGE object, if a filename is passed it assumes it is a GESLA
-# type object
-tg = coast.TIDEGAUGE()
-# specify the data read as a High Low Water dataset
-tg.dataset = tg.read_bodc_to_xarray(fn_bodc, date_start, date_end)
-tg.plot_timeseries()
-
-
 #%% Alternatively load in data obtained using the Environment Agency (England)
 #  API. These are only accessible for the last 28 days. This does not require
 # an API key.
@@ -157,3 +138,47 @@
 # (the default) station.
 eg.dataset = eg.read_EA_API_to_xarray(ndays=1, stationId="E70124")
 eg.plot_timeseries()
+
+#%%  Additionally, alternative data streams can be read in and similarly
+# processed. For example the BODC processed data from the UK Tidegauge network.
+# Data name: UK Tide Gauge Network, processed data.
+# Source: https://www.bodc.ac.uk/
+
+#%% Tide table extrema can be extracted from the timeseries
+#        Finds high and low water for a given variable.
+#        Returns in a new TIDEGAUGE object with similar data format to
+#        a TIDETABLE.
+# There are two methods for extracting extrema. A neightbour comparison method, 
+# (method="comp") and a cubic spline fitting method (method="cubic")
+
+# Load and plot BODC processed data
+fn_bodc = 'example_files/tide_gauges/LIV2010.txt'
+
+# Set the start and end dates
+date_start = np.datetime64('2020-10-13 20:00')
+date_end = np.datetime64('2020-10-13 21:00')
+
+# Initiate a TIDEGAUGE object, if a filename is passed it assumes it is a GESLA
+# type object
+del tg
+tg = coast.TIDEGAUGE()
+tg.dataset = tg.read_bodc_to_xarray(fn_bodc, date_start, date_end)
+
+# Use comparison of neighbourhood method (method="comp" is assumed)
+extrema_comp = tg.find_high_and_low_water('sea_level', method="comp")
+
+# Use cubic spline fitting method
+extrema_cubc = tg.find_high_and_low_water('sea_level', method="cubic")
+
+# Plot to show the difference between maxima find methods for high tide example.
+plt.figure()
+plt.plot(tg.dataset.time, tg.dataset.sea_level, 'k')
+plt.scatter(extrema_comp.dataset.time_highs.values, extrema_comp.dataset.sea_level_highs, marker='o', c='g')
+plt.scatter(extrema_cubc.dataset.time_highs.values, extrema_cubc.dataset.sea_level_highs, marker='+', c='g')
+plt.xlim([date_start, date_end])
+plt.ylim([7.75, 8.0])
+plt.legend(['Time Series','Maxima by comparison','Maxima by cubic spline'])
+plt.title('Tide Gauge Optima at Gladstone')
+plt.show()
+
+
diff --git a/requirements.txt b/requirements.txt
index 39ed8ee8..1c6347c6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 .
-PyYAML==5.3.1
+PyYAML==5.4
 oyaml==1.0
 pytest==6.0.2
 numpy>=1.16
diff --git a/setup.py b/setup.py
index 96e1d4a1..97c3fa89 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
 
 PACKAGE = SimpleNamespace(**{
     "name": "COAsT",
-    "version": "0.4.6",
+    "version": "0.5.0",
     "description": "This is the Coast Ocean Assessment Tool",
     "url": "https://www.bodc.ac.uk",
     "download_url": "https://github.com/British-Oceanographic-Data-Centre/COAsT/",


Test patch:
diff --git a/unit_testing/unit_test.py b/unit_testing/unit_test.py
index aa6b0c42..68fd4a1f 100755
--- a/unit_testing/unit_test.py
+++ b/unit_testing/unit_test.py
@@ -1,38 +1,28 @@
 ##
 """
 Script to do unit testing
-
 Written as procedural code that plods through the code snippets and tests the
 outputs or expected metadata.
-
 ***SECTIONS***
 This script is separated into Sections and Subsections, for which there are two
 counters to keep track: sec and subsec respectively.  At the beginning of each
 section, the sec counter should be incremented by 1 and the subsec counter
 should be reset to 96 (code for one below 'a'). At the beginning of each
 subsection, subsec should be incremented by one.
-
 ***OTHER FILES***
 There are two accompaniment files to this unit testing script:
     - unit_test_contents: A list of sections and subsections.
     - unit_test_guidelines: Further guidelines to creating unit tests.
-
 Run:
 ipython: cd COAsT; run unit_testing/unit_test.py  # I.e. from the git repo.
-
 Unit template:
-
-
 #-----------------------------------------------------------------------------#
 # ( ## ) Subsection title                                                     #
 #                                                                             #
-
 subsec = subsec+1
 # <Introduction>
-
 try:
     # Do a thing
-
     #TEST: <description here>
     check1 = #<Boolean>
     check2 = #<Boolean>
@@ -40,10 +30,8 @@
         print(str(sec) + chr(subsec) + " OK - ")
     else:
         print(str(sec) + chr(subsec) + " X - ")
-
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
-
 """
 
 import coast
@@ -55,6 +43,8 @@
 import os.path as path
 import logging
 import coast.general_utils as general_utils
+import coast.plot_util as plot_util
+import coast.stats_util as stats_util
 from socket import gethostname# to get hostname
 import traceback
 import xarray.ufuncs as uf
@@ -284,7 +274,7 @@
 # NEMO obejct and dataset.
 
 try:
-    harmonics = coast.NEMO(dn_files + fn_nemo_harmonics, 
+    harmonics = coast.NEMO(dn_files + fn_nemo_harmonics,
                            dn_files + fn_nemo_harmonics_dom)
     constituents = ['K1','M2','S2','K2']
     harmonics_combined = harmonics.harmonics_combine(constituents)
@@ -299,7 +289,7 @@
 
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
-    
+
 #-----------------------------------------------------------------------------#
 #%% ( 1j ) Convert harmonics to a/g and back                                  #
 #                                                                             #
@@ -324,6 +314,27 @@
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
 
+#-----------------------------------------------------------------------------#
+#%% ( 1k ) Compute e3 from SSH field                                      #
+#
+subsec = subsec+1
+try:
+    nemo_t = coast.NEMO( fn_data=dn_files+fn_nemo_grid_t_dat,
+                        fn_domain=dn_files+fn_nemo_dom, grid_ref='t-grid' )
+    
+    e3t,e3u,e3v,e3f,e3w = coast.NEMO.get_e3_from_ssh(nemo_t,True,True,True,True,True)
+    cksum = np.array([e3t.sum(),e3u.sum(),e3v.sum(),
+                      e3f.sum(),e3w.sum()])
+    # these references are based on the example file's ssh field
+    reference = np.array([8.337016e+08, 8.333972e+08, 8.344886e+08,
+                          8.330722e+08, 8.265948e+08])
+    if np.allclose(cksum, reference):
+        print(str(sec) + chr(subsec) + " OK - computed e3[t,u,v,f,w] as expected")
+    else:
+        print(str(sec) + chr(subsec) + " X - computed e3[t,u,v,f,w] not as expected")        
+except:
+    print(str(sec) + chr(subsec) + ' FAILED.\n' + traceback.format_exc())
+
 '''
 #################################################
 ## ( 2 ) Test general utility methods in COAsT ##
@@ -1253,7 +1264,7 @@
     print(str(sec) + chr(subsec) +' FAILED.')
 
 #-----------------------------------------------------------------------------#
-#%% ( 7n ) TIDEGAUGE method for finding extrema and troughs                   #
+#%% ( 7n ) TIDEGAUGE method for finding extrema and troughs, compare neighbours#
 #                                                                             #
 subsec = subsec+1
 
@@ -1264,31 +1275,75 @@
     lowestoft2 = coast.TIDEGAUGE(fn_tidegauge, date_start = date0,
                                 date_end = date1)
 
-    extrema = lowestoft2.find_high_and_low_water('sea_level', distance=40)
-
+    # Use comparison of neighbourhood method (method="comp" is assumed)
+    extrema_comp = lowestoft2.find_high_and_low_water('sea_level', distance=40)
     # Check actual maximum/minimum is in output dataset
-    check1 = np.nanmax(lowestoft2.dataset.sea_level) in extrema.dataset.sea_level_highs
-    check2 = np.nanmin(lowestoft2.dataset.sea_level) in extrema.dataset.sea_level_lows
+    check1 = np.nanmax(lowestoft2.dataset.sea_level) in extrema_comp.dataset.sea_level_highs
+    check2 = np.nanmin(lowestoft2.dataset.sea_level) in extrema_comp.dataset.sea_level_lows
     # Check new time dimensions have correct length (hardcoded here)
-    check3 = len(extrema.dataset.time_highs) == 19
-    check4 = len(extrema.dataset.time_lows) == 18
+    check3 = len(extrema_comp.dataset.time_highs) == 19
+    check4 = len(extrema_comp.dataset.time_lows) == 18
 
     # Attempt a plot
     f = plt.figure()
     plt.plot(lowestoft2.dataset.time, lowestoft2.dataset.sea_level)
-    plt.scatter(extrema.dataset.time_highs, extrema.dataset.sea_level_highs, c='g')
-    plt.scatter(extrema.dataset.time_lows, extrema.dataset.sea_level_lows, c='r')
+    plt.scatter(extrema_comp.dataset.time_highs.values, extrema_comp.dataset.sea_level_highs, marker='o', c='g')
+    plt.scatter(extrema_comp.dataset.time_lows.values,  extrema_comp.dataset.sea_level_lows, marker='o', c='r')
+
     plt.legend(['Time Series','Maxima','Minima'])
     plt.title('Tide Gauge Optima at Lowestoft')
     f.savefig(dn_fig + 'tidegauge_optima.png')
 
     if check1 and check2 and check3 and check4:
-        print(str(sec) + chr(subsec) + " OK - Tidegauge extrema found")
+        print(str(sec) + chr(subsec) + " OK - Tidegauge local extrema found")
+    else:
+        print(str(sec) + chr(subsec) + " X - Tidegauge local extrema")
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+
+#-----------------------------------------------------------------------------#
+#%% ( 7o ) TIDEGAUGE method for finding extrema and troughs, fit cubic spline #
+#                                                                             #
+subsec = subsec+1
+
+# Load and process BODC processed data
+try:
+    # Set the start and end dates
+    date_start = np.datetime64('2020-10-12 23:59')
+    date_end = np.datetime64('2020-10-14 00:01')
+
+    # Initiate a TIDEGAUGE object, if a filename is passed it assumes it is a GESLA
+    # type object
+    tg = coast.TIDEGAUGE()
+    # specify the data read as a High Low Water dataset
+    tg.dataset = tg.read_bodc_to_xarray(fn_tidegauge2, date_start, date_end)
+
+    # Use cubic spline fitting method
+    extrema_cubc = tg.find_high_and_low_water('sea_level', method="cubic")
+
+    # Check actual maximum/minimum is in output dataset
+    check1 = np.isclose( extrema_cubc.dataset.sea_level_highs,[7.77432795, 7.91244559])
+    check2 = np.isclose( extrema_cubc.dataset.sea_level_lows,[2.63479458, 2.54599355])
+                                
+    # Attempt a plot
+    f = plt.figure()
+    plt.plot(tg.dataset.time, tg.dataset.sea_level)
+    plt.scatter(extrema_cubc.dataset.time_highs.values, extrema_cubc.dataset.sea_level_highs, marker='o', c='g')
+    plt.scatter(extrema_cubc.dataset.time_lows.values,  extrema_cubc.dataset.sea_level_lows, marker='o', c='r')
+
+    plt.legend(['Time Series','Maxima','Minima'])
+    plt.title('Tide Gauge Optima at Gladstone, fitted cubic spline')
+    f.savefig(dn_fig + 'tidegauge_optima.png')
+
+    if check1.all() and check2.all():
+        print(str(sec) + chr(subsec) + " OK - Tidegauge cubic extrema found")
     else:
-        print(str(sec) + chr(subsec) + " X - Tidegauge extrema")
+        print(str(sec) + chr(subsec) + " X - Tidegauge cubic extrema")
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
 
+
 '''
 ###############################################################################
 ## ( 8 ) Isobath Contour Methods                                            ##
@@ -1437,140 +1492,404 @@
 
 except:
     print(str(sec) + chr(subsec) + ' FAILED.\n' + traceback.format_exc())
+    
+'''
+#################################################
+## ( 10 ) PROFILE Methods                     ##
+#################################################
+'''
+sec = sec+1
+subsec = 96
+
+#-----------------------------------------------------------------------------#
+# ( 10a ) Load EN4 data                                                       #
+#                                                                             #
+
+subsec = subsec+1
+# Create PROFILE object and read EN4 example data file
+
+try:
+    profiles = coast.PROFILE()
+    profiles.read_EN4(fn_EN4)
+
+    #TEST: Check some data
+    check1 = profiles.dataset.dims['z_dim'] == 400
+    check2 = profiles.dataset.longitude[11].values == 9.89777
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - EN4 Data read, PROFILE created")
+    else:
+        print(str(sec) + chr(subsec) + " X - Problem with EN4 reading")
 
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+
+#-----------------------------------------------------------------------------#
+# ( 10b ) Plot locations on map                                               #
+#                                                                             #
+
+subsec = subsec+1
+# Plot profile locations on a map
+
+try:
+    f,a = profiles.plot_map()
+    f.savefig(dn_fig + 'profiles_map.png')
+    print(str(sec) + chr(subsec) + " OK - Profiles map plot saved")
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+#-----------------------------------------------------------------------------#
+# ( 10c ) Plot ts diagram                                                     #
+#                                                                             #
+
+subsec = subsec+1
+# Plot ts diagram
+
+try:
+    f,a = profiles.plot_ts_diagram(10)
+    f.savefig(dn_fig + 'profile_ts_diagram.png')
+    print(str(sec) + chr(subsec) + " OK - Profiles ts diagram plot saved")
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+#-----------------------------------------------------------------------------#
+# ( 10d ) Plot temperature profile                                            #
+#                                                                             #
+
+subsec = subsec+1
+# Plot ts diagram
+
+try:
+    f,a = profiles.plot_profile(var='potential_temperature',profile_indices=[10])
+    f.savefig(dn_fig + 'profile_temperature_diagram.png')
+    print(str(sec) + chr(subsec) + " OK - Profiles temperature plot saved")
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
 #%%
 '''
-###############################################################################
-## ( 10 ) Example script testing                                              ##
-###############################################################################
+#################################################
+## ( 11 ) PLOTTING UTILITY Methods             ##
+#################################################
 '''
 sec = sec+1
 subsec = 96
 
-print(str(sec) + ". Example script testing")
-print("++++++++++++++++++++++++")
-#
 #-----------------------------------------------------------------------------#
-#%% ( 10a ) Example script testing                                               #
+# ( 11a ) Scatter with fit                                                     #
 #                                                                             #
+
 subsec = subsec+1
-# Test machine name (to check for file access) in order to test additional scripts.
-example_script_flag = True if 'livljobs' in gethostname().lower() else False
+# <Introduction>
 
 try:
-    # Do a thing
-    from example_scripts import altimetry_tutorial # This runs on example_files
-    from example_scripts import tidegauge_tutorial # This runs on example_files
-    from example_scripts import tidetable_tutorial # This runs on example_files
-    from example_scripts import export_to_netcdf_tutorial # This runs on example_files
+    # Plot an idealised dataset
+    x = np.arange(0,50)
+    y = np.arange(0,50)/1.5
+    f,a = plot_util.scatter_with_fit(x,y)
+    a.set_title('Test: Scatter_with_fit()')
+
+    f.savefig(os.path.join(dn_fig, "scatter_with_fit_test.png"))
+    plt.close()
+    
+    print(str(sec) + chr(subsec) +' OK. Scatter_with_fit()')
 
-    print(str(sec) + chr(subsec) + " OK - tutorials on example_files data")
-    subsec = subsec+1
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+    
+#-----------------------------------------------------------------------------#
+# ( 11b ) Geo axes                                                             #
+#                                                                             #
 
-    if example_script_flag:
-        from example_scripts import AMM15_example_plot
-        print(str(sec) + chr(subsec) + " OK - tutorial on AMM15 data")
-        subsec = subsec+1
-        from example_scripts import ANChor_plots_of_NSea_wvel
-        print(str(sec) + chr(subsec) + " OK - tutorial on AMM60 data")
-        subsec = subsec+1
-        from example_scripts import BLZ_example_plot
-        print(str(sec) + chr(subsec) + " OK - tutorial on Belize data")
-        subsec = subsec+1
-        from example_scripts import SEAsia_R12_example_plot
-        print(str(sec) + chr(subsec) + " OK - tutorial on SEAsia data")
-        subsec = subsec+1
-        from example_scripts import WCSSP_India_example_plot
-        print(str(sec) + chr(subsec) + " OK - tutorial on WCSSP-India data")
-        subsec = subsec+1
-        from example_scripts import internal_tide_pycnocline_diagnostics
-        print(str(sec) + chr(subsec) + " OK - tutorial on internal tides")
+subsec = subsec+1
+
+try:
+    # Plot two scatters on a map
+    lonbounds = [-20,20]
+    latbounds = [30,60]
+    f,a = plot_util.create_geo_axes(lonbounds, latbounds)
+    a.set_title('Test: create_geo_axes()')
+    a.scatter([0,-10],[50,50])
+
+    f.savefig(os.path.join(dn_fig, "create_geo_axes_test.png"))
+    plt.close()
+    
+    print(str(sec) + chr(subsec) +' OK. create_geo_axes()')
+
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+    
+#-----------------------------------------------------------------------------#
+# ( 11c ) Determine colorbar extension                                         #
+#                                                                             #
+
+subsec = subsec+1
+
+try:
+    # Test on some pretend dataset
+    pretend_data = np.arange(0,50)
+    test1 = plot_util.determine_colorbar_extension(pretend_data, -50, 100)
+    test2 = plot_util.determine_colorbar_extension(pretend_data, 1, 100)
+    test3 = plot_util.determine_colorbar_extension(pretend_data, -50, 48)
+    test4 = plot_util.determine_colorbar_extension(pretend_data, 1, 48)
+
+    #TEST: <description here>
+    check1 = test1 == 'neither'
+    check2 = test2 == 'min'
+    check3 = test3 == 'max'
+    check4 = test4 == 'both'
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - determine_colorbar_extension()")
     else:
-        print("Don't forget to test on a LIVLJOBS machine")
+        print(str(sec) + chr(subsec) + " X - ")
+
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+#-----------------------------------------------------------------------------#
+# ( 11d ) Determine clim by st.dev                                             #
+#                                                                             #
+
+subsec = subsec+1
+# <Introduction>
+
+try:
+    # Data with mean = 50.51 and std = 32.2
+    # Determine clims to exclude outlier at 200
+    pretend_data = np.arange(0,100)
+    pretend_data[-1] = 200
+    clim = plot_util.determine_clim_by_standard_deviation(pretend_data, 
+                                                          n_std_dev=2)
 
     #TEST: <description here>
-    check1 = example_script_flag
-    if check1:
-        print(str(sec) + " OK - example_scripts ran on",gethostname())
+    check1 = clim[0] == -13.808889915793792
+    check2 = clim[1] == 114.82888991579378
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - determine_clim_by_std_dev")
     else:
-        print(str(sec) + " X - example_scripts failed on",gethostname())
+        print(str(sec) + chr(subsec) + " X - determine_clim_by_std_dev")
 
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
 
+#%%
 '''
 #################################################
-## ( 11 ) PROFILE Methods                     ##
+## ( 12 ) Stats Utility                        ##
 #################################################
 '''
+
 sec = sec+1
 subsec = 96
 
 #-----------------------------------------------------------------------------#
-# ( 11a ) Load EN4 data                                                       #
+# ( 12a ) find_maxima(). Test comparison and cublic spline methods            #
 #                                                                             #
 
 subsec = subsec+1
-# Create PROFILE object and read EN4 example data file
+
 
 try:
-    profiles = coast.PROFILE()
-    profiles.read_EN4(fn_EN4)
+    date0 = datetime.datetime(2007,1,15)
+    date1 = datetime.datetime(2007,1,16)
+    tg = coast.TIDEGAUGE(fn_tidegauge, date_start = date0, date_end = date1)
 
-    #TEST: Check some data
-    check1 = profiles.dataset.dims['N_LEVELS'] == 400
-    check2 = profiles.dataset.longitude[11].values == 9.89777
-    if check1 and check2:
-        print(str(sec) + chr(subsec) + " OK - EN4 Data read, PROFILE created")
+    tt,hh = stats_util.find_maxima(tg.dataset.time, tg.dataset.sea_level, method='comp')
+    check1 = np.isclose( (tt.values[0]- np.datetime64('2007-01-15T00:15:00'))/ np.timedelta64(1,'s'), 0)
+    check2 = np.isclose(hh.values[0], 1.027)
+    
+    tt,hh = stats_util.find_maxima(tg.dataset.time, tg.dataset.sea_level, method='cubic')
+    check3 = np.isclose( (tt[0] - np.datetime64('2007-01-15T00:07:49'))/ np.timedelta64(1,'s'), 0 )
+    check4 = np.isclose( hh[0], 1.0347638302097757 )
+
+
+    if check1 and check2 and check3 and check4:
+        print(str(sec) + chr(subsec) + " OK - find_maxima() worked for comparison and cublic spline methods")
     else:
-        print(str(sec) + chr(subsec) + " X - Problem with EN4 reading")
+        print(str(sec) + chr(subsec) + " X - Problem with stats_util.find_maxima()")
 
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
-    
+
+#%%
+'''
+#################################################
+## ( 13 ) MASK_MAKER                           ##
+#################################################
+'''
+
+sec = sec+1
+subsec = 96
+
+# Preparation: Create two arrays to put mask onto, one of zeros and one of ones
+# This allows us to test the additive feature.
+sci = coast.NEMO(dn_files + fn_nemo_dat, dn_files + fn_nemo_dom, grid_ref = 't-grid')
+mask00 = np.zeros((sci.dataset.dims['y_dim'], sci.dataset.dims['x_dim']))
+mask01 = np.ones((sci.dataset.dims['y_dim'], sci.dataset.dims['x_dim']))
 
 #-----------------------------------------------------------------------------#
-# ( 11b ) Plot locations on map                                               #
+# ( 13a ) Create mask by indices                                              #
 #                                                                             #
 
 subsec = subsec+1
-# Plot profile locations on a map
+# Plot ts diagram
 
 try:
-    f,a = profiles.plot_map()
-    f.savefig(dn_fig + 'profiles_map.png')
-    print(str(sec) + chr(subsec) + " OK - Profiles map plot saved")
+    mm = coast.MASK_MAKER()
+    # Draw and fill a square
+    vertices_r = [50, 150, 150, 50]
+    vertices_c = [50, 50, 150, 150]
+    filled0 = mm.fill_polygon_by_index(mask00, vertices_r, vertices_c)
+    filled1 = mm.fill_polygon_by_index(mask01, vertices_r, vertices_c, additive=True)
+
+    #TEST: Check some data
+    check1 = filled0[49,49] == 0 and filled0[51,51] == 1
+    check2 = filled1[49,49] == 1 and filled1[51,51] == 2
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - MASKS created by index")
+    else:
+        print(str(sec) + chr(subsec) + " X - Problem mask creation by index")
+
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
     
 #-----------------------------------------------------------------------------#
-# ( 11c ) Plot ts diagram                                                     #
+# ( 13b ) Create mask by lonlat                                               #
 #                                                                             #
 
 subsec = subsec+1
 # Plot ts diagram
 
 try:
-    f,a = profiles.plot_ts_diagram(10)
-    f.savefig(dn_fig + 'profile_ts_diagram.png')
-    print(str(sec) + chr(subsec) + " OK - Profiles ts diagram plot saved")
+    mm = coast.MASK_MAKER()
+    # Draw and fill a square
+    vertices_lon = [-5, -5, 5, 5]
+    vertices_lat = [40, 60, 60, 40]
+    filled0 = mm.fill_polygon_by_lonlat(mask00, sci.dataset.longitude, 
+                                        sci.dataset.latitude, vertices_lon, 
+                                        vertices_lat)
+    filled1 = mm.fill_polygon_by_lonlat(mask01, sci.dataset.longitude, 
+                                        sci.dataset.latitude, vertices_lon, 
+                                        vertices_lat, additive=True)
+
+    #TEST: Check some data
+    check1 = filled0[50,50] == 0 and filled0[50,150] == 1
+    check2 = filled1[50,50] == 1 and filled1[50,150] == 2
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - MASKS created by lonlat")
+    else:
+        print(str(sec) + chr(subsec) + " X - Problem mask creation by lonlat")
+
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
     
+#%%
+'''
+#################################################
+## ( 14 ) CLIMATOLOGY Methods                  ##
+#################################################
+'''
+sec = sec+1
+subsec = 96
+
+sci = coast.NEMO(dn_files + fn_nemo_dat, dn_files + fn_nemo_dom, grid_ref = 't-grid')
+ds = sci.dataset[['temperature','ssh']].isel(z_dim=0)
+
+
 #-----------------------------------------------------------------------------#
-# ( 11d ) Plot temperature profile                                            #
+# ( 14a ) Monthly and Seasonal Climatology                                     #
 #                                                                             #
 
 subsec = subsec+1
-# Plot ts diagram
 
 try:
-    f,a = profiles.plot_profile(var='POTM_CORRECTED',profile_indices=[10])
-    f.savefig(dn_fig + 'profile_temperature_diagram.png')
-    print(str(sec) + chr(subsec) + " OK - Profiles temperature plot saved")
+    
+    clim = coast.CLIMATOLOGY()
+    fn_out = os.path.join(dn_files, 'test_climatology.nc')
+    monthly = clim.make_climatology(ds, 'month').load()
+    seasonal = clim.make_climatology(ds, 'season', fn_out=fn_out)
+    
+    # create dataset with missing values
+    ds2 = ds.copy(deep=True)
+    ds2['temperature'][::2, :100,:100] = np.nan
+    ds2['ssh'][::2, :100,:100] = np.nan
+    seaC = clim.make_climatology(ds2, 'season', missing_values=True)
+    seaX = ds2.groupby("time.season").mean('t_dim')
+    # throws error is not close
+    xr.testing.assert_allclose(seaC,seaX)
+       
+    mn = mn = np.nanmean(ds.temperature, axis=0)
+    check1 = np.nanmax(np.abs(mn - monthly.temperature[0])) < 1e-6
+    check2 = os.path.isfile(fn_out)
+    if check1 and check2:
+        print(str(sec) + chr(subsec) + " OK - Monthly and seasonal climatology made and written to file")
+    else:
+        print(str(sec) + chr(subsec) + " X - Problem with monthly and seasonal climatology ")
+
+except AssertionError:
+    print(str(sec) + chr(subsec) + " X - Problem with computing climatology when dataset has missing values")
+except:
+    print(str(sec) + chr(subsec) +' FAILED.')
+
+#%%
+'''
+###############################################################################
+## ( N ) Example script testing                                              ##
+###############################################################################
+'''
+sec = 'N'
+subsec = 96
+
+print(str(sec) + ". Example script testing")
+print("++++++++++++++++++++++++")
+#
+#-----------------------------------------------------------------------------#
+#%% ( Na ) Example script testing                                               #
+#                                                                             #
+subsec = subsec+1
+# Test machine name (to check for file access) in order to test additional scripts.
+example_script_flag = True if 'livljobs' in gethostname().lower() else False
+
+try:
+    # Do a thing
+    from example_scripts import altimetry_tutorial # This runs on example_files
+    from example_scripts import tidegauge_tutorial # This runs on example_files
+    from example_scripts import tidetable_tutorial # This runs on example_files
+    from example_scripts import export_to_netcdf_tutorial # This runs on example_files
+
+    print(str(sec) + chr(subsec) + " OK - tutorials on example_files data")
+    subsec = subsec+1
+
+    if example_script_flag:
+        from example_scripts import AMM15_example_plot
+        print(str(sec) + chr(subsec) + " OK - tutorial on AMM15 data")
+        subsec = subsec+1
+        from example_scripts import ANChor_plots_of_NSea_wvel
+        print(str(sec) + chr(subsec) + " OK - tutorial on AMM60 data")
+        subsec = subsec+1
+        from example_scripts import BLZ_example_plot
+        print(str(sec) + chr(subsec) + " OK - tutorial on Belize data")
+        subsec = subsec+1
+        from example_scripts import SEAsia_R12_example_plot
+        print(str(sec) + chr(subsec) + " OK - tutorial on SEAsia data")
+        subsec = subsec+1
+        from example_scripts import WCSSP_India_example_plot
+        print(str(sec) + chr(subsec) + " OK - tutorial on WCSSP-India data")
+        subsec = subsec+1
+        from example_scripts import internal_tide_pycnocline_diagnostics
+        print(str(sec) + chr(subsec) + " OK - tutorial on internal tides")
+    else:
+        print("Don't forget to test on a LIVLJOBS machine")
+
+    #TEST: <description here>
+    check1 = example_script_flag
+    if check1:
+        print(str(sec) + " OK - example_scripts ran on",gethostname())
+    else:
+        print(str(sec) + " X - example_scripts failed on",gethostname())
+
 except:
     print(str(sec) + chr(subsec) +' FAILED.')
 
 #%% Close log file
 #################################################
-log_file.close()
+log_file.close()
\ No newline at end of file
diff --git a/unit_testing/unit_test_contents.txt b/unit_testing/unit_test_contents.txt
index 95acac21..3cfa1914 100755
--- a/unit_testing/unit_test_contents.txt
+++ b/unit_testing/unit_test_contents.txt
@@ -1,4 +1,3 @@
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 |    Unit Testing Contents                                                    |
 *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 
@@ -15,7 +14,8 @@
     h. Load and combine multiple files.
     i. Load and combine harmonics 
     j. Convert harmonics to a/g and back 
-    
+    k. Compute e3[t,u,v,f,w] from SSH 
+
 2. General Utility Methods in COAsT
     a. Copying a COAsT object
     b. COAsT __getitem__ returns variable
@@ -64,7 +64,8 @@
     k. Plotting multiple tidegauge locations
     l. Tidegauge time series plot
     m. TIDEGAUGE method for tabulated data
-    n. TIDEGAUGE method for finding peaks and troughs
+    n. TIDEGAUGE method for finding peaks and troughs, compare neighbours
+    o. TIDEGAUGE method for finding extrema and troughs, fit cubic spline
 
 8. Isobath Contour Methods
     a. Extract isbath contour between two points
@@ -76,12 +77,34 @@
 9. EOF methods
     a. Compute EOFs, projections and variance
     b. Compute HEOFs, projections and variance
+    
+10. Profile Methods
+    a. Load EN4 data
+    b. Plot locations on map
+    c. Plot TS diagram
+    d. Plot temperature profile
+
+11. Plotting Utility
+    a. scatter_with_fit()
+    b. create_geo_axes()
+    c. determine_colorbar_extension()
+    d. determine_clim_by_standard_deviation()
+    
+12. Stats Utility
+    a. find_maxima()
+    
+13. MASK_MASKER
+    a. Create mask by indices
+    b. Create mask by lonlat
+
+14. CLIMATOLOGY
+    a. Create monthly and seasonal climatology, write to file
 
-10. Example script testing
+N. Example script testing
     a. tutorials using example_files (altimetry and tidegauges)
     b. tutorial on AMM15 data
     c. tutorial on AMM60 data
     d. tutorial on Belize data
     e. tutorial on SEAsia data
     f. tutorial on WCSSP-India data
-    g. tutorial on internal tides
+    g. tutorial on internal tides
\ No newline at end of file