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

Issue:
micro edit
micro edit (white space) to trigger git actions


Discussion / hints:

Fix patch:
diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md
index 24a13e46..597f6ed8 100644
--- a/.github/PULL_REQUEST_TEMPLATE.md
+++ b/.github/PULL_REQUEST_TEMPLATE.md
@@ -6,8 +6,8 @@ Please check if your PR fulfills the following requirements:
 - [ ] Tests for the changes have been added (for bug fixes / features)
 - [ ] Docs have been reviewed and added / updated if needed (for bug fixes / features)
     - Some docs updates need to be made in the `COAsT-site` repo, in a separate PR. See [contributing to documentation ](https://british-oceanographic-data-centre.github.io/COAsT/docs/contributing-docs/) for details.
-- [ ] Build (`./build.sh`) was run locally and no errors reported
-- [ ] Lint (`pylint .`) has passed locally and any fixes were made for failures
+- [ ] Build (`./build.sh`) was run locally and no errors reported. NB not sure about this requirement: GitActions test this
+- [ ] Lint (`pylint .`) has passed locally and any fixes were made for failures. NB not sure about this requirement: GitActions test this with `black`
 
 
 ## Pull request type
@@ -42,7 +42,7 @@ Issue URL: #
 ### Checklist
 - [ ] Code compiles correctly
 - [ ] Created tests which fail without the change (if possible)
-- [ ] All tests passing
+- [ ] All tests passing (`pip install . && pytest unit_testing/unit_test.py -s`)
 - [ ] Extended the README / documentation, if necessary
 - [ ] Added myself / the copyright holder to the AUTHORS file
 
@@ -56,4 +56,4 @@ Issue URL: #
 
 ## Other information
 
-<!-- Any other information that is important to this PR such as screenshots of how the component looks before and after the change. -->
\ No newline at end of file
+<!-- Any other information that is important to this PR such as screenshots of how the component looks before and after the change. -->
diff --git a/CITATION.cff b/CITATION.cff
index 18b99bb7..1cd03df8 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -37,6 +37,9 @@ authors:
  - family-names: "Gorman"
    given-names: "Luke"
    orcid: "https://orcid.org/0000-0002-2414-4311"
-title: "British-Oceanographic-Data-Centre/COAsT: v3.1.2"
-version: v3.1.2
-date-released: 2022-10-20
+ - family-names: "de Mora"
+   given-names: "Lee"
+   orcid: "https://orcid.org/0000-0002-5080-3149"
+title: "British-Oceanographic-Data-Centre/COAsT: v3.2.0"
+version: v3.2.0
+date-released: 2022-11-23
diff --git a/coast/_utils/general_utils.py b/coast/_utils/general_utils.py
index 5a167903..9ba438e9 100644
--- a/coast/_utils/general_utils.py
+++ b/coast/_utils/general_utils.py
@@ -348,3 +348,23 @@ def day_of_week(date: np.datetime64 = None):
         return "Tue"
     elif val == 6:
         return "Wed"
+
+
+def nan_helper(y):
+    """Helper to handle indices and logical indices of NaNs.
+
+    Input:
+        - y, 1d numpy array, or xr.DataArray, with possible NaNs
+    Output:
+        - nans, logical indices of NaNs
+        - index, a function, with signature indices= index(logical_indices),
+          to convert logical indices of NaNs to 'equivalent' indices
+    Example:
+        >>> # linear interpolation of NaNs
+        >>> nans, x= nan_helper(y)
+        >>> y[nans]= np.interp(x(nans), x(~nans), y[~nans])
+    """
+    if isinstance(y, np.ndarray):
+        return np.isnan(y), lambda z: z.nonzero()[0]
+    else:
+        return np.isnan(y).values, lambda z: z.nonzero()[0]
diff --git a/coast/_utils/mask_maker.py b/coast/_utils/mask_maker.py
index 07a70616..0ac0f74a 100644
--- a/coast/_utils/mask_maker.py
+++ b/coast/_utils/mask_maker.py
@@ -6,17 +6,58 @@
 
 
 class MaskMaker:
+    """
+    MaskMasker is a class of methods to assist with making regional masks within COAsT.
+    Presently these masks are external to MaskMaker.
+    It constructs a gridded boolean numpy array for each region, which are stacked over a dim_mask dimension and
+    stored as an xarray object.
+
+    A typical workflow might be:
+
+        # Define vertices
+        vertices_lon = [-5, -5, 5, 5]
+        vertices_lat = [40, 60, 60, 40]
+
+        # input lat/lon as xr.DataArray or numpy arrays. Return gridded boolean mask np.array on target grid
+        filled = mm.make_region_from_vertices(
+            sci.dataset.longitude, sci.dataset.latitude, vertices_lon, vertices_lat)
+
+        # make xr.Dataset of masks from gridded mask array or list of mask arrays
+        gridded_mask = mm.make_mask_dataset(sci.dataset.longitude.values,
+                                         sci.dataset.latitude.values,
+                                         filled)
+        # quick plot
+        mm.quick_plot(gridded_mask)
+
+
+    TO DO:
+    * Sort out region naming to be consistently applied and associated with the masks E.g. defined regions, or user defined masks
+    * Create final mask as a xr.DataArray, not a xr.Dataset
+    """
+
     def __init__(self):
 
         return
 
     @staticmethod
-    def make_mask_dataset(longitude, latitude, mask_list):
+    def make_mask_dataset(longitude, latitude, mask_list, mask_names: list = None):
+        """
+        create xr.Dataset for mask with latitude and longitude coordinates. If mask_names are given
+        create a dim_mask coordinate of names
+
+        """
         if type(mask_list) is not list:
             mask_list = [mask_list]
         gridded_mask = xr.Dataset()
         gridded_mask["longitude"] = (["y_dim", "x_dim"], longitude)
         gridded_mask["latitude"] = (["y_dim", "x_dim"], latitude)
+
+        if mask_names is not None:
+            gridded_mask["region_names"] = (["dim_mask"], mask_names)
+        else:
+            gridded_mask["region_names"] = (["dim_mask"], range(len(mask_list)))
+
+        gridded_mask = gridded_mask.set_coords(["longitude", "latitude", "region_names"])
         n_masks = len(mask_list)
         nr, nc = mask_list[0].shape
         all_masks = np.zeros((n_masks, nr, nc))
@@ -76,7 +117,7 @@ def fill_polygon_by_lonlat(
 
         Returns
         -------
-        Filled 2D array
+        Filled 2D np.array
         """
         array_to_fill = np.array(array_to_fill)
         ind_2d = general_utils.nearest_indices_2d(longitude, latitude, vertices_lon, vertices_lat)
@@ -105,12 +146,12 @@ def region_def_nws_north_sea(cls, longitude, latitude, bath):
     @classmethod
     def region_def_nws_outer_shelf(cls, longitude, latitude, bath):
         """
-        Regional definition for the Outher Shelf (Northwest European Shelf)
+        Regional definition for the Outer 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]
+        vertices_lon = [-4.1, -9.5, -1, 3.171, 3.171, -3.76, -3.76, -12, -12, -12, -4]
+        vertices_lat = [50.7, 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) * (~np.isnan(bath))
         return mask
@@ -135,8 +176,8 @@ def region_def_nws_english_channel(cls, longitude, latitude, bath):
         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]
+        vertices_lon = [-3.99, -3.99, -3.5, 12, 9]
+        vertices_lat = [51, 48.8, 48, 48, 51]
         mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
         mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
         return mask
@@ -172,3 +213,52 @@ def region_def_kattegat(cls, longitude, latitude, bath):
         mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
         mask = mask * (bath < 200) * (bath > 0) * (~np.isnan(bath))
         return mask
+
+    @classmethod
+    def make_region_from_vertices(cls, longitude, latitude, vertices_lon: list, vertices_lat: list):
+        """
+        Construct mask on supplied longitude, latitude grid with input lists of lon and lat polygon vertices
+        :param longitude: np.array/xr.DataArray of longitudes on target grid
+        :param latitude: np.array/xr.DataArray of latitudes on target grid
+        :param vertices_lon: list of vertices for bounding polygon
+        :param vertices_lat: list of vertices for bounding polygon
+        :return: mask: np.array(boolean) on target grid. Ones are bound by polygon vertices
+        """
+        try:
+            longitude = longitude.values
+            latitude = longitude.values
+        except AttributeError:
+            pass
+
+        mask = cls.fill_polygon_by_lonlat(np.zeros(longitude.shape), longitude, latitude, vertices_lon, vertices_lat)
+        return mask
+
+    @classmethod
+    def quick_plot(cls, mask: xr.Dataset):
+        """
+        Plot a map of masks in the MaskMaker object
+        Add labels
+        """
+        import matplotlib.pyplot as plt
+
+        n_mask = mask.dims["dim_mask"]
+        offset = 10  # nonzero offset to make scaled-boolean-masks [0, >offset]
+        for j in range(0, n_mask, 1):
+            tt = (j + offset) * mask["mask"].isel(dim_mask=j).squeeze()
+            ff = tt.where(tt > 0).plot(
+                x="longitude", y="latitude", levels=range(offset, n_mask + offset + 1, 1), add_colorbar=False
+            )
+
+        cbar = plt.colorbar(ff)
+        cbar.ax.get_yaxis().set_ticks([])
+        for j in range(0, n_mask, 1):
+            cbar.ax.text(
+                1 + 0.5,
+                offset + (j + 0.5),
+                mask["region_names"].isel(dim_mask=j).values,
+                ha="left",
+                va="center",
+                color="red",
+            )
+        cbar.ax.get_yaxis().labelpad = 15
+        plt.title(None)
diff --git a/coast/data/coast.py b/coast/data/coast.py
index c0966701..a9df20a3 100644
--- a/coast/data/coast.py
+++ b/coast/data/coast.py
@@ -85,7 +85,8 @@ def load_single(self, file: str, chunks: Dict = None):
             chunks (Dict): Chunks to use in Dask [default None].
         """
         info(f"Loading a single file ({file} for {get_slug(self)}")
-        self.dataset = xr.open_dataset(file, chunks=chunks)
+        with xr.open_dataset(file, chunks=chunks) as xrfile:
+            self.dataset = xrfile
 
     def load_multiple(self, directory_to_files: str, chunks: Dict = None):
         """Loads multiple files from directory into dataset variable.
@@ -95,7 +96,8 @@ def load_multiple(self, directory_to_files: str, chunks: Dict = None):
             chunks (Dict): Chunks to use in Dask [default None].
         """
         info(f"Loading a directory ({directory_to_files}) for {get_slug(self)}")
-        self.dataset = xr.open_mfdataset(directory_to_files, chunks=chunks, parallel=True, combine="by_coords")
+        with xr.open_mfdataset(directory_to_files, chunks=chunks, parallel=True, combine="by_coords") as files:
+            self.dataset = files
 
     def load_dataset(self, dataset: xr.Dataset):
         """Loads a dataset.
diff --git a/coast/data/tidegauge.py b/coast/data/tidegauge.py
index 14846f2c..fa54e2e9 100644
--- a/coast/data/tidegauge.py
+++ b/coast/data/tidegauge.py
@@ -1136,13 +1136,16 @@ def obs_operator(self, gridded, time_interp="nearest"):
         )
 
         # Interpolate model onto obs times
-        print("Interpolating in time...", flush=True)
-        extracted = extracted.rename({"time": "t_dim"})
-        extracted = extracted.interp(t_dim=ds.time.values, method=time_interp)
+        if "t_dim" in gridded.dims:
+            print("Interpolating in time...", flush=True)
+            extracted = extracted.rename_vars(
+                {"time": "t_dim"}
+            )  # make variable name match dimension name for interpolation
+            extracted = extracted.interp(t_dim=ds.time.values, method=time_interp)
+            extracted = extracted.rename_vars({"t_dim": "time"})  # restore variable name
 
         # Put interp_dist into dataset
-        extracted["interp_dist"] = interp_dist
-        extracted = extracted.rename_vars({"t_dim": "time"})
+        extracted["interp_dist"] = (("id_dim"), interp_dist)
 
         tg_out = Tidegauge()
         tg_out.dataset = extracted
diff --git a/coast/diagnostics/profile_analysis.py b/coast/diagnostics/profile_analysis.py
index 8ad50522..c815e1bc 100644
--- a/coast/diagnostics/profile_analysis.py
+++ b/coast/diagnostics/profile_analysis.py
@@ -210,9 +210,10 @@ def mask_means(cls, profile, mask_indices):
         dataset = profile.dataset
         # Get the number of masks provided
         n_masks = mask_indices.dims["dim_mask"]
+        mask_names = []
 
-        # Loop over maskal arrays. Assign mean to mask and seasonal means
-        debug(f"Calculating maskal averages..")
+        # Loop over masked arrays. Assign mean to mask and seasonal means
+        debug(f"Calculating masked averages..")
 
         # Loop over masks and index the data
         for mm in range(0, n_masks):
@@ -222,6 +223,8 @@ def mask_means(cls, profile, mask_indices):
                 if mm == 0:
                     ds_average = xr.Dataset()
                 continue
+            else:
+                mask_names.append(mask_indices["region_names"][mm].values)
 
             # Get actual profile data for this mask
             mask_data = dataset.isel(id_dim=mask_ind)
@@ -246,7 +249,12 @@ def mask_means(cls, profile, mask_indices):
             else:
                 ds_average = xr.concat((ds_average, ds_average_tmp), dim="dim_mask")
 
-        return ds_average
+        if mask_names != []:
+            ds_average["region_names"] = (["dim_mask"], mask_names)
+        else:
+            ds_average["region_names"] = (["dim_mask"], range(ds_average.dims["dim_mask"]))
+
+        return ds_average.set_coords(["region_names"])
 
     @classmethod
     def difference(cls, profile1, profile2, absolute_diff=True, square_diff=True):
diff --git a/coast/diagnostics/tidegauge_analysis.py b/coast/diagnostics/tidegauge_analysis.py
index 2367bafb..2fdb30ed 100644
--- a/coast/diagnostics/tidegauge_analysis.py
+++ b/coast/diagnostics/tidegauge_analysis.py
@@ -23,11 +23,11 @@ def match_missing_values(cls, data_array1, data_array2, fill_value=np.nan):
         Will match any missing values between two tidegauge_multiple datasets.
         Where missing values (defined by fill_value) are found in either dataset
         they are also placed in the corresponding location in the other dataset.
-        Returns two new tidegaugeMultiple objects containing the new
-        ssh data. Datasets must contain ssh variables and only ssh will be
-        masked.
+        Returns two new tidegauge objects containing only the new
+        masked data arrays.
         """
 
+        # Ensuring arrays are dimensionally aligned
         if data_array2.dims[0] == "t_dim":
             data_array2 = data_array2.transpose()
         if data_array1.dims[0] == "t_dim":
@@ -173,7 +173,7 @@ def calculate_non_tidal_residuals(
         Calculate non tidal residuals by subtracting values in data_array_tide
         from data_array_ssh. You may optionally apply a filter to the non
         tidal residual data by setting apply_filter = True. This uses the
-        scipy.signal.savgol_filter function, which you ay pass window_length
+        scipy.signal.savgol_filter function, which you may pass window_length
         and poly_order.
         """
 
@@ -181,11 +181,20 @@ def calculate_non_tidal_residuals(
         ntr = data_array_ssh - data_array_tide
         n_port = data_array_ssh.sizes["id_dim"]
 
+        if ntr.dims[0] == "t_dim" and ntr.dims[1] == "id_dim":
+            ntr = ntr.transpose()
+
         # NTR: Apply filter if wanted
         if apply_filter:
             for pp in range(n_port):
-                ntr[pp, :] = signal.savgol_filter(ntr[pp, :], 25, 3)
-
+                y = ntr[pp, :]
+                try:
+                    # interpolate over nans
+                    nans, x = general_utils.nan_helper(y)  # location interior nans
+                    y[nans] = np.interp(x(nans), x(~nans), y[~nans])  # interpolate
+                    ntr[pp, :] = signal.savgol_filter(y, window_length, polyorder)
+                except:
+                    ntr[pp, :] = np.nan
         # Create output Tidegauge object and return
         ds_coords = data_array_ssh.coords.to_dataset()
         ds_coords["ntr"] = ntr
diff --git a/environment.yml b/environment.yml
index 5a502731..88991eab 100644
--- a/environment.yml
+++ b/environment.yml
@@ -10,9 +10,9 @@ dependencies:
   - gsw==3.4.0
   - pip:
     - numpy==1.22.3
-    - dask==2022.3.0
-    - dask[complete]==2022.3.0
-    - xarray==2022.3.0
+    - dask>=2022.3.0
+    - dask[complete]>=2022.3.0
+    - xarray>=2022.3.0
     - matplotlib>=3.5.3
     - netCDF4==1.5.8
     - scipy>=1.8.0
diff --git a/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb b/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb
new file mode 100644
index 00000000..504bca24
--- /dev/null
+++ b/example_scripts/notebooks/runnable_notebooks/general/mask_maker_tutorial.ipynb
@@ -0,0 +1,491 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5eca7994-6fa1-44e1-b95c-fc8a0fecf7bd",
+   "metadata": {},
+   "source": [
+    "A demonstration of the MaskMaker class to build and use regional masking\n",
+    "\n",
+    "MaskMasker is a class of methods to assist with making regional masks within COAsT.\n",
+    "Presently the mask generated are external to MaskMaker.\n",
+    "Masks are constructed as gridded boolean numpy array for each region, which are stacked over a dim_mask dimension. \n",
+    "The mask arrays are generated on a supplied horizontal grid. The masks are then stored in xarray objects along with regions names.\n",
+    "\n",
+    "Examples are given working with Gridded and Profile data."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "14277e0d-4dbc-4e0f-b3a2-6853dca66d46",
+   "metadata": {},
+   "source": [
+    "### Relevant imports and filepath configuration"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4773751-3544-4ebd-a795-cfe128b70743",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import coast\n",
+    "import numpy as np\n",
+    "from os import path\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.colors as colors  # colormap fiddling\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "780605fd-ae53-4ec5-b7fd-80b2a2ee07ea",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set some paths\n",
+    "root = \"./\"\n",
+    "dn_files = root + \"./example_files/\"\n",
+    "fn_nemo_grid_t_dat = dn_files + \"nemo_data_T_grid_Aug2015.nc\"\n",
+    "fn_nemo_dom = dn_files + \"coast_example_nemo_domain.nc\"\n",
+    "config_t = root + \"./config/example_nemo_grid_t.json\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5d3f6987-f05d-4a54-a932-e4bbf84becb1",
+   "metadata": {},
+   "source": [
+    "### Loading data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7677050c-775d-4172-9561-61c3c89aa77b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a Gridded object and load in the data:\n",
+    "nemo = coast.Gridded(fn_nemo_grid_t_dat, fn_nemo_dom, config=config_t)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c4572cc6",
+   "metadata": {},
+   "source": [
+    "# Initialise MaskMaker and define target grid\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21b40170",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mm = coast.MaskMaker()\n",
+    "\n",
+    "# Define Regional Masks\n",
+    "regional_masks = []\n",
+    "\n",
+    "# Define convenient aliases based on nemo data\n",
+    "lon = nemo.dataset.longitude.values\n",
+    "lat = nemo.dataset.latitude.values\n",
+    "bathy = nemo.dataset.bathymetry.values\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "91dc27a4",
+   "metadata": {},
+   "source": [
+    "# Use MaskMaker to define new regions\n",
+    "\n",
+    "MaskMaker can build a stack of boolean masks in an xarray dataset for regional analysis. Regions can be supplied by providing vertices coordiates to the `make_region_from_vertices` method. (Vertices coordinates can be passed as xarray objects or as numpy arrays).\n",
+    "The method returns a numpy array of booleans."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b151cd2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Draw and fill a square\n",
+    "vertices_lon = [-5, -5, 5, 5]\n",
+    "vertices_lat = [40, 60, 60, 40]\n",
+    "\n",
+    "# input lat/lon as xr.DataArray\n",
+    "filled1 = mm.make_region_from_vertices(nemo.dataset.longitude, nemo.dataset.latitude, vertices_lon, vertices_lat)\n",
+    "# input lat/lon as np.ndarray\n",
+    "filled2 = mm.make_region_from_vertices(\n",
+    "    nemo.dataset.longitude.values, nemo.dataset.latitude.values, vertices_lon, vertices_lat\n",
+    ")\n",
+    "\n",
+    "check = (filled1 == filled2).all()\n",
+    "print(f\"numpy array outputs are the same? {check}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3c6ed915",
+   "metadata": {},
+   "source": [
+    "The boolean numpy array can then be converted to an xarray object using `make_mask_dataset()` for improved interactions with other xarray objects. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c978714c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_xr = mm.make_mask_dataset(nemo.dataset.longitude.values, nemo.dataset.latitude.values, filled1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2189f08d",
+   "metadata": {},
+   "source": [
+    "# Use MaskMaker for predefined regions\n",
+    "\n",
+    "The NWS has a number of predefined regions. These are numpy boolean arrays as functions of the specified latitude, longitude and bathymetry. They can be appended into a list of arrays, which can be similarly converted into an xarray object."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96c64f90",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "masks_list = []\n",
+    "\n",
+    "# Add regional mask for whole domain\n",
+    "masks_list.append(np.ones(lon.shape))\n",
+    "\n",
+    "# Add regional mask for English Channel\n",
+    "masks_list.append(mm.region_def_nws_north_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_outer_shelf(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_norwegian_trench(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_nws_english_channel(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_south_north_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_off_shelf(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_irish_sea(lon, lat, bathy))\n",
+    "masks_list.append(mm.region_def_kattegat(lon, lat, bathy))\n",
+    "\n",
+    "masks_names = [\"whole domain\", \"north sea\", \"outer shelf\", \"norwegian trench\",\n",
+    "                \"english_channel\", \"southern north sea\", \"off shelf\",\n",
+    "                \"irish sea\", \"kattegat\",]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ed5d7ddf",
+   "metadata": {},
+   "source": [
+    "As before the numpy arrays (here as a list) can be converted into an xarray dataset where each mask is separated along the `dim_mask` dimension"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c43565af",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_xr = mm.make_mask_dataset(lon, lat, masks_list, masks_names)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "15bb0838",
+   "metadata": {
+    "pycharm": {
+     "name": "#%% md\n"
+    }
+   },
+   "source": [
+    "Inspect mask xarray object structure\n",
+    "\n",
+    "``mask_xr``"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68a07b35",
+   "metadata": {},
+   "source": [
+    "## Plot masks\n",
+    "\n",
+    "Inspect the mask with a `quick_plot()` method."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9daf7c40",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mm.quick_plot(mask_xr)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "454bdeb4",
+   "metadata": {},
+   "source": [
+    "NB overlapping regions are not given special treatment, the layers are blindly superimposed on each other. E.g. as demonstrated with \"Norwegian Trench\" and \"off shelf\", or \"whole domain\" and any other region."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c1217563",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.subplot(2,2,1)\n",
+    "mm.quick_plot(mask_xr.sel(dim_mask=[0,3]))\n",
+    "\n",
+    "plt.subplot(2,2,2)\n",
+    "mm.quick_plot(mask_xr.sel(dim_mask=[1,2,4,5,6,7,8]))\n",
+    "\n",
+    "plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7e4a3a6f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Show overlap\n",
+    "mask_xr.mask.sum(dim='dim_mask').plot(levels=(1,2,3,4))\n",
+    "\n",
+    "# Save if required\n",
+    "#plt.savefig('tmp.png')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dcbd7b05",
+   "metadata": {},
+   "source": [
+    "# Regional analysis with Profile data\n",
+    "\n",
+    "Apply the regional masks to average SST"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d34185a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Read EN4 data into profile object\n",
+    "fn_prof = path.join(dn_files, \"coast_example_en4_201008.nc\")\n",
+    "fn_cfg_prof = path.join(\"config\",\"example_en4_profiles.json\")\n",
+    "profile = coast.Profile(config=fn_cfg_prof)\n",
+    "profile.read_en4( fn_prof )\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "181cdef2",
+   "metadata": {},
+   "source": [
+    "Then we use `ProfileAnalysis.determine_mask_indices()` to figure out which profiles in a Profile object lie within each regional mask:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f57e1a1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis = coast.ProfileAnalysis()\n",
+    "mask_indices = analysis.determine_mask_indices(profile, mask_xr)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "49ef3d40",
+   "metadata": {},
+   "source": [
+    "This returns an object called `mask_indices`, which is required to pass to `ProfileAnalysis.mask_means()`. This routine will return a new xarray dataset containing averaged data for each region:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "74900bf2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "profile_mask_means = analysis.mask_means(profile, mask_indices)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69f6cdb0",
+   "metadata": {},
+   "source": [
+    "This routine operates over all variables in the `profile` object. It calculates means by region preserving depth information (`profile_mean_*`) and also averaging over depth information (`all_mean_*`). The variables are returned with these prefixes accordingly. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4ab174e8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "profile_mask_means"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4e70c0df",
+   "metadata": {},
+   "source": [
+    "Notices that the number of mask dimensions is not necessarily preserved between the mask and the mask averaged variables. This happens if, for example, there are no profiles in one of the mask regions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c756a11c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "check1 = mask_indices.dims[\"dim_mask\"] == profile_mask_means.dims[\"dim_mask\"]\n",
+    "print(check1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5a033719",
+   "metadata": {},
+   "source": [
+    "The mean profiles can be visualised or further processed (notice the Irish Sea region is missing because there were no profiles in the example dataset)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "749b6e8a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for count_region in range(profile_mask_means.sizes['dim_mask']):    \n",
+    "    plt.plot( \n",
+    "            profile_mask_means.profile_mean_temperature.isel(dim_mask=count_region),\n",
+    "            profile_mask_means.profile_mean_depth.isel(dim_mask=count_region),\n",
+    "            label=profile_mask_means.region_names[count_region].values,\n",
+    "            marker=\".\", linestyle='none')\n",
+    "\n",
+    "plt.ylim([10,1000])\n",
+    "plt.yscale(\"log\")\n",
+    "plt.gca().invert_yaxis()\n",
+    "plt.xlabel('temperature'); plt.ylabel('depth')\n",
+    "plt.legend()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3b8544cc",
+   "metadata": {},
+   "source": [
+    "# Regional analysis with Gridded data\n",
+    "\n",
+    "Apply the regional masks to average SST. This is done manually as there are not yet COAsT methods to broadcast the operations across all variables."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4b72009d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Syntax: xr.where(if <first>, then <2nd>, else <3rd>) \n",
+    "mask_SST = xr.where( mask_xr.mask, nemo.dataset.temperature.isel(z_dim=0), np.NaN)\n",
+    "\n",
+    "# Take the mean over space for each region\n",
+    "mask_mean_SST = mask_SST.mean(dim=\"x_dim\").mean(dim=\"y_dim\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8257a975",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Inspect the processed data\n",
+    "mask_mean_SST.plot()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ad05f3e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot timeseries per region\n",
+    "\n",
+    "for count_region in range(mask_mean_SST.sizes['dim_mask']):\n",
+    "    \n",
+    "    plt.plot( \n",
+    "        mask_mean_SST.isel(dim_mask=count_region),\n",
+    "        label=mask_mean_SST.region_names[count_region].values,\n",
+    "        marker=\".\", linestyle='none')\n",
+    "\n",
+    "plt.xlabel('time'); plt.ylabel('SST')\n",
+    "plt.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73fe919e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/example_scripts/notebooks/runnable_notebooks/gridded/potential_energy_tutorial.ipynb b/example_scripts/notebooks/runnable_notebooks/gridded/potential_energy_tutorial.ipynb
index 5d684259..f5dc1d40 100644
--- a/example_scripts/notebooks/runnable_notebooks/gridded/potential_energy_tutorial.ipynb
+++ b/example_scripts/notebooks/runnable_notebooks/gridded/potential_energy_tutorial.ipynb
@@ -5,7 +5,7 @@
    "id": "5eca7994-6fa1-44e1-b95c-fc8a0fecf7bd",
    "metadata": {},
    "source": [
-    "A demonstration to calculate the Potential Energy Anomaly.\n"
+    "A demonstration to calculate the Potential Energy Anomaly and demonstrate regional masking with MaskMaker\n"
    ]
   },
   {
@@ -27,7 +27,8 @@
     "import numpy as np\n",
     "import os\n",
     "import matplotlib.pyplot as plt\n",
-    "import matplotlib.colors as colors  # colormap fiddling"
+    "import matplotlib.colors as colors  # colormap fiddling\n",
+    "import xarray as xr"
    ]
   },
   {
@@ -61,7 +62,7 @@
    "outputs": [],
    "source": [
     "# Create a Gridded object and load in the data:\n",
-    "nemo_t = coast.Gridded(fn_nemo_grid_t_dat, fn_nemo_dom, config=config_t)"
+    "nemo = coast.Gridded(fn_nemo_grid_t_dat, fn_nemo_dom, config=config_t)"
    ]
   },
   {
@@ -85,12 +86,20 @@
    "outputs": [],
    "source": [
     "# Compute a vertical max to exclude depths below 200m\n",
-    "Zd_mask, kmax, Ikmax = nemo_t.calculate_vertical_mask(200.)\n",
+    "Zd_mask, kmax, Ikmax = nemo.calculate_vertical_mask(200.)\n",
     "\n",
     "# Initiate a stratification diagnostics object\n",
-    "strat = coast.GriddedStratification(nemo_t)"
+    "strat = coast.GriddedStratification(nemo)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d342e9b0",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -99,7 +108,7 @@
    "outputs": [],
    "source": [
     "# calculate PEA for unmasked depths\n",
-    "strat.calc_pea(nemo_t, Zd_mask)"
+    "strat.calc_pea(nemo, Zd_mask)"
    ]
   },
   {
@@ -130,7 +139,191 @@
    "id": "a8b2bf5b",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "strat.dataset"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c4572cc6",
+   "metadata": {},
+   "source": [
+    "# Use MaskMaker to define regions and do regional analysis\n",
+    "\n",
+    "MaskMaker can build a stack of boolean masks in an xarray dataset for regional analysis. For the NWS we can use some built-in regions.  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96c64f90",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mm = coast.MaskMaker()\n",
+    "\n",
+    "# Define Regional Masks\n",
+    "regional_masks = []\n",
+    "\n",
+    "# Define convenient aliases based on nemo data\n",
+    "lon = nemo.dataset.longitude.values\n",
+    "lat = nemo.dataset.latitude.values\n",
+    "bathy = nemo.dataset.bathymetry.values\n",
+    "\n",
+    "# Add regional mask for whole domain\n",
+    "regional_masks.append(np.ones(lon.shape))\n",
+    "\n",
+    "# Add regional mask for English Channel\n",
+    "regional_masks.append(mm.region_def_nws_north_sea(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_nws_outer_shelf(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_nws_norwegian_trench(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_nws_english_channel(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_south_north_sea(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_off_shelf(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_irish_sea(lon, lat, bathy))\n",
+    "regional_masks.append(mm.region_def_kattegat(lon, lat, bathy))\n",
+    "\n",
+    "region_names = [\"whole domain\", \"north sea\", \"outer shelf\", \"norwegian trench\",\n",
+    "                \"english_channel\", \"southern north sea\", \"off shelf\",\n",
+    "                \"irish sea\", \"kattegat\",]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ed5d7ddf",
+   "metadata": {},
+   "source": [
+    "Convert this list of masks into a dataset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c43565af",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_list = mm.make_mask_dataset(lon, lat, regional_masks, region_names)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15bb0838",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_list"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68a07b35",
+   "metadata": {},
+   "source": [
+    "Inspect the mask with a `quick_plot()` method. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9daf7c40",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mm.quick_plot(mask_list)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "454bdeb4",
+   "metadata": {},
+   "source": [
+    "NB overlapping regions are not given special treatment, the layers are blindly superimposed on each other. E.g. as demonstrated with \"Norwegian Trench\" and \"off shelf\", or \"whole domain\" and any other region."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c1217563",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "plt.subplot(2,2,1)\n",
+    "mm.quick_plot(mask_list.sel(dim_mask=[0,3]))\n",
+    "\n",
+    "plt.subplot(2,2,2)\n",
+    "mm.quick_plot(mask_list.sel(dim_mask=[1,2,4,5,6,7,8]))\n",
+    "\n",
+    "\n",
+    "\n",
+    "plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7e4a3a6f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Show overlap\n",
+    "mask_list.mask.sum(dim='dim_mask').plot( levels=(1,2,3,4))\n",
+    "\n",
+    "# Save if required\n",
+    "#plt.savefig('tmp.png')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3b8544cc",
+   "metadata": {},
+   "source": [
+    "# Regional analysis\n",
+    "\n",
+    "Average stratification object over regions using the mask"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4b72009d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_means = (strat.dataset*mask_list.mask).mean(dim='x_dim').mean(dim='y_dim')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ac2a67a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_means"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ad05f3e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot timeseries per region\n",
+    "\n",
+    "for count_region in range(mask_means.dims['dim_mask']):\n",
+    "    \n",
+    "    plt.plot( \n",
+    "        mask_means.PEA.isel(dim_mask=count_region),\n",
+    "        label=mask_means.region_names[count_region].values,\n",
+    "        marker=\".\", linestyle='none')\n",
+    "\n",
+    "plt.xlabel('time'); plt.ylabel('PEA')\n",
+    "plt.legend()"
+   ]
   }
  ],
  "metadata": {
@@ -149,7 +342,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.8.10"
   }
  },
  "nbformat": 4,
diff --git a/example_scripts/notebooks/runnable_notebooks/profile/0._profile_introduction.ipynb b/example_scripts/notebooks/runnable_notebooks/profile/0._profile_introduction.ipynb
index 2dd2e4c7..6d740619 100644
--- a/example_scripts/notebooks/runnable_notebooks/profile/0._profile_introduction.ipynb
+++ b/example_scripts/notebooks/runnable_notebooks/profile/0._profile_introduction.ipynb
@@ -1103,4 +1103,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}
\ No newline at end of file
diff --git a/example_scripts/notebooks/runnable_notebooks/tidegauge/tidegauge_validation_tutorial.ipynb b/example_scripts/notebooks/runnable_notebooks/tidegauge/tidegauge_validation_tutorial.ipynb
index 210034a1..a53d6bb1 100644
--- a/example_scripts/notebooks/runnable_notebooks/tidegauge/tidegauge_validation_tutorial.ipynb
+++ b/example_scripts/notebooks/runnable_notebooks/tidegauge/tidegauge_validation_tutorial.ipynb
@@ -233,13 +233,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "stn_id=1  # pick a gauge station\n",
+    "stn_id=26  # pick a gauge station\n",
     "\n",
     "plt.subplot(2,1,1)\n",
     "obs.dataset.ssh[stn_id].plot()  # rename time dimension to enable automatic x-axis labelling\n",
     "\n",
     "plt.subplot(2,1,2)\n",
-    "obs.dataset.ssh[stn_id].rename({'t_dim':'time'}).plot()  # rename time dimension to enable automatic x-axis labelling\n",
+    "obs.dataset.ssh[stn_id].swap_dims({'t_dim':'time'}).plot()  # rename time dimension to enable automatic x-axis labelling\n",
     "plt.tight_layout()\n",
     "plt.show()"
    ]
@@ -370,7 +370,7 @@
    "id": "ad8d0be3-7941-40ba-af6f-ef747ae2912d",
    "metadata": {},
    "source": [
-    "However, the new `tidegauge_from_model` will the same number of time entries as the `obs` data (rather than the model data). So, for a more useful demonstration we trim the observed gauge data so it better matches the model data."
+    "However, the new `tidegauge_from_model` will the same number of time entries as the `obs` data, rather than the model data (so this will include lots of empty values). So, for a more useful demonstration we trim the observed gauge data so it better matches the model data."
    ]
   },
   {
@@ -407,6 +407,24 @@
     "model_timeseries.dataset"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "91976034",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stn_id=26  # pick a gauge station\n",
+    "\n",
+    "plt.subplot(2,1,1)\n",
+    "model_timeseries.dataset.ssh.isel(id_dim=stn_id).plot()  # rename time dimension to enable automatic x-axis labelling\n",
+    "\n",
+    "plt.subplot(2,1,2)\n",
+    "obs.dataset.ssh.isel(id_dim=stn_id).swap_dims({'t_dim':'time'}).plot()  # rename time dimension to enable automatic x-axis labelling\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "88736e30",
@@ -444,9 +462,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# In this case, transpose the interpolated dataset\n",
-    "model_timeseries.dataset = model_timeseries.dataset.transpose()\n",
-    "\n",
     "# This routine searches for missing values in each dataset and applies them\n",
     "# equally to each corresponding dataset\n",
     "obs_new, model_new = tganalysis.match_missing_values(obs.dataset.ssh, model_timeseries.dataset.ssh)"
@@ -581,7 +596,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "d6439a17",
+   "id": "bfdcaad6",
    "metadata": {},
    "source": [
     "Get new TidegaugeMultiple objects containing non tidal residuals:"
@@ -594,11 +609,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ntr_mod = tganalysis.calculate_non_tidal_residuals(model_new.dataset.ssh, tide_mod.dataset.reconstructed)\n",
-    "ntr_obs = tganalysis.calculate_non_tidal_residuals(obs_new.dataset.ssh, tide_obs.dataset.reconstructed)\n",
+    "ntr_mod = tganalysis.calculate_non_tidal_residuals(model_new.dataset.ssh, tide_mod.dataset.reconstructed, apply_filter=False)\n",
+    "ntr_obs = tganalysis.calculate_non_tidal_residuals(obs_new.dataset.ssh, tide_obs.dataset.reconstructed, apply_filter=True, window_length=10, polyorder=2)\n",
     "\n",
     "# Take a look\n",
-    "ntr_mod.dataset"
+    "ntr_obs.dataset"
    ]
   },
   {
@@ -627,6 +642,14 @@
     "dx0.dataset\n"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "deeefe20",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "markdown",
    "id": "ca8b8ae6",
@@ -734,7 +757,7 @@
    "id": "1c344448",
    "metadata": {},
    "source": [
-    "Note that the non-tidal residual is a noisy timeseries (computed as a difference between two timeseries) so peaks do not necessarily correspond to peaks in total water level. For this reason `time_over_threshold_*` can be useful. Below we see that about 6% (y-axis) of the observations of non-tidal residual exceed 20cm (x-axis). Again, threshold statistics are more useful with a longer record."
+    "Note that the non-tidal residual is a noisy timeseries (computed as a difference between two timeseries) so peaks do not necessarily correspond to peaks in total water level. For this reason `time_over_threshold_*` can be useful. Below we see that about 28% (y-axis) of the observations of non-tidal residual exceed 20cm (x-axis). Again, threshold statistics are more useful with a longer record."
    ]
   },
   {
@@ -747,8 +770,8 @@
     "normalised_event_count = 100 * thresh_obs.time_over_threshold_ntr.isel(id_dim=stn_id) / thresh_obs.time_over_threshold_ntr.isel(id_dim=stn_id).max()\n",
     "\n",
     "plt.plot( thresh_obs.threshold, normalised_event_count )\n",
-    "plt.ylim([0,20])\n",
-    "plt.xlim([0.1,0.4])\n",
+    "plt.ylim([20,40])\n",
+    "plt.xlim([0.0,0.4])\n",
     "plt.title(f\"time over threshold\")\n",
     "plt.xlabel('Analysis threshold (m)')\n",
     "plt.ylabel('normalised event count (%)')\n",
@@ -803,7 +826,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "400086e9",
+   "id": "30d1e18c",
    "metadata": {},
    "outputs": [],
    "source": []
@@ -830,4 +853,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
\ No newline at end of file
+}
diff --git a/setup.py b/setup.py
index cf06a4bd..70a8fcdb 100644
--- a/setup.py
+++ b/setup.py
@@ -2,11 +2,10 @@
 from sys import argv
 from types import SimpleNamespace
 
-
 PACKAGE = SimpleNamespace(
     **{
         "name": "COAsT",
-        "version": "3.1.2",
+        "version": "3.2.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/test_general_utils.py b/unit_testing/test_general_utils.py
index 90773012..99851c26 100644
--- a/unit_testing/test_general_utils.py
+++ b/unit_testing/test_general_utils.py
@@ -7,6 +7,7 @@
 from coast._utils import general_utils
 import unittest
 import numpy as np
+import xarray as xr
 import pytz
 import datetime
 import unit_test_files as files
@@ -43,3 +44,18 @@ def test_bst_to_gmt(self):
             "2020-10-11T11:00:00"
         )
         self.assertTrue(check1, msg="check1")
+
+    def test_nan_helper(self):
+        y = np.array([np.NaN, 1, 1, np.NaN, np.NaN, 7, 2, np.NaN, 0])
+        y_xr = xr.DataArray(y)
+        # numpy array
+        nans, x = general_utils.nan_helper(y)
+        y[nans] = np.interp(x(nans), x(~nans), y[~nans])
+        check1 = all(y == np.array([1.0, 1.0, 1.0, 3.0, 5.0, 7.0, 2.0, 1.0, 0.0]))
+        # xarray
+        nans, x = general_utils.nan_helper(y_xr)
+        y_xr[nans] = np.interp(x(nans), x(~nans), y_xr[~nans])
+        check2 = all(y_xr.values == np.array([1.0, 1.0, 1.0, 3.0, 5.0, 7.0, 2.0, 1.0, 0.0]))
+
+        self.assertTrue(check1, msg="check1")
+        self.assertTrue(check2, msg="check2")
diff --git a/unit_testing/test_maskmaker_methods.py b/unit_testing/test_maskmaker_methods.py
index 56e0cd85..047d174c 100644
--- a/unit_testing/test_maskmaker_methods.py
+++ b/unit_testing/test_maskmaker_methods.py
@@ -6,6 +6,7 @@
 import unittest
 import coast
 import numpy as np
+import matplotlib.pyplot as plt
 import unit_test_files as files
 
 
@@ -49,3 +50,42 @@ def test_fill_polygon_by_lonlat(self):
 
         self.assertTrue(check1, "check1")
         self.assertTrue(check2, "check2")
+
+    def test_make_region_from_vertices(self):
+        sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
+        mm = coast.MaskMaker()
+        # Draw and fill a square
+        vertices_lon = [-5, -5, 5, 5]
+        vertices_lat = [40, 60, 60, 40]
+        # input lat/lon as xr.DataArray
+        filled1 = mm.make_region_from_vertices(sci.dataset.longitude, sci.dataset.latitude, vertices_lon, vertices_lat)
+        # input lat/lon as np.ndarray
+        filled2 = mm.make_region_from_vertices(
+            sci.dataset.longitude.values, sci.dataset.latitude.values, vertices_lon, vertices_lat
+        )
+
+        # TEST: Check some data
+        check1 = filled1[50, 50] == 0 and filled1[50, 150] == 1
+        check2 = filled2[50, 50] == 0 and filled2[50, 150] == 1
+
+        self.assertTrue(check1, "check1")
+        self.assertTrue(check2, "check2")
+
+    def test_make_mask_dataset_and_quick_plot(self):
+        sci = coast.Gridded(files.fn_nemo_dat, files.fn_nemo_dom, config=files.fn_config_t_grid)
+        mm = coast.MaskMaker()
+        # Draw and fill a square
+        vertices_lon = [-5, -5, 5, 5]
+        vertices_lat = [40, 60, 60, 40]
+        # input lat/lon as xr.DataArray
+        filled = mm.make_region_from_vertices(sci.dataset.longitude, sci.dataset.latitude, vertices_lon, vertices_lat)
+
+        mask_list = mm.make_mask_dataset(sci.dataset.longitude.values, sci.dataset.latitude.values, filled)
+
+        check = np.isclose(mask_list.mask.values.sum(), 27300)
+        self.assertTrue(check, "check")
+
+        with self.subTest("MaskMaker quick plot"):
+            mm.quick_plot(mask_list)
+            plt.savefig(files.dn_fig + "maskmaker_quick_plot.png")
+            plt.close("all")
diff --git a/unit_testing/test_profile_methods.py b/unit_testing/test_profile_methods.py
index 6931e036..06080790 100644
--- a/unit_testing/test_profile_methods.py
+++ b/unit_testing/test_profile_methods.py
@@ -116,8 +116,11 @@ def test_compare_processed_profile_with_model(self):
             # Do average differences for each region
             mask_means = pa.mask_means(difference, mask_indices)
 
-            check1 = np.isclose(mask_means.all_mean_diff_temperature.values[0], -0.8169331445866469)
+            check1 = mask_indices.dims["dim_mask"] == mask_means.dims["dim_mask"]  # since profiles are in mask regions
+            check2 = np.isclose(mask_means.all_mean_diff_temperature.values[0], -0.8169331445866469)
+
             self.assertTrue(check1, "check1")
+            self.assertTrue(check2, "check2")
 
         with self.subTest("Surface/Bottom averaging"):
             surface = 5
diff --git a/unit_testing/unit_test_contents.txt b/unit_testing/unit_test_contents.txt
index f58fb1ea..5d250880 100755
--- a/unit_testing/unit_test_contents.txt
+++ b/unit_testing/unit_test_contents.txt
@@ -18,6 +18,7 @@
       c. copy_coast_object
       d. day_of_week
       e. getitem
+      f. nan_helper
 
 3. test_gridded_harmonics
       a. combine_and_convert_harmonics
@@ -86,6 +87,8 @@
 16. test_maskmaker_methods
       a. fill_polygon_by_index
       b. fill_polygon_by_lonlat
+      c. make_mask_dataset_and_quick_plot
+      d. make_region_from_vertices
 
 17. test_climatology
       a. monthly_and_seasonal_climatology