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

Issue:
restructured coord_var to be consistent between gridded/indexed
Fix for coord_var consistency. Moved out of static vairables and into dataset instead.


Discussion / hints:

Fix patch:
diff --git a/coast/__init__.py b/coast/__init__.py
index 09f6c868..a74fb7c7 100644
--- a/coast/__init__.py
+++ b/coast/__init__.py
@@ -15,5 +15,15 @@
 from . import general_utils
 from . import plot_util
 from . import crps_util
+from .index import Indexed
+from .track import Track
+from .lagrangian import Lagrangian
+from .oceanparcels import Oceanparcels
+from .glider import Glider
+from .argos import Argos
+from .altimetry import Altimetry
+from .timeseries import Timeseries
+from .tidegauge import Tidegauge
+from .profile import Profile
 from .contour import Contour, ContourF, ContourT
 from .eof import compute_eofs, compute_hilbert_eofs
diff --git a/coast/altimetry.py b/coast/altimetry.py
index 18f8c245..59a9120f 100644
--- a/coast/altimetry.py
+++ b/coast/altimetry.py
@@ -1,12 +1,15 @@
-from .coast import Coast
+"""Altimetry class"""
+from .track import Track
 import numpy as np
 import xarray as xr
-from .logging_util import get_slug, debug, info
 import sklearn.metrics as metrics
 from . import general_utils, plot_util, crps_util
+from .logging_util import get_slug, debug, info, warn, warning
+from typing import Union
+from pathlib import Path
 
 
-class Altimetry(Coast):  # TODO All abstract methods should be implemented
+class Altimetry(Track):
     """
     An object for reading, storing and manipulating altimetry data.
     Currently the object contains functionality for reading altimetry netCDF
@@ -17,9 +20,9 @@ class Altimetry(Coast):  # TODO All abstract methods should be implemented
     * Date Format Overview *
 
         1. A single dimension (time).
-        2. Three coordinates: time, latitude, longitude. All lie on the 't_dim'
+        2. Three coordinates: time, latitude, longitude. All lie on the 'time'
            dimension.
-        3. Observed variable DataArrays on the t_dim dimension.
+        3. Observed variable DataArrays on the time dimension.
 
     There are currently no naming conventions for the variables however
     examples from the CMEMS database include sla_filtered, sla_unfiltered and
@@ -49,33 +52,71 @@ class Altimetry(Coast):  # TODO All abstract methods should be implemented
 
     """
 
-    ##############################################################################
-    ###                ~ Initialisation and File Reading ~                     ###
-    ##############################################################################
+    def __init__(self, file_path: str = None, multiple=False, config: Union[Path, str] = None):
+        """Initialization and file reading.
 
-    def __init__(self, file=None, chunks: dict = None, multiple=False):
+        Args:
+            file_path (str): path to data file
+            multiple (boolean): True if reading multiple files otherwise False
+            config (Union[Path, str]): path to json config file.
+        """
         debug(f"Creating a new {get_slug(self)}")
-        if file is not None:
-            self.read_cmems(file, chunks, multiple)
+        super().__init__(config)
+
+        if file_path is None:
+            warn("Object created but no file or directory specified: \n" "{0}".format(str(self)), UserWarning)
         else:
-            self.dataset = None
-        debug(f"{get_slug(self)} initialised")
-        return  # TODO Super __init__ should be called at some point
+            self.read_cmems(file_path, multiple)
+            self.apply_config_mappings()
+
+        print(f"{get_slug(self)} initialised")
 
-    def read_cmems(self, file, chunks, multiple):
-        """Reads altimetry data from a CMEMS netcdf file. Calls Coast.init()
-        to make use of its load methods"""
-        super().__init__(file, chunks, multiple)
-        self.dataset = self.dataset.rename_dims(self.dim_mapping)
-        # self.dataset.attrs = {}
+    def read_cmems(self, file_path: str, multiple) -> None:
+        """Read file.
+
+        Args:
+            file_path (str): path to data file
+            multiple (boolean): True if reading multiple files otherwise False
+        """
+        self.load(file_path, self.chunks, multiple)
 
-    def set_dimension_mapping(self):
-        self.dim_mapping = {"time": "t_dim"}  # TODO Object attributes should be set in __init__
-        debug(f"{get_slug(self)} dim_mapping set to {self.dim_mapping}")
+    def load(self, file_or_dir, chunks: dict = None, multiple=False) -> None:
+        """
+        Loads a file into a object's dataset variable using xarray
+
+        Args:
+            file_or_dir (str) : file name or directory to multiple files.
+            chunks (dict)     : Chunks to use in Dask [default None]
+            multiple (bool)   : If true, load in multiple files from directory.
+                                If false load a single file [default False]
+        """
+        if multiple:
+            self.load_multiple(file_or_dir, chunks)
+        else:
+            self.load_single(file_or_dir, chunks)
 
-    def set_variable_mapping(self):
-        self.var_mapping = None  # TODO Object attributes should be set in __init__
-        debug(f"{get_slug(self)} var_mapping set to {self.var_mapping}")
+    def __getitem__(self, name: str):
+        return self.dataset[name]
+
+    def load_single(self, file, chunks: dict = None) -> None:
+        """Loads a single file into object's dataset variable.
+
+        Args:
+            file (str) : file name or directory to multiple files.
+            chunks (dict) : Chunks to use in Dask [default None]
+        """
+        self.dataset = xr.open_dataset(file, chunks=chunks)
+
+    def load_multiple(self, directory_to_files, chunks: dict = None) -> None:
+        """Loads multiple files from directory into dataset variable.
+
+        Args:
+            directory_to_files (str) : directory path to multiple files.
+            chunks (dict) : Chunks to use in Dask [default None]
+        """
+        self.dataset = xr.open_mfdataset(
+            directory_to_files, chunks=chunks, parallel=True, combine="by_coords"
+        )  # , compat='override')
 
     def subset_indices_lonlat_box(self, lonbounds, latbounds):
         """Generates array indices for data which lies in a given lon/lat box.
@@ -88,9 +129,9 @@ def subset_indices_lonlat_box(self, lonbounds, latbounds):
 
         return: Indices corresponding to datapoints inside specified box
         """
-        debug(f"Subsetting {get_slug(self)} indices in {lonbounds}, {latbounds}")
-        lon = self.dataset.longitude.copy()
-        lat = self.dataset.latitude
+        print(f"Subsetting {get_slug(self)} indices in {lonbounds}, {latbounds}")
+        lon = self.dataset.longitude.copy().values
+        lat = self.dataset.latitude.values
         lon[lon > 180] = lon[lon > 180] - 360
         lon[lon < -180] = lon[lon < -180] + 360
         ff1 = (lon > lonbounds[0]).astype(int)  # FIXME This should fail? We can just treat bools as ints here...
@@ -100,12 +141,10 @@ def subset_indices_lonlat_box(self, lonbounds, latbounds):
         indices = np.where(ff1 * ff2 * ff3 * ff4)
         return indices[0]
 
-    ##############################################################################
-    ###                ~            Plotting             ~                     ###
-    ##############################################################################
+    """======================= Plotting ======================="""
 
     def quick_plot(self, color_var_str: str = None):
-        """"""
+        """Quick plot"""
 
         if color_var_str is not None:
             color_var = self.dataset[color_var_str]
@@ -118,9 +157,7 @@ def quick_plot(self, color_var_str: str = None):
         info("Plot ready, displaying!")
         return fig, ax
 
-    ##############################################################################
-    ###                ~        Model Comparison         ~                     ###
-    ##############################################################################
+    """======================= MODEL comparison routines ======================="""
 
     def obs_operator(self, model, mod_var_name: str, time_interp="nearest", model_mask=None):
         """
@@ -138,7 +175,7 @@ def obs_operator(self, model, mod_var_name: str, time_interp="nearest", model_ma
         Parameters
         ----------
         model : model object (e.g. NEMO)
-        mod_var_name: variable name string to use from model object
+        mod_var: variable name string to use from model object
         time_interp: time interpolation method (optional, default: 'nearest')
             This can take any string scipy.interpolate would take. e.g.
             'nearest', 'linear' or 'cubic'
@@ -152,12 +189,15 @@ def obs_operator(self, model, mod_var_name: str, time_interp="nearest", model_ma
         Adds a DataArray to self.dataset, containing interpolated values.
         """
 
-        debug(f'Interpolating {get_slug(model)} "{mod_var_name}" with time_interp "{time_interp}"')
+        print(f'Interpolating {get_slug(model)} "{mod_var_name}" with time_interp "{time_interp}"')
 
         # Determine mask
         if model_mask == "bathy":
             model_mask = model.dataset.bathymetry.values == 0
 
+        # Get data arrays
+        mod_var_array = model.dataset[mod_var_name]
+
         # Get data arrays
         mod_var = model.dataset[mod_var_name]
 
@@ -235,23 +275,23 @@ def crps(
         if create_new_object:
             new_object = Altimetry()
             new_dataset = self.dataset[["longitude", "latitude", "time"]]
-            new_dataset["crps"] = ("t_dim", crps_list)
-            new_dataset["crps_n_model_pts"] = ("t_dim", n_model_pts)
-            new_dataset["crps_contains_land"] = ("t_dim", contains_land)
+            new_dataset["crps"] = (("t_dim"), crps_list)
+            new_dataset["crps_n_model_pts"] = (("t_dim"), n_model_pts)
+            new_dataset["crps_contains_land"] = (("t_dim"), contains_land)
             new_object.dataset = new_dataset
             return new_object
         else:
-            self.dataset["crps"] = ("t_dim", crps_list)
-            self.dataset["crps_n_model_pts"] = ("t_dim", n_model_pts)
-            self.dataset["crps_contains_land"] = ("t_dim", contains_land)
+            self.dataset["crps"] = (("t_dim"), crps_list)
+            self.dataset["crps_n_model_pts"] = (("t_dim"), n_model_pts)
+            self.dataset["crps_contains_land"] = (("t_dim"), contains_land)
 
     def difference(self, var_str0: str, var_str1: str, date0=None, date1=None):
         """Difference two variables defined by var_str0 and var_str1 between
         two dates date0 and date1. Returns xr.DataArray"""
         var0 = self.dataset[var_str0]
         var1 = self.dataset[var_str1]
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1).values
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1).values
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
         diff = var0 - var1
         return xr.DataArray(diff, dims="t_dim", name="error", coords={"time": self.dataset.time})
 
@@ -260,8 +300,8 @@ def absolute_error(self, var_str0, var_str1, date0=None, date1=None):
         between two dates date0 and date1. Return xr.DataArray"""
         var0 = self.dataset[var_str0]
         var1 = self.dataset[var_str1]
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1).values
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1).values
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
         adiff = np.abs(var0 - var1)
         return xr.DataArray(adiff, dims="t_dim", name="absolute_error", coords={"time": self.dataset.time})
 
@@ -270,8 +310,8 @@ def mean_absolute_error(self, var_str0, var_str1, date0=None, date1=None):
         var_str1 between two dates date0 and date1. Return xr.DataArray"""
         var0 = self.dataset[var_str0]
         var1 = self.dataset[var_str1]
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1).values
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1).values
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
         mae = metrics.mean_absolute_error(var0, var1)
         return mae
 
@@ -280,21 +320,21 @@ def root_mean_square_error(self, var_str0, var_str1, date0=None, date1=None):
         var_str1 between two dates date0 and date1. Return xr.DataArray"""
         var0 = self.dataset[var_str0]
         var1 = self.dataset[var_str1]
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1).values
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1).values
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
         rmse = metrics.mean_squared_error(var0, var1)
         return np.sqrt(rmse)
 
     def time_mean(self, var_str, date0=None, date1=None):
         """Time mean of variable var_str between dates date0, date1"""
         var = self.dataset[var_str]
-        var = general_utils.dataarray_time_slice(var, date0, date1)
+        var = general_utils.data_array_time_slice(var, date0, date1)
         return np.nanmean(var)
 
     def time_std(self, var_str, date0=None, date1=None):
         """Time st. dev of variable var_str between dates date0 and date1"""
         var = self.dataset[var_str]
-        var = general_utils.dataarray_time_slice(var, date0, date1)
+        var = general_utils.data_array_time_slice(var, date0, date1)
         return np.nanstd(var)
 
     def time_correlation(self, var_str0, var_str1, date0=None, date1=None, method="pearson"):
@@ -304,8 +344,8 @@ def time_correlation(self, var_str0, var_str1, date0=None, date1=None, method="p
         var1 = self.dataset[var_str1]
         var0 = var0.rename("var1")
         var1 = var1.rename("var2")
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1)
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1)
+        var0 = general_utils.data_array_time_slice(var0, date0, date1)
+        var1 = general_utils.data_array_time_slice(var1, date0, date1)
         pdvar = xr.merge((var0, var1))
         pdvar = pdvar.to_dataframe()
         corr = pdvar.corr(method=method)
@@ -318,8 +358,8 @@ def time_covariance(self, var_str0, var_str1, date0=None, date1=None):
         var1 = self.dataset[var_str1]
         var0 = var0.rename("var1")
         var1 = var1.rename("var2")
-        var0 = general_utils.dataarray_time_slice(var0, date0, date1)
-        var1 = general_utils.dataarray_time_slice(var1, date0, date1)
+        var0 = general_utils.data_array_time_slice(var0, date0, date1)
+        var1 = general_utils.data_array_time_slice(var1, date0, date1)
         pdvar = xr.merge((var0, var1))
         pdvar = pdvar.to_dataframe()
         cov = pdvar.cov()
@@ -333,28 +373,28 @@ def basic_stats(self, var_str0, var_str1, date0=None, date1=None, create_new_obj
         then this method returns a new ALTIMETRY object containing statistics,
         otherwise variables are saved to the dateset inside this object."""
 
-        difference = self.difference(var_str0, var_str1, date0, date1)
-        absolute_error = self.absolute_error(var_str0, var_str1, date0, date1)
-        mean_absolute_error = self.mean_absolute_error(var_str0, var_str1, date0, date1)
-        root_mean_square_error = self.root_mean_square_error(var_str0, var_str1, date0, date1)
-        correlation = self.time_correlation(var_str0, var_str1, date0, date1)
-        covariance = self.time_covariance(var_str0, var_str1, date0, date1)
+        diff = self.difference(var_str0, var_str1, date0, date1)
+        ae = self.absolute_error(var_str0, var_str1, date0, date1)
+        mae = self.mean_absolute_error(var_str0, var_str1, date0, date1)
+        rmse = self.root_mean_square_error(var_str0, var_str1, date0, date1)
+        corr = self.time_correlation(var_str0, var_str1, date0, date1)
+        cov = self.time_covariance(var_str0, var_str1, date0, date1)
 
         if create_new_object:
             new_object = Altimetry()
             new_dataset = self.dataset[["longitude", "latitude", "time"]]
-            new_dataset["absolute_error"] = absolute_error
-            new_dataset["error"] = difference
-            new_dataset["mae"] = mean_absolute_error
-            new_dataset["rmse"] = root_mean_square_error
-            new_dataset["corr"] = correlation
-            new_dataset["cov"] = covariance
+            new_dataset["absolute_error"] = ae
+            new_dataset["error"] = diff
+            new_dataset["mae"] = mae
+            new_dataset["rmse"] = rmse
+            new_dataset["corr"] = corr
+            new_dataset["cov"] = cov
             new_object.dataset = new_dataset
             return new_object
         else:
-            self.dataset["absolute_error"] = absolute_error
-            self.dataset["error"] = difference
-            self.dataset["mae"] = mean_absolute_error
-            self.dataset["rmse"] = root_mean_square_error
-            self.dataset["corr"] = correlation
-            self.dataset["cov"] = covariance
+            self.dataset["absolute_error"] = ae
+            self.dataset["error"] = diff
+            self.dataset["mae"] = mae
+            self.dataset["rmse"] = rmse
+            self.dataset["corr"] = corr
+            self.dataset["cov"] = cov
diff --git a/coast/argos.py b/coast/argos.py
new file mode 100644
index 00000000..3f55110d
--- /dev/null
+++ b/coast/argos.py
@@ -0,0 +1,49 @@
+"""Argos class"""
+from .index import Indexed
+import numpy as np
+import xarray as xr
+import pandas as pd
+from .logging_util import get_slug, debug
+from typing import Union
+from pathlib import Path
+
+
+class Argos(Indexed):
+    """Class for reading Argos CSV formatted data files into an xarray object"""
+
+    def __init__(self, file_path: str = None, config: Union[Path, str] = None):
+        """Init Argos data object
+
+        Args:
+            file_path (str) : path of data fil
+            config (Path or str) : configuration file
+        """
+        debug(f"Creating a new {get_slug(self)}")
+        super().__init__(config)
+        self.read_data(file_path)
+        self.apply_config_mappings()
+        print(f"{get_slug(self)} initialised")
+
+    def read_data(self, file_path: str) -> None:
+        """Read the data file
+
+        Expected format and variable names are
+
+            DATIM,LAT,LON,SST,SST_F,PSST,PSST_F,PPRES,PPRES_F,BEN
+
+        xarray dataset to have dimension as time and coordinates as time, latitude and longitude
+
+        Args:
+            file_path (str) : path of data file
+        """
+        df = pd.read_csv(file_path, na_values=" ")
+
+        self.dataset = xr.Dataset()
+        var_names = list(df.columns)
+        for name in var_names:
+            if name == "DATIM":
+                time = np.array(pd.to_datetime(df[name]))
+                self.dataset[name] = xr.DataArray(list(time), dims=["time"])
+                self.dataset = self.dataset.set_coords(name)
+            else:
+                self.dataset[name] = xr.DataArray(list(np.array(df[name])), dims=["time"])
diff --git a/coast/config/config_parser.py b/coast/config/config_parser.py
index 9e6220ec..d6520b83 100644
--- a/coast/config/config_parser.py
+++ b/coast/config/config_parser.py
@@ -72,7 +72,6 @@ def _get_code_processing_object(json_content: dict) -> CodeProcessing:
         """
         dataset_json = json_content[ConfigKeys.CODE_PROCESSING]
         return CodeProcessing(
-            coord_variables=dataset_json[ConfigKeys.COO_VAR],
             delete_variables=dataset_json[ConfigKeys.DEL_VAR],
             not_grid_variables=dataset_json[ConfigKeys.NO_GR_VAR],
         )
@@ -89,7 +88,18 @@ def _get_datafile_object(json_content: dict, data_file_type: str) -> Union[Datas
         dataset_var = dataset_json[ConfigKeys.VAR_MAP]
         dataset_dim = dataset_json[ConfigKeys.DIM_MAP]
 
+        try:
+            dataset_keep_all_vars = dataset_json[ConfigKeys.KEEP_ALL_VARS]
+        except KeyError:
+            dataset_keep_all_vars = "False"
+
         if data_file_type is ConfigKeys.DATASET:
-            return Dataset(variable_map=dataset_var, dimension_map=dataset_dim)
+            dataset_coord_vars = dataset_json[ConfigKeys.COO_VAR]
+            return Dataset(
+                variable_map=dataset_var,
+                dimension_map=dataset_dim,
+                coord_var=dataset_coord_vars,
+                keep_all_vars=dataset_keep_all_vars,
+            )
         elif data_file_type is ConfigKeys.DOMAIN:
-            return Domain(variable_map=dataset_var, dimension_map=dataset_dim)
+            return Domain(variable_map=dataset_var, dimension_map=dataset_dim, keep_all_vars=False)
diff --git a/coast/config/config_structure.py b/coast/config/config_structure.py
index 8515bec3..75e987b9 100644
--- a/coast/config/config_structure.py
+++ b/coast/config/config_structure.py
@@ -27,6 +27,7 @@ class ConfigKeys:
     NO_GR_VAR = "not_grid_vars"
     COO_VAR = "coord_vars"
     DEL_VAR = "delete_vars"
+    KEEP_ALL_VARS = "keep_all_vars"
 
 
 @dataclass(frozen=True)
@@ -36,10 +37,13 @@ class DataFile:
     Args:
         variable_map (dict): dict containing mapping for variable names.
         dimension_map (dict): dict containing mapping for dimension names.
+        keep_all_vars (boolean): True if xarray is to retain all data file variables
+                                  otherwise False i.e keep only those in the json config file variable mappings.
     """
 
     variable_map: dict
     dimension_map: dict
+    keep_all_vars: bool = False
 
 
 @dataclass(frozen=True)
@@ -48,19 +52,22 @@ class CodeProcessing:
 
     Args:
         not_grid_variables (list): A list of variables not belonging to the grid.
-        coord_variables (list):  A list of coordinate variables.
         delete_variables (list):  A list of variables to drop from the dataset.
     """
 
     not_grid_variables: list
-    coord_variables: list
     delete_variables: list
 
 
 @dataclass(frozen=True)
 class Dataset(DataFile):
-    """Dataclass holding config attributes for Dataset datafiles. Extends DataFile."""
+    """Dataclass holding config attributes for Dataset datafiles. Extends DataFile.
 
+    Args:
+        coord_var (list): list of dataset coordinate variables to apply once dataset is loaded
+    """
+
+    coord_var: list = None
     pass
 
 
diff --git a/coast/general_utils.py b/coast/general_utils.py
index 3f7f1aa7..9a20f5ca 100755
--- a/coast/general_utils.py
+++ b/coast/general_utils.py
@@ -121,7 +121,7 @@ def compare_angles(a1, a2, degrees=True):
     return diff
 
 
-def cart2polar(x, y, degrees=True):
+def cartesian_to_polar(x, y, degrees=True):
     """
     # Conversion of cartesian to polar coordinate system
     # Output theta is in radians
@@ -133,7 +133,7 @@ def cart2polar(x, y, degrees=True):
     return r, theta
 
 
-def polar2cart(r, theta, degrees=True):
+def polar_to_cartesian(r, theta, degrees=True):
     """
     # Conversion of polar to cartesian coordinate system
     # Input theta must be in radians
@@ -283,7 +283,7 @@ def nearest_indices_2d(mod_lon, mod_lat, new_lon, new_lat, mask=None):
     return ind_x, ind_y
 
 
-def dataarray_time_slice(data_array, date0, date1):
+def data_array_time_slice(data_array, date0, date1):
     """Takes an xr.DataArray object and returns a new object with times
     sliced between dates date0 and date1. date0 and date1 may be a string or
     datetime type object."""
diff --git a/coast/glider.py b/coast/glider.py
new file mode 100644
index 00000000..f5ae63dc
--- /dev/null
+++ b/coast/glider.py
@@ -0,0 +1,37 @@
+"""Glider class"""
+from .index import Indexed
+import xarray as xr
+from .logging_util import get_slug, debug, info, warn, warning
+from typing import Union
+from pathlib import Path
+
+
+class Glider(Indexed):
+    """Glider class for reading in glider data (netcdf format) into an xarray object."""
+
+    def __init__(self, file_path=None, config: Union[Path, str] = None):
+        """Initialization and file reading.
+
+        Args:
+            file_path (str): path to data file
+            config (Union[Path, str]): path to json config file.
+        """
+        debug(f"Creating a new {get_slug(self)}")
+        super().__init__(config)
+
+        if file_path is None:
+            warn("Object created but no file or directory specified: \n" "{0}".format(str(self)), UserWarning)
+        else:
+            self.load_single(file_path)
+            self.apply_config_mappings()
+
+        print(f"{get_slug(self)} initialised")
+
+    def load_single(self, file_path, chunks: dict = None) -> None:
+        """Loads a single file into object's dataset variable.
+
+        Args:
+            file_path (str): path to data file
+            chunks (dict): chunks
+        """
+        self.dataset = xr.open_dataset(file_path, chunks=chunks)
diff --git a/coast/gridded.py b/coast/gridded.py
index adb4e1a0..1d68bd7e 100644
--- a/coast/gridded.py
+++ b/coast/gridded.py
@@ -132,7 +132,7 @@ def merge_domain_into_dataset(self, dataset_domain):
 
         # Reset & set specified coordinates
         self.dataset = self.dataset.reset_coords()
-        for var in self.config.code_processing.coord_variables:
+        for var in self.config.dataset.coord_var:
             try:
                 self.dataset = self.dataset.set_coords(var)
             except ValueError as err:
@@ -380,7 +380,7 @@ def interpolate_in_time(model_array, new_times, interp_method="nearest", extrapo
 
         return interpolated
 
-    def construct_density(self, EOS="EOS10"):
+    def construct_density(self, eos="EOS10"):
 
         """
             Constructs the in-situ density using the salinity, temperture and
@@ -396,7 +396,7 @@ def construct_density(self, EOS="EOS10"):
 
         Parameters
         ----------
-        EOS : equation of state, optional
+        eos : equation of state, optional
             DESCRIPTION. The default is 'EOS10'.
 
 
@@ -406,10 +406,10 @@ def construct_density(self, EOS="EOS10"):
         adds attribute NEMO.dataset.density
 
         """
-        debug(f'Constructing in-situ density for {get_slug(self)} with EOS "{EOS}"')
+        debug(f'Constructing in-situ density for {get_slug(self)} with EOS "{eos}"')
         try:
-            if EOS != "EOS10":
-                raise ValueError(str(self) + ": Density calculation for " + EOS + " not implemented.")
+            if eos != "EOS10":
+                raise ValueError(str(self) + ": Density calculation for " + eos + " not implemented.")
             if self.grid_ref != "t-grid":
                 raise ValueError(
                     str(self)
@@ -507,7 +507,7 @@ def copy_domain_vars_to_dataset(self, dataset_domain, grid_vars):
             except:  # FIXME Catch specific exception(s)
                 pass  # TODO Should we log something here?
 
-    def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=None, out_obj=None):
+    def differentiate(self, in_var_str, config_path=None, dim="z_dim", out_var_str=None, out_obj=None):
         """
         Derivatives are computed in x_dim, y_dim, z_dim (or i,j,k) directions
         wrt lambda, phi, or z coordinates (with scale factor in metres not degrees).
@@ -522,7 +522,7 @@ def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=Non
         1) d(grid_t)/dz --> grid_w
 
         Returns  an object (with the appropriate target grid_ref) containing
-        derivative (out_varstr) as xr.DataArray
+        derivative (out_var_str) as xr.DataArray
 
         This is hardwired to expect:
         1) depth_0 and e3_0 fields exist
@@ -541,30 +541,30 @@ def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=Non
 
         # For f(z)=-z. Compute df/dz = -1. Surface value is set to zero
         nemo_t.dataset['depth4D'],_ = xr.broadcast( nemo_t.dataset['depth_0'], nemo_t.dataset['temperature'] )
-        nemo_w_4 = nemo_t.differentiate( 'depth4D', dim='z_dim', out_varstr='dzdz' )
+        nemo_w_4 = nemo_t.differentiate( 'depth4D', dim='z_dim', out_var_str='dzdz' )
 
         Provide an existing target NEMO object and target variable name:
-        nemo_w_1 = nemo_t.differentiate( 'temperature', dim='z_dim', out_varstr='dTdz', out_obj=nemo_w_1 )
+        nemo_w_1 = nemo_t.differentiate( 'temperature', dim='z_dim', out_var_str='dTdz', out_obj=nemo_w_1 )
 
 
         Parameters
         ----------
-        in_varstr : str, name of variable to differentiate
+        in_var_str : str, name of variable to differentiate
         config_path : str, path to the w grid config file
         dim : str, dimension to operate over. E.g. {'z_dim', 'y_dim', 'x_dim', 't_dim'}
-        out_varstr : str, (optional) name of the target xr.DataArray
+        out_var_str : str, (optional) name of the target xr.DataArray
         out_obj : exiting NEMO obj to store xr.DataArray (optional)
 
         """
-        # import xarray as xr
+        import xarray as xr
 
         new_units = ""
 
-        # Check in_varstr exists in self.
-        if hasattr(self.dataset, in_varstr):
-            # self.dataset[in_varstr] exists
+        # Check in_var_str exists in self.
+        if hasattr(self.dataset, in_var_str):
+            # self.dataset[in_var_str] exists
 
-            var = self.dataset[in_varstr]  # for convenience
+            var = self.dataset[in_var_str]  # for convenience
 
             nt = var.sizes["t_dim"]
             nz = var.sizes["z_dim"]
@@ -586,9 +586,9 @@ def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=Non
                             "filename_domain={} is empty?".format(self.filename_domain),
                         )
 
-                # Check is out_varstr is defined, else create it
-                if out_varstr is None:
-                    out_varstr = in_varstr + "_dz"
+                # Check is out_var_str is defined, else create it
+                if out_var_str is None:
+                    out_var_str = in_var_str + "_dz"
 
                 # Create new DataArray with the same dimensions as the parent
                 # Crucially have a coordinate value that is appropriate to the target location.
@@ -598,12 +598,12 @@ def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=Non
                 diff = xr.concat([blank, var.diff(dim)], dim)
                 diff_ndim, e3w_ndim = xr.broadcast(diff, out_obj.dataset.e3_0.squeeze())
                 # Compute the derivative
-                out_obj.dataset[out_varstr] = -diff_ndim / e3w_ndim
+                out_obj.dataset[out_var_str] = -diff_ndim / e3w_ndim
 
                 # Assign attributes
                 new_units = var.units + "/" + out_obj.dataset.depth_0.units
                 # Convert to a xr.DataArray and return
-                out_obj.dataset[out_varstr].attrs = {"units": new_units, "standard_name": out_varstr}
+                out_obj.dataset[out_var_str].attrs = {"units": new_units, "standard_name": out_var_str}
 
                 # Return in object.
                 return out_obj
@@ -612,7 +612,7 @@ def differentiate(self, in_varstr, config_path=None, dim="z_dim", out_varstr=Non
                 warn("Not ready for that combination of grid ({}) and " "derivative ({})".format(self.grid_ref, dim))
                 return None
         else:
-            warn(f"{in_varstr} does not exist in {get_slug(self)} dataset")
+            warn(f"{in_var_str} does not exist in {get_slug(self)} dataset")
             return None
 
     def apply_doodson_x0_filter(self, var_str):
diff --git a/coast/index.py b/coast/index.py
new file mode 100644
index 00000000..5b853a9c
--- /dev/null
+++ b/coast/index.py
@@ -0,0 +1,64 @@
+"""Index class."""
+from dask import array
+from dask.distributed import Client
+from .config import config_parser
+from .logging_util import get_slug, debug, info, warn, warning
+from typing import Union
+from pathlib import Path
+from ast import literal_eval
+from .coast import Coast
+
+
+def setup_dask_client(workers: int = 2, threads: int = 2, memory_limit_per_worker: str = "2GB"):
+    Client(n_workers=workers, threads_per_worker=threads, memory_limit=memory_limit_per_worker)
+
+
+class Indexed(Coast):
+    def __init__(self, config: Union[Path, str] = None):
+        """Configuration init.
+
+        Args:
+            config (Union[Path, str]): path to json config file.
+        """
+        debug(f"Creating a new {get_slug(self)}")
+
+        self.dataset = None
+        self.chunks = None
+        self.var_mapping = None
+        self.dim_mapping = None
+        self.coord_vars = None
+        self.keep_all_vars = False
+
+        if config:
+            self.json_config = config_parser.ConfigParser(config)
+            self.chunks = self.json_config.config.chunks
+            self.dim_mapping = self.json_config.config.dataset.dimension_map
+            self.var_mapping = self.json_config.config.dataset.variable_map
+            self.coord_vars = self.json_config.config.dataset.coord_var
+            self.keep_all_vars = literal_eval(self.json_config.config.dataset.keep_all_vars)
+
+    def apply_config_mappings(self) -> None:
+        """Applies json configuration and mappings"""
+        # We iterate through each mapping one by one which enables us to catch those variables that do not exist
+
+        if self.dim_mapping is not None:
+            for k in self.dim_mapping:
+                try:
+                    self.dataset = self.dataset.rename_dims({k: self.dim_mapping[k]})
+                except ValueError as e:
+                    debug(f"Warning: {str(e)}")
+
+        if self.var_mapping is not None:
+            keep_vars = []
+            for k in self.var_mapping:
+                try:
+                    self.dataset = self.dataset.rename_vars({k: self.var_mapping[k]})
+                    keep_vars.append(self.var_mapping[k])
+                except ValueError as e:
+                    debug(f"Warning: {str(e)}")
+
+            if not self.keep_all_vars:
+                self.dataset = self.dataset[keep_vars]
+
+        if self.coord_vars is not None:
+            self.dataset = self.dataset.set_coords(self.coord_vars)
diff --git a/coast/internal_tide.py b/coast/internal_tide.py
index 65034769..4acdb5bb 100644
--- a/coast/internal_tide.py
+++ b/coast/internal_tide.py
@@ -114,12 +114,12 @@ def construct_pycnocline_vars(self, nemo_t: xr.Dataset, nemo_w: xr.Dataset, stra
         debug(f"Constructing pycnocline variables for {get_slug(self)}")
         # Construct in-situ density if not already done
         if not hasattr(nemo_t.dataset, "density"):
-            nemo_t.construct_density(EOS="EOS10")
+            nemo_t.construct_density(eos="EOS10")
 
         # Construct stratification if not already done. t-pts --> w-pts
         if not hasattr(nemo_w.dataset, "rho_dz"):
             nemo_w = nemo_t.differentiate(
-                "density", dim="z_dim", out_varstr="rho_dz", out_obj=nemo_w
+                "density", dim="z_dim", out_var_str="rho_dz", out_obj=nemo_w
             )  # TODO These kwargs don't appear to exist
 
         # Define the spatial dimensional size and check the dataset and domain arrays are the same size in
diff --git a/coast/lagrangian.py b/coast/lagrangian.py
new file mode 100644
index 00000000..0a98550d
--- /dev/null
+++ b/coast/lagrangian.py
@@ -0,0 +1,10 @@
+"""Lagrangian class"""
+from .index import Indexed
+
+
+class Lagrangian(Indexed):
+    """Parent class for subclasses OCEANPARCELS ...
+    Common methods ....
+    """
+
+    pass
diff --git a/coast/oceanparcels.py b/coast/oceanparcels.py
new file mode 100644
index 00000000..b81cf5a4
--- /dev/null
+++ b/coast/oceanparcels.py
@@ -0,0 +1,36 @@
+"""Oceanparcels class for reading ocean parcels data."""
+from .lagrangian import Lagrangian
+import xarray as xr
+from .logging_util import get_slug, debug, info, warn, warning
+from typing import Union
+from pathlib import Path
+
+
+class Oceanparcels(Lagrangian):
+    """Reading ocean parcels data (netcdf format) into an xarray object."""
+
+    def __init__(self, file_path: str = None, config: Union[Path, str] = None):
+        """Initialization and file reading.
+
+        Args:
+           file_path (str): path to data file
+           config (Union[Path, str]): path to json config file.
+        """
+        debug(f"Creating a new {get_slug(self)}")
+        super().__init__(config)
+
+        if file_path is None:
+            warn("Object created but no file or directory specified: \n" "{0}".format(str(self)), UserWarning)
+        else:
+            self.load_single(file_path)
+            self.apply_config_mappings()
+
+        print(f"{get_slug(self)} initialised")
+
+    def load_single(self, file_path) -> None:
+        """Loads a single file into object's dataset variable.
+
+        Args:
+            file_path (str): path to data file
+        """
+        self.dataset = xr.open_dataset(file_path, chunks=self.chunks)
diff --git a/coast/profile.py b/coast/profile.py
index 9b61e79e..8688ef3e 100755
--- a/coast/profile.py
+++ b/coast/profile.py
@@ -1,13 +1,17 @@
+"""Profile Class"""
+from .index import Indexed
 import numpy as np
 import xarray as xr
 from . import general_utils, plot_util
-from .coast import Coast
 import matplotlib.pyplot as plt
 import glob
 import datetime
+from .logging_util import get_slug, debug, info, warn, warning
+from typing import Union
+from pathlib import Path
 
 
-class Profile(Coast):
+class Profile(Indexed):
     """
     OBSERVATION type class for storing data from a CTD Profile (or similar
     down and up observations). The structure of the class is based on data from
@@ -20,18 +24,32 @@ class Profile(Coast):
                      files.
     """
 
-    ##############################################################################
-    ###                ~ Initialisation and File Reading ~                     ###
-    ##############################################################################
+    def __init__(self, file_path: str = None, multiple=False, config: Union[Path, str] = None):
+        """Initialization and file reading.
 
-    def __init__(self):
-        self.dataset = None
-        return
-
-    def read_en4(self, fn_en4, multiple=False, chunks={}):
+        Args:
+            file_path (str): path to data file
+            multiple (boolean): True if reading multiple files otherwise False
+            config (Union[Path, str]): path to json config file.
         """
-        Reads a single or multiple EN4 netCDF files into the COAsT profile
-        data structure.
+        debug(f"Creating a new {get_slug(self)}")
+        super().__init__(config)
+
+        if file_path is None:
+            warn("Object created but no file or directory specified: \n" "{0}".format(str(self)), UserWarning)
+        else:
+            self.read_en4(file_path, self.chunks, multiple)
+            self.apply_config_mappings()
+
+        print(f"{get_slug(self)} initialised")
+
+    def read_en4(self, fn_en4, chunks: dict = {}, multiple=False) -> None:
+        """Reads a single or multiple EN4 netCDF files into the COAsT profile data structure.
+
+        Args:
+            fn_en4 (str): path to data file
+            chunks (dict): chunks
+            multiple (boolean): True if reading multiple files otherwise False
         """
         if not multiple:
             self.dataset = xr.open_dataset(fn_en4, chunks=chunks)
@@ -62,36 +80,7 @@ def read_en4(self, fn_en4, multiple=False, chunks={}):
                 else:
                     self.dataset = xr.concat((self.dataset, data_tmp), dim="N_PROF")
 
-        rename_vars = {
-            "LATITUDE": "latitude",
-            "LONGITUDE": "longitude",
-            "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",
-            "QC_FLAGS_PROFILES": "qc_flags_profiles",
-            "QC_FLAGS_LEVELS": "qc_flags_levels",
-        }
-        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           ~                     ###
-    ##############################################################################
+    """======================= Manipulate ======================="""
 
     def subset_indices_lonlat_box(self, lonbounds, latbounds):
         """Generates array indices for data which lies in a given lon/lat box.
@@ -108,9 +97,7 @@ def subset_indices_lonlat_box(self, lonbounds, latbounds):
         )
         return ind
 
-    ##############################################################################
-    ###                ~            Plotting             ~                     ###
-    ##############################################################################
+    """======================= Plotting ======================="""
 
     def plot_profile(self, var: str, profile_indices=None):
 
@@ -129,7 +116,7 @@ def plot_profile(self, var: str, profile_indices=None):
         plt.xlabel(var + "(" + self.dataset[var].units + ")")
         plt.ylabel("Depth (" + self.dataset.depth.units + ")")
         plt.grid()
-        return fig, ax  # TODO ax might not be defined
+        return fig, ax
 
     def plot_map(self, profile_indices=None, var_str=None, depth_index=None):
 
@@ -144,7 +131,6 @@ def plot_map(self, profile_indices=None, var_str=None, depth_index=None):
             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="potential_temperature", var_s="practical_salinity"):
@@ -157,7 +143,4 @@ def plot_ts_diagram(self, profile_index, var_t="potential_temperature", var_s="p
 
         return fig, ax
 
-
-##############################################################################
-###                ~        Model Comparison         ~                     ###
-##############################################################################
+    """======================= Model Comparison ======================="""
diff --git a/coast/tidegauge.py b/coast/tidegauge.py
new file mode 100755
index 00000000..ef527fe4
--- /dev/null
+++ b/coast/tidegauge.py
@@ -0,0 +1,1290 @@
+"""Tide Gauge class"""
+from .timeseries import Timeseries
+import numpy as np
+import xarray as xr
+import matplotlib.pyplot as plt
+import pandas as pd
+import glob
+import re
+import pytz
+import sklearn.metrics as metrics
+from . import general_utils, plot_util, crps_util, stats_util
+from .logging_util import get_slug, debug, error, info
+from typing import Union
+from pathlib import Path
+
+
+class Tidegauge(Timeseries):
+    """
+    An object for reading, storing and manipulating tide gauge data.
+    Functionality available for reading and organisation of GESLA files.
+    (Source: https://www.gesla.org/).  However, any fixed time series data can
+    be used if in the correct format.
+
+    The data format used for this object is as follows:
+
+    *Data Format Overview*
+
+        1. Data for a single tide gauge is stored in an xarray Dataset object.
+           This can be accessed using TIDEGAUGE.dataset.
+        2. The dataset has a single dimension: time.
+        3. Latitude/Longitude and other single values parameters are stored as
+           attributes or single float variables.
+        4. Time is a coordinate variable and time dimension.
+        5. Data variables are stored along the time dimension.
+        6. The attributes: site_name, latitude, longitude are expected. If they
+            are missing functionality may be reduced.
+
+
+    *Methods Overview*
+
+        *Initialisation and File Reading*
+        -> __init__: Can be initialised with a GESLA file or empty.
+        -> obs_operator: Interpolates model data to time series locations
+           and times (not yet implemented).
+        -> read_gesla_to_xarray_v3: Reads a format version 3.0
+           GESLA file to an xarray Dataset.
+        -> read_gesla_header_v3: Reads the header of a version 3
+           GESLA file.
+        -> read_gesla_data_v3: Reads data from a version 3 GESLA
+           file.
+        -> create_multiple_tidegauge: Creates multiple tide gauge objects
+           objects from a list of filenames or directory and returns them
+           in a list.
+
+        *Plotting*
+        -> plot_on_map: Plots location of TIDEGAUGE object on map.
+        -> plot_timeseries: Plots a specified time series.
+
+        *Model Comparison*
+        -> obs_operator(): For interpolating model data to this object.
+        -> cprs(): Calculates the CRPS between a model and obs variable.
+        -> difference(): Differences two specified variables
+        -> absolute_error(): Absolute difference, two variables
+        -> mean_absolute_error(): MAE between two variables
+        -> root_mean_square_error(): RMSE between two variables
+        -> time_mean(): Mean of a variable in time
+        -> time_std(): St. Dev of a variable in time
+        -> time_correlation(): Correlation between two variables
+        -> time_covariance(): Covariance between two variables
+        -> basic_stats(): Calculates multiple of the above metrics.
+
+        *Analysis*
+        -> resample_mean(): For resampling data in time using averaging
+        -> apply_doodson_xo_filter(): Remove tidal signal using Doodson XO
+        -> find_high_and_low_water(): Find maxima and minima of time series
+    """
+
+    def __init__(self, file_path: str = None, date_start=None, date_end=None, config: Union[Path, str] = None):
+        """
+        Initialise TIDEGAUGE object either as empty (no arguments) or by
+        reading GESLA data from a directory between two datetime objects.
+
+        Example usage:
+        --------------
+        # Read tide gauge data for data in January 1990
+        date0 = datetime.datetime(1990,1,1)
+        date1 = datetime.datetime(1990,2,1)
+        tg = coast.TIDEGAUGE(<'path_to_file'>, date0, date1)
+
+        # Access the data
+        tg.dataset
+
+        Parameters
+        ----------
+        file_path (list of str) : Filename to read from directory.
+        date_start (datetime) : Start date for data read. Optional
+        date_end (datetime) : end date for data read. Optional
+        config (Path or str) : configuration file
+
+        Returns
+        -------
+        Self
+        """
+        debug(f"Creating a new ..... {get_slug(self)}")
+        super().__init__(config)
+
+        # If file list is supplied, read files from directory
+        if file_path is None:
+            self.dataset = None
+        else:
+            self.dataset = self.read_gesla_to_xarray_v3(file_path, date_start, date_end)
+            self.apply_config_mappings()
+
+        print(f"{get_slug(self)} initialised")
+
+    ############ tide gauge methods ##############################################
+    @classmethod
+    def read_gesla_to_xarray_v3(cls, fn_gesla, date_start=None, date_end=None):
+        """
+        For reading from a single GESLA2 (Format version 3.0) file into an
+        xarray dataset. Formatting according to Woodworth et al. (2017).
+        Website: https://www.gesla.org/
+        If no data lies between the specified dates, a dataset is still created
+        containing information on the tide gauge, but the time dimension will
+        be empty.
+        Parameters
+        ----------
+        fn_gesla (str) : path to gesla tide gauge file
+        date_start (datetime) : start date for returning data
+        date_end (datetime) : end date for returning data
+
+        Returns
+        -------
+        xarray.Dataset object.
+        """
+        debug(f'Reading "{fn_gesla}" as a GESLA file with {get_slug(cls)}')  # TODO Maybe include start/end dates
+        try:
+            header_dict = cls.read_gesla_header_v3(fn_gesla)
+            dataset = cls.read_gesla_data_v3(fn_gesla, date_start, date_end)
+        except:
+            raise Exception("Problem reading GESLA file: " + fn_gesla)
+        # Attributes
+        dataset["longitude"] = header_dict["longitude"]
+        dataset["latitude"] = header_dict["latitude"]
+        del header_dict["longitude"]
+        del header_dict["latitude"]
+
+        dataset.attrs = header_dict
+
+        return dataset
+
+    @staticmethod
+    def read_gesla_header_v3(fn_gesla):
+        """
+        Reads header from a GESLA file (format version 3.0).
+
+        Parameters
+        ----------
+        fn_gesla (str) : path to gesla tide gauge file
+
+        Returns
+        -------
+        dictionary of attributes
+        """
+        debug(f'Reading GESLA header from "{fn_gesla}"')
+        fid = open(fn_gesla)
+
+        # Read lines one by one (hopefully formatting is consistent)
+        fid.readline()  # Skip first line
+        # Geographical stuff
+        site_name = fid.readline().split()[3:]
+        site_name = "_".join(site_name)
+        country = fid.readline().split()[2:]
+        country = "_".join(country)
+        contributor = fid.readline().split()[2:]
+        contributor = "_".join(contributor)
+        # Coordinates
+        latitude = float(fid.readline().split()[2])
+        longitude = float(fid.readline().split()[2])
+        coordinate_system = fid.readline().split()[3]
+        # Dates
+        start_date = fid.readline().split()[3:5]
+        start_date = " ".join(start_date)
+        start_date = pd.to_datetime(start_date)
+        end_date = fid.readline().split()[3:5]
+        end_date = " ".join(end_date)
+        end_date = pd.to_datetime(end_date)
+        time_zone_hours = float(fid.readline().split()[4])
+        # Other
+        fid.readline()  # Datum
+        fid.readline()  # Instrument
+        precision = float(fid.readline().split()[2])
+        null_value = float(fid.readline().split()[3])
+
+        debug(f'Read done, close file "{fn_gesla}"')
+        fid.close()
+        # Put all header info into an attributes dictionary
+        header_dict = {
+            "site_name": site_name,
+            "country": country,
+            "contributor": contributor,
+            "latitude": latitude,
+            "longitude": longitude,
+            "coordinate_system": coordinate_system,
+            "original_start_date": start_date,
+            "original_end_date": end_date,
+            "time_zone_hours": time_zone_hours,
+            "precision": precision,
+            "null_value": null_value,
+        }
+        return header_dict
+
+    @staticmethod
+    def read_gesla_data_v3(fn_gesla, date_start=None, date_end=None, header_length: int = 32):
+        """
+        Reads observation data from a GESLA file (format version 3.0).
+
+        Parameters
+        ----------
+        fn_gesla (str) : path to gesla tide gauge file
+        date_start (datetime) : start date for returning data
+        date_end (datetime) : end date for returning data
+        header_length (int) : number of lines in header (to skip when reading)
+
+        Returns
+        -------
+        xarray.Dataset containing times, sealevel and quality control flags
+        """
+        # Initialise empty dataset and lists
+        debug(f'Reading GESLA data from "{fn_gesla}"')
+        dataset = xr.Dataset()
+        time = []
+        sea_level = []
+        qc_flags = []
+        # Open file and loop until EOF
+        with open(fn_gesla) as file:
+            line_count = 1
+            for line in file:
+                # Read all data. Date boundaries are set later.
+                if line_count > header_length:
+                    working_line = line.split()
+                    if working_line[0] != "#":
+                        time.append(working_line[0] + " " + working_line[1])
+                        sea_level.append(float(working_line[2]))
+                        qc_flags.append(int(working_line[3]))
+
+                line_count = line_count + 1
+            debug(f'Read done, close file "{fn_gesla}"')
+
+        # Convert time list to datetimes using pandas
+        time = np.array(pd.to_datetime(time))
+
+        # Return only values between stated dates
+        start_index = 0
+        end_index = len(time)
+        if date_start is not None:
+            date_start = np.datetime64(date_start)
+            start_index = np.argmax(time >= date_start)
+        if date_end is not None:
+            date_end = np.datetime64(date_end)
+            end_index = np.argmax(time > date_end)
+        time = time[start_index:end_index]
+        sea_level = sea_level[start_index:end_index]
+        qc_flags = qc_flags[start_index:end_index]
+
+        # Set null values to nan
+        sea_level = np.array(sea_level)
+        qc_flags = np.array(qc_flags)
+        sea_level[qc_flags == 5] = np.nan
+
+        # Assign arrays to Dataset
+        dataset["sea_level"] = xr.DataArray(sea_level, dims=["time"])
+        dataset["qc_flags"] = xr.DataArray(qc_flags, dims=["time"])
+        dataset = dataset.assign_coords(time=("time", time))
+
+        # Assign local dataset to object-scope dataset
+        return dataset
+
+    @classmethod
+    def create_multiple_tidegauge(cls, file_list, date_start=None, date_end=None, config: Union[Path, str] = None):
+        """
+        Reads multiple GESLA tide gauge files from file_list (can include
+        wildcards) and return them in a list. date_start and date_end should
+        be datetime like objects. For a lot of files/data, this may take a
+        while.
+
+        Example usage:
+        --------------
+            # Read all data in directory in January 1990
+            date0 = datetime.datetime(1990,1,1)
+            date1 = datetime.datetime(1990,2,1)
+            tg = coast.TIDEGAUGE('gesla_directory/*', date0, date1)
+        Returns
+        -------
+        List of TIDEGAUGE objects.
+        """
+        # If single string is given then put into a single element list
+        if type(file_list) is str:
+            file_list = [file_list]
+
+        # Check file_list for wildcards and make list of files to read
+        file_to_read = []
+        for file in file_list:
+            if "*" in file:
+                wildcard_list = glob.glob(file)
+                file_to_read = file_to_read + wildcard_list
+            else:
+                file_to_read.append(file)
+
+        # Loop over files to read and read them into datasets
+        tidegauge_list = []
+        for file in file_to_read:
+            try:
+                new_object = Tidegauge(file, date_start, date_end, config)
+                tidegauge_list.append(new_object)
+            except FileNotFoundError as fne:
+                Warning(str(fne))
+
+        return tidegauge_list
+
+    ############ tide table methods (HLW) #########################################
+    @classmethod
+    def read_hlw_to_xarray(cls, fn_hlw, date_start=None, date_end=None):
+        """
+        For reading from a file of tidetable High and Low Waters (HLW) data into an
+        xarray dataset. File contains high water and low water heights and times
+
+        If no data lies between the specified dates, a dataset is still created
+        containing information on the tide gauge, but the time dimension will
+        be empty.
+
+        The data takes the form:
+        LIVERPOOL (GLADSTONE DOCK)    TZ: UT(GMT)/BST     Units: METRES    Datum: Chart Datum
+        01/10/2020  06:29    1.65
+        01/10/2020  11:54    9.01
+        01/10/2020  18:36    1.87
+        ...
+
+        Parameters
+        ----------
+        fn_hlw (str) : path to tabulated High Low Water file
+        date_start (datetime) : start date for returning data
+        date_end (datetime) : end date for returning data
+
+        Returns
+        -------
+        xarray.Dataset object.
+        """
+        debug(f'Reading "{fn_hlw}" as a HLW file with {get_slug(cls)}')  # TODO Maybe include start/end dates
+        try:
+            header_dict = cls.read_HLW_header(fn_hlw)
+            dataset = cls.read_HLW_data(fn_hlw, header_dict, date_start, date_end)
+            if header_dict["field"] == "TZ:UT(GMT)/BST":
+                debug("Read in as BST, stored as UTC")
+            elif header_dict["field"] == "TZ:GMTonly":
+                debug("Read and store as GMT/UTC")
+            else:
+                debug("Not expecting that timezone")
+
+        except:
+            raise Exception("Problem reading HLW file: " + fn_hlw)
+
+        dataset.attrs = header_dict
+
+        return dataset
+
+    @staticmethod
+    def read_HLW_header(filnam):
+        """
+        Reads header from a HWL file.
+
+        The data takes the form:
+        LIVERPOOL (GLADSTONE DOCK)    TZ: UT(GMT)/BST     Units: METRES    Datum: Chart Datum
+        01/10/2020  06:29    1.65
+        01/10/2020  11:54    9.01
+        01/10/2020  18:36    1.87
+        ...
+
+        Parameters
+        ----------
+        filnam (str) : path to file
+
+        Returns
+        -------
+        dictionary of attributes
+        """
+        debug(f'Reading HLW header from "{filnam}" ')
+        fid = open(filnam)
+
+        # Read lines one by one (hopefully formatting is consistent)
+        header = re.split(r"\s{2,}", fid.readline())
+        site_name = header[0]
+        site_name = site_name.replace(" ", "")
+
+        field = header[1]
+        field = field.replace(" ", "")
+
+        units = header[2]
+        units = units.replace(" ", "")
+
+        datum = header[3]
+        datum = datum.replace(" ", "")
+
+        debug(f'Read done, close file "{filnam}"')
+        fid.close()
+        # Put all header info into an attributes dictionary
+        header_dict = {"site_name": site_name, "field": field, "units": units, "datum": datum}
+        return header_dict
+
+    @staticmethod
+    def read_HLW_data(filnam, header_dict, date_start=None, date_end=None, header_length: int = 1):
+        """
+        Reads HLW data from a tidetable file.
+
+        Parameters
+        ----------
+        filnam (str) : path to HLW tide gauge file
+        date_start (np.datetime64) : start date for returning data.
+        date_end (np.datetime64) : end date for returning data.
+        header_length (int) : number of lines in header (to skip when reading)
+
+        Returns
+        -------
+        xarray.Dataset containing times, High and Low water values
+        """
+        import datetime
+
+        # Initialise empty dataset and lists
+        debug(f'Reading HLW data from "{filnam}"')
+        dataset = xr.Dataset()
+        time = []
+        sea_level = []
+
+        if header_dict["field"] == "TZ:UT(GMT)/BST":
+            localtime_flag = True
+        else:
+            localtime_flag = False
+
+        # Open file and loop until EOF
+        with open(filnam) as file:
+            line_count = 1
+            for line in file:
+                # Read all data. Date boundaries are set later.
+                if line_count > header_length:
+                    working_line = line.split()
+                    if working_line[0] != "#":
+                        time_str = working_line[0] + " " + working_line[1]
+                        # Read time as datetime.datetime because it can handle local timezone easily AND the unusual date format
+                        datetime_obj = datetime.datetime.strptime(time_str, "%d/%m/%Y %H:%M")
+                        if localtime_flag == True:
+                            time.append(np.datetime64(datetime_obj.astimezone()))
+                        else:
+                            time.append(np.datetime64(datetime_obj))
+                        sea_level.append(float(working_line[2]))
+                line_count = line_count + 1
+            debug(f'Read done, close file "{filnam}"')
+
+        # Return only values between stated dates
+        start_index = 0
+        end_index = len(time)
+        if date_start is not None:
+            # start_index = general_utils.nearest_datetime_ind(time, date_start)
+            date_start = np.datetime64(date_start)
+            start_index = np.argmax(time >= date_start)
+            debug(f"date_start: {date_start}. start_index: {start_index}")
+        if date_end is not None:
+            # end_index = general_utils.nearest_datetime_ind(time, date_end)
+            date_end = np.datetime64(date_end)
+            end_index = np.argmax(time > date_end)
+            debug(f"date_end: {date_end}. end_index: {end_index}")
+        time = time[start_index:end_index]
+        sea_level = sea_level[start_index:end_index]
+        debug(f"sea_level: {sea_level}")
+        # Assign arrays to Dataset
+        dataset["sea_level"] = xr.DataArray(sea_level, dims=["time"])
+        dataset = dataset.assign_coords(time=("time", time))
+        # Assign local dataset to object-scope dataset
+        return dataset
+
+    def show(self, timezone: str = None):
+        """
+        Print out the values in the xarray
+        Displays with specified timezone
+        """
+        # debug(" Saltney pred", np.datetime_as_string(Saltney_time_pred[i], unit='m', timezone=pytz.timezone('Europe/London')),". Height: {:.2f} m".format( HT.values[i] ))
+        if timezone == None:
+            for i in range(len(self.dataset.sea_level)):
+                #               debug('time:', self.dataset.time[i].values,
+                debug(
+                    "time (UTC):",
+                    general_utils.dayoweek(self.dataset.time[i].values),
+                    np.datetime_as_string(self.dataset.time[i], unit="m"),
+                    "height:",
+                    self.dataset.sea_level[i].values,
+                    "m",
+                )
+        else:  # display timezone aware times
+            for i in range(len(self.dataset.sea_level)):
+                #               debug('time:', self.dataset.time[i].values,
+                debug(
+                    "time (" + timezone + "):",
+                    general_utils.day_of_week(self.dataset.time[i].values),
+                    np.datetime_as_string(self.dataset.time[i], unit="m", timezone=pytz.timezone(timezone)),
+                    "height:",
+                    self.dataset.sea_level[i].values,
+                    "m",
+                )
+
+    def get_tide_table_times(
+        self,
+        time_guess: np.datetime64 = None,
+        time_var: str = "time",
+        measure_var: str = "sea_level",
+        method: str = "window",
+        winsize=None,
+    ):
+        """
+        Get tide times and heights from tide table.
+        input:
+        time_guess : np.datetime64 or datetime
+                assumes utc
+        time_var : name of time variable [default: 'time']
+        measure_var : name of sea_level variable [default: 'sea_level']
+
+        method =
+            window:  +/- hours window size, winsize, (int) return values in that window
+                uses additional variable winsize (int) [default 2hrs]
+            nearest_1: return only the nearest event, if in winsize [default:None]
+            nearest_2: return nearest event in future and the nearest in the past (i.e. high and a low), if in winsize [default:None]
+            nearest_HW: return nearest High Water event (computed as the max of `nearest_2`), if in winsize [default:None]
+
+        returns: xr.DataArray( measure_var, coords=time_var)
+            E.g. sea_level (m), time (utc)
+            If value is not found, it returns a NaN with time value as the
+            guess value.
+
+        """
+        # Ensure the date objects are datetime
+        if type(time_guess) is not np.datetime64:
+            debug("Convert date to np.datetime64")
+            time_guess = np.datetime64(time_guess)
+
+        if time_guess == None:
+            debug("Use today's date")
+            time_guess = np.datetime64("now")
+
+        if method == "window":
+            if winsize == None:
+                winsize = 2
+            # initialise start_index and end_index
+            start_index = 0
+            end_index = len(self.dataset[time_var])
+
+            date_start = time_guess - np.timedelta64(winsize, "h")
+            start_index = np.argmax(self.dataset[time_var].values >= date_start)
+
+            date_end = time_guess + np.timedelta64(winsize, "h")
+            end_index = np.argmax(self.dataset[time_var].values > date_end)
+
+            sea_level = self.dataset[measure_var][start_index:end_index]
+
+            return sea_level
+
+        elif method == "nearest_1":
+            dt = np.abs(self.dataset[time_var] - time_guess)
+            index = np.argsort(dt).values
+            if winsize is not None:  # if search window trucation exists
+                if np.timedelta64(dt[index[0]].values, "m").astype("int") <= 60 * winsize:  # compare in minutes
+                    debug(f"dt:{np.timedelta64(dt[index[0]].values, 'm').astype('int')}")
+                    debug(f"winsize:{winsize}")
+                    return self.dataset[measure_var][index[0]]
+                else:
+                    # return a NaN in an xr.Dataset
+                    # The rather odd trailing zero is to remove the array layer
+                    # on both time and measurement, and to match the other
+                    # alternative for a return object
+                    return xr.DataArray([np.NaN], dims=(time_var), coords={time_var: [time_guess]})[0]
+            else:  # give the closest without window search truncation
+                return self.dataset[measure_var][index[0]]
+
+        elif method == "nearest_2":
+            index = np.argsort(np.abs(self.dataset[time_var] - time_guess)).values
+            nearest_2 = self.dataset[measure_var][index[0 : 1 + 1]]  # , self.dataset.time[index[0:1+1]]
+            return nearest_2
+
+        elif method == "nearest_HW":
+            index = np.argsort(np.abs(self.dataset[time_var] - time_guess)).values
+            # return self.dataset.sea_level[ index[np.argmax( self.dataset.sea_level[index[0:1+1]]] )] #, self.dataset.time[index[0:1+1]]
+            nearest_2 = self.dataset[measure_var][index[0 : 1 + 1]]  # , self.dataset.time[index[0:1+1]]
+            return nearest_2[nearest_2.argmax()]
+
+        else:
+            debug("Not expecting that option / method")
+
+    ############ environment.data.gov.uk gauge methods ###########################
+    @classmethod
+    def read_EA_API_to_xarray(
+        cls, ndays: int = 5, date_start: np.datetime64 = None, date_end: np.datetime64 = None, stationId="E70124"
+    ):
+        """
+        load gauge data via environment.data.gov.uk EA API
+        Either loads last ndays, or from date_start:date_end
+
+        API Source:
+        https://environment.data.gov.uk/flood-monitoring/doc/reference
+
+        Details of available tidal stations are recovered with:
+        https://environment.data.gov.uk/flood-monitoring/id/stations?type=TideGauge
+        Recover the "stationReference" for the gauge of interest and pass as
+        stationId:str. The default stationId="E70124" is Liverpool.
+
+        INPUTS:
+            ndays : int. Extact the last ndays from now.
+            date_start : datetime. UTC format string "yyyy-MM-dd" E.g 2020-01-05
+            date_end : datetime
+            stationId : int. Station id. Also referred to as stationReference in
+             EA API. Default value is for Liverpool.
+        OUTPUT:
+            sea_level, time : xr.Dataset
+        """
+        import requests, json
+
+        cls.ndays = ndays
+        cls.date_start = date_start
+        cls.date_end = date_end
+        cls.stationId = stationId  # EA id: stationReference
+
+        # %% Obtain and process header information
+        info("load station info")
+        url = "https://environment.data.gov.uk/flood-monitoring/id/stations/" + cls.stationId + ".json"
+        try:
+            request_raw = requests.get(url)
+            header_dict = json.loads(request_raw.content)
+        except ValueError:
+            debug(f"Failed request for station {cls.stationId}")
+            return
+
+        try:
+            header_dict["site_name"] = header_dict["items"]["label"]
+            header_dict["latitude"] = header_dict["items"]["lat"]
+            header_dict["longitude"] = header_dict["items"]["long"]
+        except:
+            info(f"possible missing some header info: site_name,latitude,longitude")
+        try:
+            # Define url call with parameter from station info
+            htmlcall_stationId = header_dict["items"]["measures"]["@id"] + "/readings?"
+        except:
+            debug(f"problem defining the parameter to read")
+
+        # %% Construct API request for data recovery
+        info("load station data")
+        if (cls.date_start == None) & (cls.date_end == None):
+            info(f"GETting ndays= {cls.ndays} of data")
+            url = (
+                htmlcall_stationId
+                + "since="
+                + (np.datetime64("now") - np.timedelta64(ndays, "D")).item().strftime("%Y-%m-%dT%H:%M:%SZ")
+            )
+            debug(f"url request: {url}")
+        else:
+            # Check date_start and date_end are timetime objects
+            if (type(cls.date_start) is np.datetime64) & (type(cls.date_end) is np.datetime64):
+                info(f"GETting data from {cls.date_start} to {cls.date_end}")
+                startdate = cls.date_start.item().strftime("%Y-%m-%d")
+                enddate = cls.date_end.item().strftime("%Y-%m-%d")
+                url = htmlcall_stationId + "startdate=" + startdate + "&enddate=" + enddate
+                debug(f"url request: {url}")
+
+            else:
+                debug("Expecting date_start and date_end as datetime objects")
+
+        # %% Get the data
+        try:
+            request_raw = requests.get(url)
+            request = json.loads(request_raw.content)
+            debug(f"EA API request: {request_raw.text}")
+        except ValueError:
+            debug(f"Failed request: {request_raw}")
+            return
+
+        # %% Process timeseries data
+        dataset = xr.Dataset()
+        time = []
+        sea_level = []
+        nvals = len(request["items"])
+        time = np.array([np.datetime64(request["items"][i]["dateTime"]) for i in range(nvals)])
+        sea_level = np.array([request["items"][i]["value"] for i in range(nvals)])
+
+        # %% Assign arrays to Dataset
+        dataset["sea_level"] = xr.DataArray(sea_level, dims=["time"])
+        dataset = dataset.assign_coords(time=("time", time))
+        dataset.attrs = header_dict
+        debug(f"EA API request headers: {header_dict}")
+        # debug(f"EA API request 1st time: {time[0]} and value: {sea_level[0]}")
+
+        # Assign local dataset to object-scope dataset
+        return dataset
+
+    ############ BODC tide gauge methods ##############################################
+    @classmethod
+    def read_bodc_to_xarray(cls, fn_bodc, date_start=None, date_end=None):
+        """
+        For reading from a single BODC (processed) file into an
+        xarray dataset.
+        If no data lies between the specified dates, a dataset is still created
+        containing information on the tide gauge, but the time dimension will
+        be empty.
+
+        Data name: UK Tide Gauge Network, processed data.
+        Source: https://www.bodc.ac.uk/
+        See data notes from source for description of QC flags.
+
+        The data takes the form:
+            Port:              P234
+            Site:              Liverpool, Gladstone Dock
+            Latitude:          53.44969
+            Longitude:         -3.01800
+            Start Date:        01AUG2020-00.00.00
+            End Date:          31AUG2020-23.45.00
+            Contributor:       National Oceanography Centre, Liverpool
+            Datum information: The data refer to Admiralty Chart Datum (ACD)
+            Parameter code:    ASLVBG02 = Surface elevation (unspecified datum) of the water body by bubbler tide gauge (second sensor)
+              Cycle    Date      Time    ASLVBG02   Residual
+             Number yyyy mm dd hh mi ssf         f          f
+                 1) 2020/08/01 00:00:00     5.354M     0.265M
+                 2) 2020/08/01 00:15:00     5.016M     0.243M
+                 3) 2020/08/01 00:30:00     4.704M     0.241M
+                 4) 2020/08/01 00:45:00     4.418M     0.255M
+                 5) 2020/08/01 01:00:00     4.133      0.257
+                 ...
+
+        Parameters
+        ----------
+        fn_bodc (str) : path to bodc tide gauge file
+        date_start (datetime) : start date for returning data
+        date_end (datetime) : end date for returning data
+
+        Returns
+        -------
+        xarray.Dataset object.
+        """
+        debug(f'Reading "{fn_bodc}" as a BODC file with {get_slug(cls)}')  # TODO Maybe include start/end dates
+        try:
+            header_dict = cls.read_bodc_header(fn_bodc)
+            dataset = cls.read_bodc_data(fn_bodc, date_start, date_end)
+        except:
+            raise Exception("Problem reading BODC file: " + fn_bodc)
+        # Attributes
+        dataset["longitude"] = header_dict["longitude"]
+        dataset["latitude"] = header_dict["latitude"]
+        del header_dict["longitude"]
+        del header_dict["latitude"]
+
+        dataset.attrs = header_dict
+
+        return dataset
+
+    @staticmethod
+    def read_bodc_header(fn_bodc):
+        """
+        Reads header from a BODC file (format version 3.0).
+
+        Parameters
+        ----------
+        fn_bodc (str) : path to bodc tide gauge file
+
+        Returns
+        -------
+        dictionary of attributes
+        """
+        debug(f'Reading BODC header from "{fn_bodc}"')
+        fid = open(fn_bodc)
+
+        # Read lines one by one (hopefully formatting is consistent)
+        # Geographical stuff
+        header_dict = {}
+        header = True
+        for line in fid:
+            if ":" in line and header == True:
+                (key, val) = line.split(":")
+                key = key.lower().strip().replace(" ", "_")
+                val = val.lower().strip().replace(" ", "_")
+                header_dict[key] = val
+                debug(f"Header key: {key} and value: {val}")
+            else:
+                # debug('No colon')
+                header = False
+        header_dict["site_name"] = header_dict["site"]  # duplicate as standard name
+        debug(f'Read done, close file "{fn_bodc}"')
+        fid.close()
+
+        header_dict["latitude"] = float(header_dict["latitude"])
+        header_dict["longitude"] = float(header_dict["longitude"])
+
+        return header_dict
+
+    @staticmethod
+    def read_bodc_data(fn_bodc, date_start=None, date_end=None, header_length: int = 11):
+        """
+        Reads observation data from a BODC file.
+
+        Parameters
+        ----------
+        fn_bodc (str) : path to bodc tide gauge file
+        date_start (datetime) : start date for returning data
+        date_end (datetime) : end date for returning data
+        header_length (int) : number of lines in header (to skip when reading)
+
+        Returns
+        -------
+        xarray.Dataset containing times, sealevel and quality control flags
+        """
+        # Initialise empty dataset and lists
+        debug(f'Reading BODC data from "{fn_bodc}"')
+        dataset = xr.Dataset()
+        time = []
+        sea_level = []
+        qc_flags = []
+        residual = []
+        # Open file and loop until EOF
+        with open(fn_bodc) as file:
+            line_count = 1
+            for line in file:
+                # Read all data. Date boundaries are set later.
+                if line_count > header_length:
+                    try:
+                        working_line = line.split()
+                        time_str = working_line[1] + " " + working_line[2]  # Empty lines cause trouble
+                        sea_level_str = working_line[3]
+                        residual_str = working_line[4]
+                        if sea_level_str[-1].isalpha():
+                            qc_flag_str = sea_level_str[-1]
+                            sea_level_str = sea_level_str.replace(qc_flag_str, "")
+                            residual_str = residual_str.replace(qc_flag_str, "")
+                        elif residual_str[-1].isalpha():  # sometimes residual has a
+                            # flag when elevation does not
+                            qc_flag_str = residual_str[-1]
+                            sea_level_str = sea_level_str.replace(qc_flag_str, "")
+                            residual_str = residual_str.replace(qc_flag_str, "")
+                        else:
+                            qc_flag_str = ""
+                        # debug(line_count-header_length, residual_str, float(residual_str))
+                        # debug(working_line, sea_level_str, qc_flag_str)
+                        time.append(time_str)
+                        qc_flags.append(qc_flag_str)
+                        sea_level.append(float(sea_level_str))
+                        residual.append(float(residual_str))
+                    except:
+                        debug(f"{file} probably empty line at end. Breaks split()")
+                line_count = line_count + 1
+            debug(f'Read done, close file "{fn_bodc}"')
+
+        # Convert time list to datetimes using pandas
+        time = np.array(pd.to_datetime(time))
+
+        # Return only values between stated dates
+        start_index = 0
+        end_index = len(time)
+        if date_start is not None:
+            date_start = np.datetime64(date_start)
+            start_index = np.argmax(time >= date_start)
+        if date_end is not None:
+            date_end = np.datetime64(date_end)
+            end_index = np.argmax(time > date_end)
+        time = time[start_index:end_index]
+        sea_level = sea_level[start_index:end_index]
+        qc_flags = qc_flags[start_index:end_index]
+
+        # Set null values to nan
+        sea_level = np.array(sea_level)
+        qc_flags = np.array(qc_flags)
+        # sea_level[qc_flags==5] = np.nan
+
+        # Assign arrays to Dataset
+        dataset["sea_level"] = xr.DataArray(sea_level, dims=["time"])
+        dataset["qc_flags"] = xr.DataArray(qc_flags, dims=["time"])
+        dataset = dataset.assign_coords(time=("time", time))
+
+        # Assign local dataset to object-scope dataset
+        return dataset
+
+    ##############################################################################
+    ###                ~            Plotting             ~                     ###
+    ##############################################################################
+
+    def plot_on_map(self):
+        """
+        Show the location of a tidegauge on a map.
+
+        Example usage:
+        --------------
+        # For a TIDEGAUGE object tg
+        tg.plot_map()
+
+        """
+
+        debug(f"Plotting tide gauge locations for {get_slug(self)}")
+
+        title = "Location: " + self.dataset.attrs["site_name"]
+        X = self.dataset.longitude
+        Y = self.dataset.latitude
+        fig, ax = plot_util.geo_scatter(X, Y, title=title)
+        ax.set_xlim((X - 10, X + 10))
+        ax.set_ylim((Y - 10, Y + 10))
+        return fig, ax
+
+    @classmethod
+    def plot_on_map_multiple(cls, tidegauge_list, color_var_str=None):
+        """
+        Show the location of a tidegauge on a map.
+
+        Example usage:
+        --------------
+        # For a TIDEGAUGE object tg
+        tg.plot_map()
+
+        """
+
+        debug(f"Plotting tide gauge locations for {get_slug(cls)}")
+
+        X = []
+        Y = []
+        C = []
+        for tg in tidegauge_list:
+            X.append(tg.dataset.longitude)
+            Y.append(tg.dataset.latitude)
+            if color_var_str is not None:
+                C.append(tg.dataset[color_var_str].values)
+
+        title = ""
+        if color_var_str is 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
+
+    def plot_timeseries(self, var_list=["sea_level"], date_start=None, date_end=None, plot_line=False):
+        """
+        Quick plot of time series stored within object's dataset
+        Parameters
+        ----------
+        date_start (datetime) : Start date for plotting
+        date_end (datetime) : End date for plotting
+        var_list  (str)  : List of variables to plot. Default: just sea_level
+        plot_line (bool) : If true, draw line between markers
+
+        Returns
+        -------
+        matplotlib figure and axes objects
+        """
+        debug(f"Plotting timeseries for {get_slug(self)}")
+        fig = plt.figure(figsize=(10, 10))
+        # Check input is a list (even for one variable)
+        if type(var_list) is str:
+            var_list = [var_list]
+
+        for var_str in var_list:
+            dim_str = self.dataset[var_str].dims[0]
+            x = np.array(self.dataset[dim_str])
+            y = np.array(self.dataset[var_str])
+
+            # Use only values between stated dates
+            start_index = 0
+            end_index = len(x)
+            if date_start is not None:
+                date_start = np.datetime64(date_start)
+                start_index = np.argmax(x >= date_start)
+            if date_end is not None:
+                date_end = np.datetime64(date_end)
+                end_index = np.argmax(x > date_end)
+            x = x[start_index:end_index]
+            y = y[start_index:end_index]
+
+            # Plot lines first if needed
+            if plot_line:
+                plt.plot(x, y, c=[0.5, 0.5, 0.5], linestyle="--", linewidth=0.5)
+
+            ax = plt.scatter(x, y, s=10)
+
+        plt.grid()
+        plt.xticks(rotation=45)
+        plt.legend(var_list)
+        # Title and axes
+        plt.xlabel("Date")
+        plt.title("Site: " + self.dataset.site_name)
+
+        return fig, ax
+
+    ##############################################################################
+    ###                ~        Model Comparison         ~                     ###
+    ##############################################################################
+
+    def obs_operator(self, model, mod_var_name: str, time_interp="nearest", model_mask=None):
+        """
+        Interpolates a model array (specified using a model object and variable
+        string) to TIDEGAUGE location and times. Takes the nearest model grid
+        cell to the tide gauge.
+
+        Parameters
+        ----------
+        model : MODEL object (e.g. NEMO)
+        model_var_name (str) : Name of variable (inside MODEL) to interpolate.
+        time_interp (str) : type of scipy time interpolation (e.g. linear)
+        model_mask : Mask to apply to model data in geographical interpolation
+                     of model. For example, use to ignore land points.
+                     If None, no mask is applied. If 'bathy', model variable
+                     (bathymetry==0) is used. Custom 2D mask arrays can be
+                     supplied.
+
+        Returns
+        -------
+        Saves interpolated array to TIDEGAUGE.dataset
+        """
+        # Determine mask
+        if model_mask == "bathy":
+            model_mask = model.dataset.bathymetry.values == 0
+
+        # Get data arrays
+        mod_var_array = model.dataset[mod_var_name]
+
+        # Depth interpolation -> for now just take 0 index
+        if "z_dim" in mod_var_array.dims:
+            mod_var_array = mod_var_array.isel(z_dim=0).squeeze()
+
+        # Cast lat/lon to numpy arrays
+        obs_lon = np.array([self.dataset.longitude])
+        obs_lat = np.array([self.dataset.latitude])
+
+        interpolated = model.interpolate_in_space(mod_var_array, obs_lon, obs_lat, mask=model_mask)
+
+        interpolated = model.interpolate_in_time(interpolated, self.dataset.time)
+
+        # Store interpolated array in dataset
+        new_var_name = "interp_" + mod_var_name
+        self.dataset[new_var_name] = interpolated.drop(["longitude", "latitude"])
+        return
+
+    def crps(
+        self,
+        model_object,
+        model_var_name,
+        obs_var_name: str = "sea_level",
+        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.
+        This method specifically performs a single-observation neighbourhood-
+        forecast method.
+
+        Parameters
+        ----------
+        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.
+        time_interp (str)    : Type of time interpolation to use (s)
+        create_new_obj (bool): If True, save output to new TIDEGAUGE obj.
+                               Otherwise, save to this obj.
+
+        Returns
+        -------
+        xarray.Dataset containing times, sealevel and quality control flags
+
+        Example Useage
+        -------
+        # Compare modelled 'sossheig' with 'sea_level' using CRPS
+        crps = altimetry.crps(nemo, 'sossheig', 'sea_level')
+        """
+
+        mod_var = model_object.dataset[model_var_name]
+        obs_var = self.dataset[obs_var_name]
+
+        crps_list, n_model_pts, contains_land = crps_util.crps_sonf_fixed(
+            mod_var,
+            self.dataset.longitude.values,
+            self.dataset.latitude.values,
+            obs_var.values,
+            obs_var.time.values,
+            nh_radius,
+            time_interp,
+        )
+        if create_new_obj:
+            new_object = Tidegauge()
+            new_dataset = self.dataset[["longitude", "latitude", "time"]]
+            new_dataset["crps"] = (("time"), crps_list)
+            new_dataset["crps_n_model_pts"] = (("time"), n_model_pts)
+            new_object.dataset = new_dataset
+            return new_object
+        else:
+            self.dataset["crps"] = (("time"), crps_list)
+            self.dataset["crps_n_model_pts"] = (("time"), n_model_pts)
+
+    def difference(self, var_str0: str, var_str1: str, date0=None, date1=None):
+        """Difference two variables defined by var_str0 and var_str1 between
+        two dates date0 and date1. Returns xr.DataArray"""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
+        diff = var0 - var1
+        return xr.DataArray(diff, dims="time", name="error", coords={"time": self.dataset.time})
+
+    def absolute_error(self, var_str0, var_str1, date0=None, date1=None):
+        """Absolute difference two variables defined by var_str0 and var_str1
+        between two dates date0 and date1. Return xr.DataArray"""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
+        adiff = np.abs(var0 - var1)
+        return xr.DataArray(adiff, dims="time", name="absolute_error", coords={"time": self.dataset.time})
+
+    def mean_absolute_error(self, var_str0, var_str1, date0=None, date1=None):
+        """Mean absolute difference two variables defined by var_str0 and
+        var_str1 between two dates date0 and date1. Return xr.DataArray"""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
+        mae = metrics.mean_absolute_error(var0, var1)
+        return mae
+
+    def root_mean_square_error(self, var_str0, var_str1, date0=None, date1=None):
+        """Root mean square difference two variables defined by var_str0 and
+        var_str1 between two dates date0 and date1. Return xr.DataArray"""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = general_utils.data_array_time_slice(var0, date0, date1).values
+        var1 = general_utils.data_array_time_slice(var1, date0, date1).values
+        rmse = metrics.mean_squared_error(var0, var1)
+        return np.sqrt(rmse)
+
+    def time_mean(self, var_str, date0=None, date1=None):
+        """Time mean of variable var_str between dates date0, date1"""
+        var = self.dataset[var_str]
+        var = general_utils.data_array_time_slice(var, date0, date1)
+        return np.nanmean(var)
+
+    def time_std(self, var_str, date0=None, date1=None):
+        """Time st. dev of variable var_str between dates date0 and date1"""
+        var = self.dataset[var_str]
+        var = general_utils.data_array_time_slice(var, date0, date1)
+        return np.nanstd(var)
+
+    def time_correlation(self, var_str0, var_str1, date0=None, date1=None, method="pearson"):
+        """Time correlation between two variables defined by var_str0,
+        var_str1 between dates date0 and date1. Uses Pandas corr()."""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = var0.rename("var1")
+        var1 = var1.rename("var2")
+        var0 = general_utils.data_array_time_slice(var0, date0, date1)
+        var1 = general_utils.data_array_time_slice(var1, date0, date1)
+        pdvar = xr.merge((var0, var1))
+        pdvar = pdvar.to_dataframe()
+        corr = pdvar.corr(method=method)
+        return corr.iloc[0, 1]
+
+    def time_covariance(self, var_str0, var_str1, date0=None, date1=None):
+        """Time covariance between two variables defined by var_str0,
+        var_str1 between dates date0 and date1. Uses Pandas corr()."""
+        var0 = self.dataset[var_str0]
+        var1 = self.dataset[var_str1]
+        var0 = var0.rename("var1")
+        var1 = var1.rename("var2")
+        var0 = general_utils.data_array_time_slice(var0, date0, date1)
+        var1 = general_utils.data_array_time_slice(var1, date0, date1)
+        pdvar = xr.merge((var0, var1))
+        pdvar = pdvar.to_dataframe()
+        cov = pdvar.cov()
+        return cov.iloc[0, 1]
+
+    def basic_stats(self, var_str0, var_str1, date0=None, date1=None, create_new_object=True):
+        """Calculates a selection of statistics for two variables defined by
+        var_str0 and var_str1, between dates date0 and date1. This will return
+        their difference, absolute difference, mean absolute error, root mean
+        square error, correlation and covariance. If create_new_object is True
+        then this method returns a new TIDEGAUGE object containing statistics,
+        otherwise variables are saved to the dateset inside this object."""
+
+        diff = self.difference(var_str0, var_str1, date0, date1)
+        ae = self.absolute_error(var_str0, var_str1, date0, date1)
+        mae = self.mean_absolute_error(var_str0, var_str1, date0, date1)
+        rmse = self.root_mean_square_error(var_str0, var_str1, date0, date1)
+        corr = self.time_correlation(var_str0, var_str1, date0, date1)
+        cov = self.time_covariance(var_str0, var_str1, date0, date1)
+
+        if create_new_object:
+            new_object = Tidegauge()
+            new_dataset = self.dataset[["longitude", "latitude", "time"]]
+            new_dataset["absolute_error"] = ae
+            new_dataset["error"] = diff
+            new_dataset["mae"] = mae
+            new_dataset["rmse"] = rmse
+            new_dataset["corr"] = corr
+            new_dataset["cov"] = cov
+            new_object.dataset = new_dataset
+            return new_object
+        else:
+            self.dataset["absolute_error"] = ae
+            self.dataset["error"] = diff
+            self.dataset["mae"] = mae
+            self.dataset["rmse"] = rmse
+            self.dataset["corr"] = corr
+            self.dataset["cov"] = cov
+
+    ##############################################################################
+    ###                ~            Analysis             ~                     ###
+    ##############################################################################
+
+    def resample_mean(self, var_str: str, time_freq: str, **kwargs):
+        """Resample a TIDEGAUGE variable in time by calculating the mean
+            of all data points at a given frequency.
+
+        Parameters
+        ----------
+        var_str (str)    : Variable name to resample
+        time_freq (str)  : Time frequency. e.g. '1H' for hourly, '1D' for daily
+                           Can also be a timedelta object. See Pandas resample
+                           method for more info.
+        **kwargs (other) : Other arguments to pass to xarray.Dataset.resample
+        (http://xarray.pydata.org/en/stable/generated/xarray.Dataset.resample.html)
+
+        Returns
+        -------
+        New variable (var_str_freq) and dimension (time_freq) in tg.dataset
+        """
+        # Define new variable and dimension names
+        var = self.dataset[var_str]
+        new_var_str = var_str + "_" + time_freq
+        new_dim_str = "time_" + time_freq
+
+        # Resample using xarray.resample
+        resampled = var.resample(time=time_freq, **kwargs).mean()
+
+        # Rename dimensions and variables. Put into original dataset.
+        resampled = resampled.rename({"time": new_dim_str})
+        resampled = resampled.rename(new_var_str)
+        self.dataset[new_var_str] = resampled
+
+    def apply_doodson_x0_filter(self, var_str):
+        """Applies doodson X0 filter to a specified TIDEGAUGE variable
+        Input ius expected to be hourly. Use resample_mean to average data
+        to hourly frequency."""
+        filtered = stats_util.doodson_x0_filter(self.dataset[var_str], ax=0)
+        self.dataset[var_str + "_dx0"] = (("time_1H"), filtered)
+
+    def find_high_and_low_water(self, var_str, method="comp", **kwargs):
+        """
+        Finds high and low water for a given variable.
+        Returns in a new TIDEGAUGE object with similar data format to
+        a TIDETABLE.
+
+        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.
+        '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
+        y = self.dataset[var_str]
+
+        time_max, values_max = stats_util.find_maxima(x, y, method=method, **kwargs)
+        time_min, values_min = stats_util.find_maxima(x, -y, method=method, **kwargs)
+
+        new_dataset = xr.Dataset()
+        new_dataset.attrs = self.dataset.attrs
+        new_dataset[var_str + "_highs"] = ("time_highs", values_max)
+        new_dataset[var_str + "_lows"] = ("time_lows", -values_min)
+        new_dataset["time_highs"] = ("time_highs", time_max)
+        new_dataset["time_lows"] = ("time_lows", time_min)
+
+        new_object = Tidegauge()
+        new_object.dataset = new_dataset
+
+        return new_object
diff --git a/coast/timeseries.py b/coast/timeseries.py
new file mode 100644
index 00000000..2e069fd4
--- /dev/null
+++ b/coast/timeseries.py
@@ -0,0 +1,10 @@
+"""Timeseries Class"""
+from .index import Indexed
+
+
+class Timeseries(Indexed):
+    """Parent class for Tidegauge and other timeseries type datasets
+    Common methods ...
+    """
+
+    pass
diff --git a/coast/track.py b/coast/track.py
new file mode 100644
index 00000000..b13c5ccc
--- /dev/null
+++ b/coast/track.py
@@ -0,0 +1,10 @@
+"""Track class"""
+from .index import Indexed
+
+
+class Track(Indexed):
+    """Parent class for subclasses Altimetry
+    Common methods ....
+    """
+
+    pass
diff --git a/config/example_altimetry.json b/config/example_altimetry.json
new file mode 100644
index 00000000..a3d44b25
--- /dev/null
+++ b/config/example_altimetry.json
@@ -0,0 +1,18 @@
+{
+	"type": "indexed",
+	"dimensionality": 1,
+	"chunks": {"time":1000}
+	,
+	"dataset": {
+		"dimension_map": {
+			"time": "t_dim"
+		},
+		"variable_map": {
+			"ocean_tide": "ocean_tide_standard_name"
+		},
+		"coord_vars": []
+	},
+	"processing_flags": [
+		"example_flag"
+	]
+}
\ No newline at end of file
diff --git a/config/example_argos.json b/config/example_argos.json
new file mode 100644
index 00000000..6cef3e4b
--- /dev/null
+++ b/config/example_argos.json
@@ -0,0 +1,23 @@
+{
+	"type": "indexed",
+	"dimensionality": 1,
+	"chunks": {},
+	"dataset": {
+		"dimension_map": {
+			"time": "t_time"
+		},
+		"variable_map_comment":"mapping from external variable name to standard variable names see https://pcmdi.llnl.gov/mips/cmip3/variableList.html",
+		"variable_map": {
+			"DATIM": "time",
+			"LAT": "latitude",
+			"LON": "longitude",
+			"SST": "sea_surface_temperature",
+			"PPRES": "pressure"
+		},
+		"coord_vars": ["time","latitude","longitude"],
+		"keep_all_vars": "False"
+	},
+	"processing_flags": [
+		"example_flag"
+	]
+}
\ No newline at end of file
diff --git a/config/example_en4_profiles.json b/config/example_en4_profiles.json
new file mode 100644
index 00000000..38c330b5
--- /dev/null
+++ b/config/example_en4_profiles.json
@@ -0,0 +1,29 @@
+{
+	"type": "indexed",
+	"dimensionality": 2,
+	"chunks": {},
+	"dataset": {
+		"dimension_map": {
+			"N_PROF": "profile",
+			"N_LEVELS": "z_dim"
+		},
+		"variable_map": {
+			"LATITUDE": "latitude",
+			"LONGITUDE": "longitude",
+			"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"
+		},
+		"coord_vars": ["time","latitude","longitude"],
+		"keep_all_vars": "False"
+	},
+	"processing_flags": [
+		"example_flag"
+	]
+}
diff --git a/config/example_glider_ego.json b/config/example_glider_ego.json
new file mode 100644
index 00000000..eb06ee28
--- /dev/null
+++ b/config/example_glider_ego.json
@@ -0,0 +1,28 @@
+{
+    "type": "indexed",
+    "dimensionality": 1,
+    "chunks": {},
+    "dataset": {
+        "dimension_map": {
+            "ex0": "N_HISTORY",
+            "ex1": "N_TRANS_SYSTEM",
+            "ex2": "N_POSITIONING_SYSTEM",
+            "N_PARAM": "N_PARAM_STANDARD",
+            "ex4": "N_DERIVATION",
+            "TIME": "TIME_STANDARD",
+            "ex6": "TIME_GPS"
+        },
+        "variable_map": {
+            "PLATFORM_TYPE": "PLATFORM_TYPE",
+            "PLATFORM_MAKER": "PLATFORM_MAKER",
+            "LATITUDE_GPS": "LATITUDE_GPS",
+            "ex68": "LATITUDE_GPS_QC",
+            "LONGITUDE_GPS": "LONGITUDE_GPS",
+            "TEMP": "TEMP"
+        },
+        "coord_vars": []
+    },
+    "processing_flags": [
+        "example_flag"
+    ]
+}
diff --git a/config/example_glider_ego2.json b/config/example_glider_ego2.json
new file mode 100644
index 00000000..0eb4c23b
--- /dev/null
+++ b/config/example_glider_ego2.json
@@ -0,0 +1,125 @@
+{
+    "type": "indexed",
+    "dimensionality": 1,
+    "chunks": [],
+    "dataset": {
+        "dimension_map": {
+            "ex0": "N_HISTORY",
+            "ex1": "N_TRANS_SYSTEM",
+            "ex2": "N_POSITIONING_SYSTEM",
+            "N_PARAM": "N_PARAM_STANDARD",
+            "ex4": "N_DERIVATION",
+            "TIME": "TIME_STANDARD",
+            "ex6": "TIME_GPS"
+        },
+        "variable_map": {
+            "PLATFORM_TYPE": "PLATFORM_TYPE",
+            "PLATFORM_MAKER": "PLATFORM_MAKER",
+            "ex2": "WMO_INST_TYPE",
+            "ex3": "PLATFORM_FAMILY",
+            "ex4": "GLIDER_SERIAL_NO",
+            "ex5": "GLIDER_OWNER",
+            "ex6": "HISTORY_INSTITUTION",
+            "ex7": "HISTORY_STEP",
+            "ex8": "HISTORY_SOFTWARE",
+            "ex9": "HISTORY_SOFTWARE_RELEASE",
+            "ex10": "HISTORY_REFERENCE",
+            "ex11": "HISTORY_DATE",
+            "ex12": "HISTORY_ACTION",
+            "ex13": "HISTORY_PARAMETER",
+            "ex14": "HISTORY_PREVIOUS_VALUE",
+            "ex15": "HISTORY_START_TIME_INDEX",
+            "ex16": "HISTORY_STOP_TIME_INDEX",
+            "ex17": "HISTORY_QCTEST",
+            "ex18": "TRANS_SYSTEM",
+            "ex19": "TRANS_SYSTEM_ID",
+            "ex20": "TRANS_FREQUENCY",
+            "ex21": "POSITIONING_SYSTEM",
+            "ex22": "FIRMWARE_VERSION_NAVIGATION",
+            "ex23": "FIRMWARE_VERSION_SCIENCE",
+            "ex24": "MANUAL_VERSION",
+            "ex25": "STANDARD_FORMAT_ID",
+            "ex26": "DAC_FORMAT_ID",
+            "ex27": "PROJECT_NAME",
+            "ex28": "DATA_CENTRE",
+            "ex29": "PI_NAME",
+            "ex30": "ANOMALY",
+            "ex31": "BATTERY_TYPE",
+            "ex32": "BATTERY_PACKS",
+            "ex33": "SPECIAL_FEATURES",
+            "ex34": "OPERATING_INSTITUTION",
+            "ex35": "CUSTOMIZATION",
+            "ex36": "DEPLOYMENT_START_DATE",
+            "ex37": "DEPLOYMENT_START_LATITUDE",
+            "ex38": "DEPLOYMENT_START_LONGITUDE",
+            "ex39": "DEPLOYMENT_START_QC",
+            "ex40": "DEPLOYMENT_PLATFORM",
+            "ex41": "DEPLOYMENT_CRUISE_ID",
+            "ex42": "DEPLOYMENT_REFERENCE_STATION_ID",
+            "ex43": "DEPLOYMENT_END_DATE",
+            "ex44": "DEPLOYMENT_END_LATITUDE",
+            "ex45": "DEPLOYMENT_END_LONGITUDE",
+            "ex46": "DEPLOYMENT_END_QC",
+            "ex47": "DEPLOYMENT_END_STATUS",
+            "ex48": "DEPLOYMENT_OPERATOR",
+            "ex49": "SENSOR",
+            "ex50": "SENSOR_MAKER",
+            "ex51": "SENSOR_MODEL",
+            "ex52": "SENSOR_SERIAL_NO",
+            "ex53": "SENSOR_UNITS",
+            "ex54": "SENSOR_ACCURACY",
+            "ex55": "SENSOR_RESOLUTION",
+            "ex56": "DERIVATION_PARAMETER",
+            "ex57": "DERIVATION_EQUATION",
+            "ex58": "DERIVATION_COEFFICIENT",
+            "ex59": "DERIVATION_COMMENT",
+            "ex60": "DERIVATION_DATE",
+            "ex61": "TIME_QC",
+            "ex62": "JULD",
+            "ex63": "JULD_QC",
+            "ex64": "POSITION_QC",
+            "ex65": "TIME_GPS_QC",
+            "ex66": "POSITION_GPS_QC",
+            "LATITUDE_GPS": "LATITUDE_GPS",
+            "ex68": "LATITUDE_GPS_QC",
+            "LONGITUDE_GPS": "LONGITUDE_GPS",
+            "ex70": "LONGITUDE_GPS_QC",
+            "ex71": "ALTITUDE",
+            "ex72": "ALTITUDE_QC",
+            "ex73": "ALTITUDE_UNCERTAINTY",
+            "ex74": "LATITUDE_QC",
+            "ex75": "LATITUDE_UNCERTAINTY",
+            "ex76": "LONGITUDE_QC",
+            "ex77": "LONGITUDE_UNCERTAINTY",
+            "ex78": "GLIDER_PITCH",
+            "ex79": "GLIDER_PITCH_QC",
+            "ex80": "GLIDER_PITCH_UNCERTAINTY",
+            "ex81": "GLIDER_ROLL",
+            "ex82": "GLIDER_ROLL_QC",
+            "ex83": "GLIDER_ROLL_UNCERTAINTY",
+            "ex84": "BBP700",
+            "ex85": "BBP700_QC",
+            "ex86": "BBP700_UNCERTAINTY",
+            "ex87": "CDOM",
+            "ex88": "CDOM_QC",
+            "ex89": "CDOM_UNCERTAINTY",
+            "ex90": "CHLA",
+            "ex91": "CHLA_QC",
+            "ex92": "CHLA_UNCERTAINTY",
+            "ex93": "CNDC",
+            "ex94": "CNDC_QC",
+            "ex95": "CNDC_UNCERTAINTY",
+            "ex96": "PRES_QC",
+            "ex97": "PRES_UNCERTAINTY",
+            "TEMP": "TEMP",
+            "ex99": "TEMP_QC",
+            "ex100": "TEMP_UNCERTAINTY",
+            "ex101": "PHASE",
+            "ex102": "PHASE_NUMBER",
+            "ex103": "POSITIONING_METHOD"
+        }
+    },
+    "processing_flags": [
+        "example_flag"
+    ]
+}
diff --git a/config/example_gridded.json b/config/example_gridded.json
index 4aa94763..6fe579e3 100644
--- a/config/example_gridded.json
+++ b/config/example_gridded.json
@@ -40,7 +40,13 @@
 			"ex12": "vozocrtx",
 			"ex13": "sossheig",
 			"ex14": "zos"
-		}
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
 	},
 	"domain": {
 		"dimension_map": {
@@ -95,12 +101,6 @@
 			"ln_sco",
 			"ln_isfcav"
 		],
-		"coord_vars": [
-			"longitude",
-			"latitude",
-			"time",
-			"depth_0"
-		],
 		"delete_vars": [
 			"nav_lat",
 			"nav_lon",
diff --git a/config/example_indexed.json b/config/example_indexed.json
index 113124d2..a37c46ff 100644
--- a/config/example_indexed.json
+++ b/config/example_indexed.json
@@ -28,7 +28,8 @@
 			"ex12": "vozocrtx",
 			"ex13": "sossheig",
 			"ex14": "zos"
-		}
+		},
+		"coord_vars": []
 	},
 	"processing_flags": [
 		"example_flag"
diff --git a/config/example_nemo_grid_f.json b/config/example_nemo_grid_f.json
new file mode 100644
index 00000000..65899560
--- /dev/null
+++ b/config/example_nemo_grid_f.json
@@ -0,0 +1,76 @@
+{
+	"type": "gridded",
+	"dimensionality": 3,
+	"chunks": {},
+	"grid_ref": {
+		"f-grid": [
+			"glamf",
+			"gphif",
+			"e1f",
+			"e2f",
+			"e3f_0",
+			"depthf_0"
+		]
+	},
+	"dataset": {
+		"dimension_map": {
+			"time_counter": "t_dim",
+			"deptht": "z_dim",
+			"depthu": "z_dim",
+			"depthv": "z_dim",
+			"y": "y_dim",
+			"x": "x_dim",
+			"x_grid_T": "x_dim",
+			"y_grid_T": "y_dim"
+		},
+		"variable_map": {
+			"time_counter": "time"
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
+	},
+	"domain": {
+		"dimension_map": {
+			"t": "t_dim0",
+			"x": "x_dim",
+			"y": "y_dim",
+			"z": "z_dim"
+		},
+		"variable_map": {
+			"time_counter": "time0",
+			"glamf": "longitude",
+			"gphif": "latitude",
+			"e1f": "e1",
+			"e2f": "e2",
+			"ff_f": "ff",
+			"e3f_0": "e3_0",
+			"depthf_0": "depth_0",
+			"bottom_level": "bottom_level"
+		}
+	},
+	"static_variables": {
+		"not_grid_vars": [
+			"jpiglo",
+			"jpjglo",
+			"jpkglo",
+			"jperio",
+			"ln_zco",
+			"ln_zps",
+			"ln_sco",
+			"ln_isfcav"
+		],
+		"delete_vars": [
+			"nav_lat",
+			"nav_lon",
+			"deptht"
+		]
+	},
+	"processing_flags": [
+		"example_flag1",
+		"example_flag2"
+	]
+}
diff --git a/config/example_nemo_grid_t.json b/config/example_nemo_grid_t.json
new file mode 100644
index 00000000..49fb345a
--- /dev/null
+++ b/config/example_nemo_grid_t.json
@@ -0,0 +1,84 @@
+{
+	"type": "gridded",
+	"dimensionality": 3,
+	"chunks": {},
+	"grid_ref": {
+		"t-grid": [
+			"glamt",
+			"gphit",
+			"e1t",
+			"e2t",
+			"e3t_0",
+			"deptht_0",
+			"tmask"
+		]
+	},
+	"dataset": {
+		"dimension_map": {
+			"time_counter": "t_dim",
+			"deptht": "z_dim",
+			"y": "y_dim",
+			"x": "x_dim",
+			"x_grid_T": "x_dim",
+			"y_grid_T": "y_dim"
+		},
+		"variable_map": {
+			"time_counter": "time",
+			"votemper": "temperature",
+			"thetao": "temperature",
+			"temp": "temperature",
+			"toce": "temperature",
+			"so": "salinity",
+			"vosaline": "salinity",
+			"soce": "salinity",
+			"sossheig": "ssh",
+			"zos": "ssh"
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
+	},
+	"domain": {
+		"dimension_map": {
+			"t": "t_dim0",
+			"x": "x_dim",
+			"y": "y_dim",
+			"z": "z_dim"
+		},
+		"variable_map": {
+			"time_counter": "time0",
+			"glamt": "longitude",
+			"gphit": "latitude",
+			"e1t": "e1",
+			"e2t": "e2",
+			"e3t_0": "e3_0",
+			"tmask":"mask",
+			"deptht_0": "depth_0",
+			"bottom_level": "bottom_level"
+		}
+	},
+	"static_variables": {
+		"not_grid_vars": [
+			"jpiglo",
+			"jpjglo",
+			"jpkglo",
+			"jperio",
+			"ln_zco",
+			"ln_zps",
+			"ln_sco",
+			"ln_isfcav"
+		],
+		"delete_vars": [
+			"nav_lat",
+			"nav_lon",
+			"deptht"
+		]
+	},
+	"processing_flags": [
+		"example_flag1",
+		"example_flag2"
+	]
+}
diff --git a/config/example_nemo_grid_u.json b/config/example_nemo_grid_u.json
new file mode 100644
index 00000000..0d2a1716
--- /dev/null
+++ b/config/example_nemo_grid_u.json
@@ -0,0 +1,75 @@
+{
+	"type": "gridded",
+	"dimensionality": 3,
+	"chunks": {},
+	"grid_ref": {
+		"u-grid": [
+			"glamu",
+			"gphiu",
+			"e1u",
+			"e2u",
+			"e3u_0",
+			"depthu_0",
+			"umask"
+		]
+	},
+	"dataset": {
+		"dimension_map": {
+			"time_counter": "t_dim",
+			"depthu": "z_dim",
+			"y": "y_dim",
+			"x": "x_dim"
+		},
+		"variable_map": {
+			"time_counter": "time",
+			"uoce": "u_velocity",
+			"vozocrtx": "u_velocity"
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
+	},
+	"domain": {
+		"dimension_map": {
+			"t": "t_dim0",
+			"x": "x_dim",
+			"y": "y_dim",
+			"z": "z_dim"
+		},
+		"variable_map": {
+			"time_counter": "time0",
+			"glamu": "longitude",
+			"gphiu": "latitude",
+			"e1u": "e1",
+			"e2u": "e2",
+			"e3u_0": "e3_0",
+			"depthu_0": "depth_0",
+			"bottom_level": "bottom_level",
+			"umask": "mask"
+		}
+	},
+	"static_variables": {
+		"not_grid_vars": [
+			"jpiglo",
+			"jpjglo",
+			"jpkglo",
+			"jperio",
+			"ln_zco",
+			"ln_zps",
+			"ln_sco",
+			"ln_isfcav"
+		],
+		"delete_vars": [
+			"nav_lat",
+			"nav_lon",
+			"deptht"
+		]
+	},
+	"processing_flags": [
+		"example_flag1",
+		"example_flag2"
+	]
+}
diff --git a/config/example_nemo_grid_v.json b/config/example_nemo_grid_v.json
new file mode 100644
index 00000000..24a4279d
--- /dev/null
+++ b/config/example_nemo_grid_v.json
@@ -0,0 +1,75 @@
+{
+	"type": "gridded",
+	"dimensionality": 3,
+	"chunks": {},
+	"grid_ref": {
+		"v-grid": [
+			"glamv",
+			"gphiv",
+			"e1v",
+			"e2v",
+			"e3v_0",
+			"depthv_0",
+			"vmask"
+		]
+	},
+	"dataset": {
+		"dimension_map": {
+			"time_counter": "t_dim",
+			"depthv": "z_dim",
+			"y": "y_dim",
+			"x": "x_dim"
+		},
+		"variable_map": {
+			"time_counter": "time",
+			"voce": "v_velocity",
+			"vomecrty": "v_velocity"
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
+	},
+	"domain": {
+		"dimension_map": {
+			"t": "t_dim0",
+			"x": "x_dim",
+			"y": "y_dim",
+			"z": "z_dim"
+		},
+		"variable_map": {
+			"time_counter": "time0",
+			"glamv": "longitude",
+			"gphiv": "latitude",
+			"e1v": "e1",
+			"e2v": "e2",
+			"e3v_0": "e3_0",
+			"depthv_0": "depth_0",
+			"bottom_level": "bottom_level",
+			"vmask": "mask"
+		}
+	},
+	"static_variables": {
+		"not_grid_vars": [
+			"jpiglo",
+			"jpjglo",
+			"jpkglo",
+			"jperio",
+			"ln_zco",
+			"ln_zps",
+			"ln_sco",
+			"ln_isfcav"
+		],
+		"delete_vars": [
+			"nav_lat",
+			"nav_lon",
+			"deptht"
+		]
+	},
+	"processing_flags": [
+		"example_flag1",
+		"example_flag2"
+	]
+}
diff --git a/config/example_nemo_grid_w.json b/config/example_nemo_grid_w.json
new file mode 100644
index 00000000..9a668201
--- /dev/null
+++ b/config/example_nemo_grid_w.json
@@ -0,0 +1,80 @@
+{
+	"type": "gridded",
+	"dimensionality": 3,
+	"chunks": {},
+	"grid_ref": {
+		"w-grid": [
+			"glamt",
+			"gphit",
+			"e1t",
+			"e2t",
+			"e3w_0",
+			"depthw_0"
+		]
+	},
+	"dataset": {
+		"dimension_map": {
+			"time_counter": "t_dim",
+			"deptht": "z_dim",
+			"depthu": "z_dim",
+			"depthv": "z_dim",
+			"y": "y_dim",
+			"x": "x_dim",
+			"x_grid_T": "x_dim",
+			"y_grid_T": "y_dim"
+		},
+		"variable_map": {
+			"time_counter": "time"
+		},
+		"coord_vars": [
+			"longitude",
+			"latitude",
+			"time",
+			"depth_0"
+		]
+	},
+	"domain": {
+		"dimension_map": {
+			"t": "t_dim0",
+			"x": "x_dim",
+			"y": "y_dim",
+			"z": "z_dim"
+		},
+		"variable_map": {
+			"time_counter": "time0",
+			"glamt": "longitude",
+			"glamf": "longitude",
+			"gphit": "latitude",
+			"gphif": "latitude",
+			"e1t": "e1",
+			"e1f": "e1",
+			"e2t": "e2",
+			"e2f": "e2",
+			"e3w_0": "e3_0",
+			"tmask": "mask",
+			"depthw_0": "depth_0",
+			"bottom_level": "bottom_level"
+		}
+	},
+	"static_variables": {
+		"not_grid_vars": [
+			"jpiglo",
+			"jpjglo",
+			"jpkglo",
+			"jperio",
+			"ln_zco",
+			"ln_zps",
+			"ln_sco",
+			"ln_isfcav"
+		],
+		"delete_vars": [
+			"nav_lat",
+			"nav_lon",
+			"deptht"
+		]
+	},
+	"processing_flags": [
+		"example_flag1",
+		"example_flag2"
+	]
+}
diff --git a/config/example_oceanparcels.json b/config/example_oceanparcels.json
new file mode 100644
index 00000000..6678249f
--- /dev/null
+++ b/config/example_oceanparcels.json
@@ -0,0 +1,26 @@
+{
+	"type": "indexed",
+	"dimensionality": 2,
+	"chunks": {},
+	"dataset": {
+		"dimension_map": {
+			"ex0": "obs",
+			"ex1": "traj"
+		},
+		"variable_map_comment":"mapping from external variable name to standard variable names see https://pcmdi.llnl.gov/mips/cmip3/variableList.html",
+		"variable_map": {
+			"ex0": "trajectory",
+			"time": "time",
+			"lat": "latitude",
+			"lon": "longitude",
+			"z": "depth",
+			"t": "t",
+			"s": "s"
+		},
+		"coord_vars": ["time","latitude","longitude"],
+		"keep_all_vars": "False"
+	},
+	"processing_flags": [
+		"example_flag"
+	]
+}
\ No newline at end of file
diff --git a/config/example_tidegauge.json b/config/example_tidegauge.json
new file mode 100644
index 00000000..9c058c50
--- /dev/null
+++ b/config/example_tidegauge.json
@@ -0,0 +1,20 @@
+{
+	"type": "indexed",
+	"dimensionality": 1,
+	"chunks": {},
+	"dataset": {
+		"dimension_map": {
+			"time": "time"
+		},
+		"variable_map": {
+			"sea_level": "sea_level",
+			"qc_flags": "qc_flags",
+			"latitude": "latitude",
+			"longitude": "longitude"
+		},
+		"coord_vars": ["time"]
+	},
+	"processing_flags": [
+		"example_flag"
+	]
+}
diff --git a/example_scripts/altimetry_tutorial.py b/example_scripts/altimetry_tutorial.py
index 40a230de..4ae1d77e 100755
--- a/example_scripts/altimetry_tutorial.py
+++ b/example_scripts/altimetry_tutorial.py
@@ -8,18 +8,22 @@
 import coast
 
 # And by defining some file paths
-fn_nemo_dat = "./example_files/COAsT_example_NEMO_data.nc"
-fn_nemo_dom = "./example_files/COAsT_example_NEMO_domain.nc"
-fn_altimetry = "./example_files/COAsT_example_altimetry_data.nc"
+dn_files = "./example_files/"
+fn_nemo_dat = dn_files + "coast_example_nemo_data.nc"
+fn_nemo_dom = dn_files + "coast_example_nemo_domain.nc"
+fn_nemo_config = "./config/example_nemo_grid_t.json"
+fn_altimetry = dn_files + "coast_example_altimetry_data.nc"
+fn_altimetry_config = "./config/example_altimetry.json"
 
 # We need to load in a NEMO object for doing NEMO things.
-nemo = coast.Nemo(fn_nemo_dat, fn_nemo_dom, grid_ref="t-grid")
+nemo = coast.Gridded(fn_nemo_dat, fn_nemo_dom, config=fn_nemo_config)
 
 # And now we can load in our Altimetry data. By default, Altimetry is set up
 # to read in CMEMS netCDF files. However, if no path is supplied, then the
 # object's dataset will be initialised as None. Custom data can then be loaded
 # if desired, as long as it follows the data formatting for Altimetry.
-altimetry = coast.Altimetry(fn_altimetry)
+# altimetry = coast.Altimetry(fn_altimetry)
+altimetry = coast.Altimetry(fn_altimetry, config=fn_altimetry_config)
 
 # Before going any further, lets just cut out the bit of the altimetry that
 # is over the model domain. This can be done using subset_indices_lonlat_box
@@ -42,7 +46,7 @@
 # using some basic metrics. The basic_stats() routine can be used for this,
 # which calculates some simple metrics including differences, RMSE and
 # correlations. NOTE: This may not be a wise choice of variables.
-stats = altimetry.basic_stats("interp_ssh", "sla_filtered")
+stats = altimetry.basic_stats("ocean_tide_standard_name", "interp_ssh")
 
 # Take a look inside stats.dataset to see all of the new variables. When using
 # basic stats, the returned object is also an Altimetry object, so all of the
@@ -53,7 +57,7 @@
 # Probability Score (CRPS). For this, we need to hand over the model object,
 # a model variable and an observed variable. We also give it a neighbourhood
 # radius in km (nh_radius).
-crps = altimetry.crps(nemo, model_var_name="ssh", obs_var_name="sla_filtered", nh_radius=20)
+crps = altimetry.crps(nemo, model_var_name="ssh", obs_var_name="ocean_tide_standard_name", nh_radius=20)
 
 # Again, take a look inside crps.dataset to see some new variables. Similarly
 # to basic_stats, create_new_object can be set to false to save output to
@@ -62,7 +66,7 @@
 # Altimetry has a ready built quick_plot() routine for taking a look at any
 # of the observed or derived quantities above. So to take a look at the
 # 'sla_filtered' variable:
-fig, ax = altimetry.quick_plot("sla_filtered")
+fig, ax = altimetry.quick_plot("ocean_tide_standard_name")
 
 # As stats and crps are also Altimetry objects, quick_plot() can also be used:
 fig, ax = crps.quick_plot("crps")
diff --git a/example_scripts/amm15_example_plot.py b/example_scripts/amm15_example_plot.py
index dcddd512..fb8352c0 100755
--- a/example_scripts/amm15_example_plot.py
+++ b/example_scripts/amm15_example_plot.py
@@ -18,7 +18,7 @@
 dir_nam = "/projectsa/NEMO/gmaya/2013p2/"
 fil_nam = "20130415_25hourm_grid_T.nc"
 dom_nam = "/projectsa/NEMO/gmaya/AMM15_GRID/amm15.mesh_mask.cs3x.nc"
-config = "/work/jelt/GitHub/COAsT/example_files/example_t_nemo_config.json"
+config = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_t.json"
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dom_nam, config=config)  # , chunks=chunks)
 chunks = {
diff --git a/example_scripts/anchor_plots_of_nsea_wvel.py b/example_scripts/anchor_plots_of_nsea_wvel.py
index 8ab2f677..eac1aad2 100755
--- a/example_scripts/anchor_plots_of_nsea_wvel.py
+++ b/example_scripts/anchor_plots_of_nsea_wvel.py
@@ -22,7 +22,7 @@
 dir_nam = "/projectsa/NEMO/jelt/AMM60_ARCHER_DUMP/AMM60smago/EXP_NSea/OUTPUT/"
 fil_nam = "AMM60_1h_20120204_20120208_NorthSea.nc"
 
-config = "/work/jelt/GitHub/COAsT/example_files/example_w_nemo_config.json"
+config = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_w.json"
 
 chunks = {
     "x_dim": 10,
diff --git a/example_scripts/blz_example_plot.py b/example_scripts/blz_example_plot.py
index 978a37fe..6aa0b642 100755
--- a/example_scripts/blz_example_plot.py
+++ b/example_scripts/blz_example_plot.py
@@ -17,11 +17,10 @@
 dir_nam = "/projectsa/accord/GCOMS1k/OUTPUTS/BLZE12_02/2015/"
 fil_nam = "BLZE12_1h_20151101_20151130_grid_T.nc"
 dom_nam = "/projectsa/accord/GCOMS1k/INPUTS/BLZE12_C1/domain_cfg.nc"
-config_t = "/work/jelt/GitHub/COAsT/example_files/example_t_nemo_config.json"
-config_u = "/work/jelt/GitHub/COAsT/example_files/example_u_nemo_config.json"
-config_v = "/work/jelt/GitHub/COAsT/example_files/example_v_nemo_config.json"
-config_w = "/work/jelt/GitHub/COAsT/example_files/example_w_nemo_config.json"
-
+config_t = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_t.json"
+config_u = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_u.json"
+config_v = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_v.json"
+config_w = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_w.json"
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dom_nam, config=config_t)
 
@@ -29,7 +28,7 @@
 
 sci_v = coast.Gridded(dir_nam + fil_nam.replace("grid_T", "grid_V"), dom_nam, config=config_v)
 
-sci_v = coast.Nemo(dir_nam + fil_nam.replace("grid_T", "grid_V"), dom_nam, grid_ref="v-grid", multiple=False)
+# sci_v = coast.Nemo(dir_nam + fil_nam.replace("grid_T", "grid_V"), dom_nam, grid_ref="v-grid", multiple=False)
 
 # create an empty w-grid object, to store stratification
 sci_w = coast.Gridded(fn_domain=dom_nam, config=config_w)
diff --git a/example_scripts/export_to_netcdf_tutorial.py b/example_scripts/export_to_netcdf_tutorial.py
index f49a0717..95b66ea2 100644
--- a/example_scripts/export_to_netcdf_tutorial.py
+++ b/example_scripts/export_to_netcdf_tutorial.py
@@ -9,12 +9,15 @@
 import coast
 
 # And by defining some file paths
-fn_nemo_dat = "./example_files/COAsT_example_NEMO_data.nc"
-fn_nemo_dom = "./example_files/COAsT_example_NEMO_domain.nc"
+fn_nemo_dat = "./example_files/coast_example_nemo_data.nc"
+fn_nemo_dom = "./example_files/coast_example_nemo_domain.nc"
+config = "./config/example_nemo_grid_t.json"
+
 ofile = "example_export_output.nc"  # The target filename for output
 
 # We need to load in a NEMO object for doing NEMO things.
-nemo = coast.Nemo(fn_nemo_dat, fn_nemo_dom, grid_ref="t-grid")
+nemo = coast.Gridded(fn_nemo_dat, fn_nemo_dom, config=config)
+
 # We can export the whole xr.DataSet to a netCDF file
 nemo.dataset.to_netcdf(ofile, mode="w", format="NETCDF4")
 # Other file formats are available. From the documentation:
diff --git a/example_scripts/indexed_tutorial.py b/example_scripts/indexed_tutorial.py
new file mode 100644
index 00000000..37e5b68f
--- /dev/null
+++ b/example_scripts/indexed_tutorial.py
@@ -0,0 +1,55 @@
+"""
+This is a demonstration script for using INDEXED classes in the COAsT
+package. This object has strict data formatting requirements, which are
+outlined for various (observational) data types.
+"""
+
+import coast
+
+dir = "/projectsa/IMMERSE/WP8/TestData/"
+
+altimetry = coast.Altimetry(
+    file_path="./example_files/COAsT_example_altimetry_data.nc", config="./config/example_altimetry.json"
+)
+
+# Inspect object
+# $ altimetry.dataset
+
+profile = coast.Profile(file_path="./example_files/EN4_example.nc", config="./config/example_en4_profiles.json")
+
+# Inspect object
+# $ profile.dataset
+
+# GLider data https://linkedsystems.uk/erddap/files/Public_Glider_Data_0711/Doombar_20210214/Doombar_553_R.nc
+# More gliders downloadable at https://www.bodc.ac.uk/data/bodc_database/gliders/
+glider = coast.Glider(file_path=dir + "Bellatrix_555_R.nc", config="./config/example_glider_ego.json")
+
+# Inspect object
+# $ glider.dataset
+
+argos = coast.Argos(file_path=dir + "ARGOS.CSV", config="./config/example_argos.json")
+
+# Inspect object
+# $ argos.dataset
+
+oceanparcels = coast.Oceanparcels(
+    file_path=dir + "OceanParcels_GFDL-ESM2M.2043r3.1.nc", config="./config/example_oceanparcels.json"
+)
+
+# Inspect object
+# $ oceanparcels.dataset
+
+import datetime
+
+date0 = datetime.datetime(2007, 1, 10)
+date1 = datetime.datetime(2007, 1, 12)
+
+tidegauge = coast.Tidegauge(
+    file_path="./example_files/tide_gauges/lowestoft-p024-uk-bodc",
+    date_start=date0,
+    date_end=date1,
+    config="./config/example_tidegauge.json",
+)
+
+# Inspect object
+# $ tidegauge.dataset
diff --git a/example_scripts/internal_tide_pycnocline_diagnostics.py b/example_scripts/internal_tide_pycnocline_diagnostics.py
index 3f889527..987f3eed 100755
--- a/example_scripts/internal_tide_pycnocline_diagnostics.py
+++ b/example_scripts/internal_tide_pycnocline_diagnostics.py
@@ -24,8 +24,8 @@
     config = "AMM60"
     dir_AMM60 = "/projectsa/COAsT/NEMO_example_data/AMM60/"
     fil_nam_AMM60 = "AMM60_1d_20100704_20100708_grid_T.nc"
-    config_t = "/work/jelt/GitHub/COAsT/example_files/example_t_nemo_config.json"
-    config_w = "/work/jelt/GitHub/COAsT/example_files/example_w_nemo_config.json"
+    config_t = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_t.json"
+    config_w = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_w.json"
     mon = "July"
     # mon = 'Feb'
 
@@ -60,7 +60,9 @@
 
     dn_fig = "unit_testing/figures/"
     fn_nemo_grid_t_dat = "nemo_data_T_grid_Aug2015.nc"
-    fn_nemo_dom = "COAsT_example_NEMO_domain.nc"
+    fn_nemo_dom = "coast_example_nemo_domain.nc"
+    config_t = "config/example_nemo_grid_t.json"
+    config_w = "config/example_nemo_grid_w.json"
 
     sci_t = coast.Gridded(dn_files + fn_nemo_grid_t_dat, dn_files + fn_nemo_dom, config=config_t, multiple=True)
 
diff --git a/example_scripts/seasia_r12_example_plot.py b/example_scripts/seasia_r12_example_plot.py
index e407042a..36eefcb7 100644
--- a/example_scripts/seasia_r12_example_plot.py
+++ b/example_scripts/seasia_r12_example_plot.py
@@ -17,7 +17,7 @@
 dir_nam = "/projectsa/COAsT/NEMO_example_data/SEAsia_R12/"
 fil_nam = "SEAsia_R12_5d_20120101_20121231_gridT.nc"
 dom_nam = "domain_cfg_ORCA12_adj.nc"
-config_t = "/work/jelt/GitHub/COAsT/example_files/example_t_nemo_config.json"
+config_t = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_t.json"
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dir_nam + dom_nam, config=config_t)
 
diff --git a/example_scripts/tidegauge_tutorial.py b/example_scripts/tidegauge_tutorial.py
index 7bffd322..554c3b36 100755
--- a/example_scripts/tidegauge_tutorial.py
+++ b/example_scripts/tidegauge_tutorial.py
@@ -11,22 +11,26 @@
 import matplotlib.pyplot as plt
 
 # And by defining some file paths
-fn_nemo_dat = "./example_files/COAsT_example_NEMO_data.nc"
-fn_nemo_dom = "./example_files/COAsT_example_NEMO_domain.nc"
-fn_tidegauge = "./example_files/tide_gauges/lowestoft-p024-uk-bodc"
-fn_tidegauge_mult = "./example_files/tide_gauges/l*"
+dn_files = "./example_files/"
+fn_nemo_dat = dn_files + "coast_example_nemo_data.nc"
+fn_nemo_dom = dn_files + "coast_example_nemo_domain.nc"
+fn_nemo_config = "./config/example_nemo_grid_t.json"
+
+fn_tidegauge = dn_files + "tide_gauges/lowestoft-p024-uk-bodc"
+fn_tidegauge_mult = dn_files + "tide_gauges/l*"
+fn_tidegauge_config = "./config/example_tidegauge.json"
 
 # We need to load in a NEMO object for doing NEMO things.
-nemo = coast.Nemo(fn_nemo_dat, fn_nemo_dom, grid_ref="t-grid")
+nemo = coast.Gridded(fn_nemo_dat, fn_nemo_dom, config=fn_nemo_config)
 
 # And now we can load in our ALTIMETRY data. By default, TideGauge is set up
 # to read in GESLA ASCII files. However, if no path is supplied, then the
 # object's dataset will be initialised as None. Custom data can then be loaded
 # if desired, as long as it follows the data formatting for TideGauge. Here
 # we load data between two specified dates:
-date0 = datetime.datetime(2007, 1, 10)
-date1 = datetime.datetime(2007, 1, 16)
-tidegauge = coast.TideGauge(fn_tidegauge, date_start=date0, date_end=date1)
+date_0 = datetime.datetime(2007, 1, 10)
+date_1 = datetime.datetime(2007, 1, 16)
+tidegauge = coast.Tidegauge(file_path=fn_tidegauge, date_start=date_0, date_end=date_1, config=fn_tidegauge_config)
 
 # Before comparing our observations to the model, we will interpolate a model
 # variable to the same time and geographical space as the tidegauge. This is
diff --git a/example_scripts/tidetable_tutorial.py b/example_scripts/tidetable_tutorial.py
index 84cf6faa..372305c2 100755
--- a/example_scripts/tidetable_tutorial.py
+++ b/example_scripts/tidetable_tutorial.py
@@ -35,7 +35,7 @@
 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()
+tg = coast.Tidegauge()
 # specify the data read as a High Low Water dataset
 tg.dataset = tg.read_hlw_to_xarray(filnam, date_start, date_end)
 # Show dataset. If timezone is specified then it is presented as requested, otherwise uses UTC
diff --git a/example_scripts/wcssp_india_example_plot.py b/example_scripts/wcssp_india_example_plot.py
index 1d88445b..bbc712a6 100644
--- a/example_scripts/wcssp_india_example_plot.py
+++ b/example_scripts/wcssp_india_example_plot.py
@@ -18,7 +18,7 @@
 dir_nam = "/projectsa/COAsT/NEMO_example_data/MO_INDIA/"
 fil_nam = "ind_1d_cat_20180101_20180105_25hourm_grid_T.nc"
 dom_nam = "domain_cfg_wcssp.nc"
-config_t = "/work/jelt/GitHub/COAsT/example_files/example_t_nemo_config.json"
+config_t = "/work/jelt/GitHub/COAsT/config/example_nemo_grid_t.json"
 
 sci_t = coast.Gridded(dir_nam + fil_nam, dir_nam + dom_nam, config=config_t)
 


Test patch:
diff --git a/tests/config/__init__.py b/tests/config/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/config/config_parser_test.py b/tests/config/config_parser_test.py
index 047c389b..e8e4705d 100644
--- a/tests/config/config_parser_test.py
+++ b/tests/config/config_parser_test.py
@@ -20,7 +20,7 @@
     "grid_ref": {},
     "dimensionality": 3,
     "chunks": [],
-    "dataset": {"variable_map": {}, "dimension_map": {}},
+    "dataset": {"variable_map": {}, "dimension_map": {}, "coord_vars": [], "keep_all_vars": "False"},
     "domain": {
         "variable_map": {},
         "dimension_map": {},
@@ -34,7 +34,7 @@
     "type": "indexed",
     "dimensionality": 3,
     "chunks": [],
-    "dataset": {"variable_map": {}, "dimension_map": {}},
+    "dataset": {"variable_map": {}, "dimension_map": {}, "coord_vars": [], "keep_all_vars": "False"},
     "processing_flags": [],
 }
 
diff --git a/unit_testing/test_index_classes_load_datasets.py b/unit_testing/test_index_classes_load_datasets.py
new file mode 100644
index 00000000..785f52fb
--- /dev/null
+++ b/unit_testing/test_index_classes_load_datasets.py
@@ -0,0 +1,148 @@
+"""
+
+Test loading various types of INDEXED datasets with and without configuration
+ files
+
+        pytest -s tests/test_index_classes_load_datasets.py
+
+This file was in directory tests/ and executed under Git Actions. However, it
+calls external data that are too large for accessible hosting. Nevertheless,
+the structure proposed here "one file per type" could lend itself well to a
+restructure unit_testing work flow. So I move it there as a potential template
+for when the unit testing is restructured.
+"""
+from coast import Altimetry, Profile, Glider, Argos, Oceanparcels, Tidegauge
+import datetime
+
+# Altimetry data (NetCDF)
+fn_altimetry = "./example_files/COAsT_example_altimetry_data.nc"
+fn_altimetry_config = "config/example_altimetry.json"
+
+# En4 profile data (NetCDF)
+fn_profile = "./example_files/EN.4.2.1.f.profiles.l09.201501.nc"
+fn_profile_config = "config/example_en4_profiles.json"
+
+# EGO glider (NetCDF). More examples see https://www.bodc.ac.uk/data/bodc_database/gliders/
+fn_glider = "./example_files/Doombar_553_R.nc"
+fn_glider_config = "config/example_glider_ego.json"
+
+# Argos drifter data (CSV format)
+fn_argos = "./example_files/ARGOS.CSV"
+fn_argos_config = "config/example_argos.json"
+
+# Ocean parcels data (NetCDF)
+fn_oceanparc = "./example_files/OceanParcels_GFDL-ESM2M.2043r3.1.nc"
+fn_oceanparc_config = "config/example_oceanparcels.json"
+
+# Tidegauge data GESLA file (format version 3.0)
+fn_tidegauge = "./example_files/tide_gauges/portellen-p202-uk-bodc"
+fn_tidegauge_config = "config/example_tidegauge.json"
+fn_tidegauge_multiple = "./example_files/tide_gauges/m*"
+
+
+def test_load_altimetry_no_config():
+    altimetry = Altimetry(file_path=fn_altimetry)
+    assert altimetry is not None
+
+
+def test_load_altimetry_config():
+    altimetry = Altimetry(file_path=fn_altimetry, config=fn_altimetry_config)
+    assert altimetry is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(altimetry.dataset.coords) == ["time", "longitude", "latitude"]
+    assert list(altimetry.dataset.data_vars) == ["ocean_tide_standard_name"]
+
+
+def test_load_profile_no_config():
+    profile = Profile(file_path=fn_altimetry)
+    assert profile is not None
+
+
+def test_load_profile_config():
+    profile = Profile(file_path=fn_profile, config=fn_profile_config)
+    assert profile is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(profile.dataset.coords) == ["latitude", "longitude", "time"]
+    assert list(profile.dataset.data_vars) == [
+        "depth",
+        "temperature",
+        "qc_potential_temperature",
+        "qc_practical_salinity",
+        "qc_depth",
+        "qc_time",
+    ]
+
+
+def test_load_glider_no_config():
+    glider = Glider(file_path=fn_glider)
+    assert glider is not None
+
+
+def test_load_glider_config():
+    glider = Glider(file_path=fn_glider, config=fn_glider_config)
+    assert glider is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(glider.dataset.coords) == ["TIME", "TIME_GPS", "LATITUDE", "LONGITUDE", "PRES"]
+    assert list(glider.dataset.data_vars) == [
+        "PLATFORM_TYPE",
+        "PLATFORM_MAKER",
+        "LATITUDE_GPS",
+        "LONGITUDE_GPS",
+        "TEMP",
+    ]
+
+
+def test_load_argos_no_config():
+    argos = Argos(file_path=fn_argos)
+    assert argos is not None
+
+
+def test_load_argos_config():
+    argos = Argos(file_path=fn_argos, config=fn_argos_config)
+    assert argos is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(argos.dataset.coords) == ["time", "latitude", "longitude"]
+    assert list(argos.dataset.data_vars) == ["sea_surface_temperature", "pressure"]
+
+
+def test_load_oceanparcels_no_config():
+    oceanparc = Oceanparcels(file_path=fn_oceanparc)
+    assert oceanparc is not None
+
+
+def test_load_oceanparcels_config():
+    oceanparc = Oceanparcels(file_path=fn_oceanparc, config=fn_oceanparc_config)
+    assert oceanparc is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(oceanparc.dataset.coords) == ["time", "latitude", "longitude"]
+    assert list(oceanparc.dataset.data_vars) == ["depth", "t", "s"]
+
+
+def test_load_tidegauge_no_config():
+    tide_gauge = Tidegauge(file_path=fn_tidegauge)
+    assert tide_gauge is not None
+
+
+def test_load_tidegauge_config():
+    date_0 = datetime.datetime(2007, 1, 10)
+    date_1 = datetime.datetime(2007, 1, 12)
+
+    tide_gauge = Tidegauge(file_path=fn_tidegauge, date_start=date_0, date_end=date_1, config=fn_tidegauge_config)
+    assert tide_gauge is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(tide_gauge.dataset.coords) == ["time"]
+    assert list(tide_gauge.dataset.data_vars) == ["h", "qc_flags", "lat_dim", "lon_dim"]
+
+
+def test_load_tidegauge_multiple_config():
+    date_0 = datetime.datetime(2007, 1, 10)
+    date_1 = datetime.datetime(2007, 1, 12)
+
+    tide_gauge_list = Tidegauge.create_multiple_tidegauge(
+        file_list=fn_tidegauge_multiple, date_start=date_0, date_end=date_1, config=fn_tidegauge_config
+    )
+
+    assert tide_gauge_list is not None
+    # assert that we have the coordinate and data variable names as specified in the json config file
+    assert list(tide_gauge_list[0].dataset.coords) == ["time"]
+    assert list(tide_gauge_list[0].dataset.data_vars) == ["h", "qc_flags", "lat_dim", "lon_dim"]
diff --git a/unit_testing/unit_test.py b/unit_testing/unit_test.py
index a6b80d2b..82080f5b 100755
--- a/unit_testing/unit_test.py
+++ b/unit_testing/unit_test.py
@@ -77,15 +77,17 @@
 fn_nemo_grid_t_dat = "nemo_data_T_grid.nc"
 fn_nemo_grid_u_dat = "nemo_data_U_grid.nc"
 fn_nemo_grid_v_dat = "nemo_data_V_grid.nc"
-fn_nemo_dat = "COAsT_example_NEMO_data.nc"
-fn_nemo_dat_subset = "COAsT_example_NEMO_subset_data.nc"
-fn_nemo_dom = "COAsT_example_NEMO_domain.nc"
-fn_altimetry = "COAsT_example_altimetry_data.nc"
+fn_nemo_dat = "coast_example_nemo_data.nc"
+fn_nemo_dat_subset = "coast_example_nemo_subset_data.nc"
+fn_nemo_dom = "coast_example_nemo_domain.nc"
+fn_altimetry = "coast_example_altimetry_data.nc"
 fn_tidegauge = dn_files + "tide_gauges/lowestoft-p024-uk-bodc"
 fn_tidegauge2 = dn_files + "tide_gauges/LIV2010.txt"
-fn_EN4 = dn_files + "EN4_example.nc"
 fn_nemo_harmonics = "coast_nemo_harmonics.nc"
 fn_nemo_harmonics_dom = "coast_nemo_harmonics_dom.nc"
+# EN4 profile data (NetCDF)
+fn_profile = dn_files + "EN4_example.nc"
+fn_profile_config = "config/example_en4_profiles.json"
 
 sec = 1
 subsec = 96  # Code for '`' (1 below 'a')
@@ -137,12 +139,12 @@
     sci_load_file = coast.Nemo()
     sci_load_file.load(dn_files + fn_nemo_dat)
     if sci_load_ds.dataset.identical(sci_load_file.dataset):
-        print(str(sec) + chr(subsec) + " OK - COAsT.load_dataset()")
+        print(str(sec) + chr(subsec) + " OK - coast.load_dataset()")
     else:
         print(
             str(sec)
             + chr(subsec)
-            + " X - COAsT.load_dataset() ERROR - not identical to dataset loaded via COAsT.load()"
+            + " X - coast.load_dataset() ERROR - not identical to dataset loaded via coast.load()"
         )
 except:
     print(str(sec) + chr(subsec) + " FAILED")
@@ -331,14 +333,14 @@
 
 """
 #################################################
-## ( 2 ) Test general utility methods in COAsT ##
+## ( 2 ) Test general utility methods in coast ##
 #################################################
 """
 sec = sec + 1
 subsec = 96
 
 # -----------------------------------------------------------------------------#
-#%% ( 2a ) Copying a COAsT object                                               #
+#%% ( 2a ) Copying a coast object                                               #
 #                                                                             #
 
 subsec = subsec + 1
@@ -346,28 +348,28 @@
 try:
     sci_copy = sci.copy()
     if sci_copy.dataset == sci.dataset:
-        print(str(sec) + chr(subsec) + " OK - Copied COAsT object ")
+        print(str(sec) + chr(subsec) + " OK - Copied coast object ")
     else:
         print(str(sec) + chr(subsec) + " X - Copy Failed ")
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
 # -----------------------------------------------------------------------------#
-#%% ( 2b ) COAsT __getitem__ returns variable                                   #
+#%% ( 2b ) coast __getitem__ returns variable                                   #
 #                                                                             #
 
 subsec = subsec + 1
 
 try:
     if sci.dataset["ssh"].equals(sci["ssh"]):
-        print(str(sec) + chr(subsec) + " OK - COAsT.__getitem__ works correctly ")
+        print(str(sec) + chr(subsec) + " OK - coast.__getitem__ works correctly ")
     else:
-        print(str(sec) + chr(subsec) + " X - Problem with COAsT.__getitem__ ")
+        print(str(sec) + chr(subsec) + " X - Problem with coast.__getitem__ ")
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
 # -----------------------------------------------------------------------------#
-#%% ( 2c ) Renaming variables inside a COAsT object                             #
+#%% ( 2c ) Renaming variables inside a coast object                             #
 #                                                                             #
 
 subsec = subsec + 1
@@ -747,11 +749,11 @@
         print(
             str(sec)
             + chr(subsec)
-            + " OK - Nemo COAsT get_subset_as_xarray extracted expected array size and "
+            + " OK - Nemo coast get_subset_as_xarray extracted expected array size and "
             + "extreme values"
         )
     else:
-        print(str(sec) + chr(subsec) + " X - Issue with Nemo COAsT get_subset_as_xarray method")
+        print(str(sec) + chr(subsec) + " X - Issue with Nemo coast get_subset_as_xarray method")
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
@@ -794,7 +796,7 @@
 try:
     altimetry = coast.Altimetry(dn_files + fn_altimetry)
     ind = altimetry.subset_indices_lonlat_box([-10, 10], [45, 60])
-    altimetry_nwes = altimetry.isel(t_dim=ind)  # nwes = northwest europe shelf
+    altimetry_nwes = altimetry.isel(time=ind)  # nwes = northwest europe shelf
     ind_x, ind_y = general_utils.nearest_indices_2d(
         sci.dataset.longitude, sci.dataset.latitude, altimetry_nwes.dataset.longitude, altimetry_nwes.dataset.latitude
     )
@@ -859,9 +861,9 @@
 subsec = subsec + 1
 # We can load altimetry data straight from a CMEMS netcdf file on initialisation
 try:
-    altimetry = coast.Altimetry(dn_files + fn_altimetry)
+    altimetry = coast.Altimetry(dn_files + fn_altimetry, config="./config/example_altimetry.json")
 
-    # Test the data has loaded using attribute comparison, as for NEMO_data
+    # Test the data has loaded using attribute comparison, as for nemo_data
     alt_attrs_ref = dict(
         [
             ("source", "Jason-1 measurements"),
@@ -934,10 +936,10 @@
 # to the existing object. Here we look at the first option.
 
 try:
-    crps = altimetry_nwes.crps(sci, "ssh", "sla_filtered")
+    crps = altimetry_nwes.crps(sci, "ssh", "ocean_tide_standard_name")
 
     # TEST: Check length of crps and that it contains values
-    check1 = crps.dataset.crps.shape[0] == altimetry_nwes.dataset.sla_filtered.shape[0]
+    check1 = crps.dataset.crps.shape[0] == altimetry_nwes.dataset.ocean_tide_standard_name.shape[0]
     check2 = False in np.isnan(crps.dataset.crps)
     if check1 and check2:
         print(str(sec) + chr(subsec) + " OK - Altimetry CRPS")
@@ -958,13 +960,13 @@
 # Here we compare an altimetry variable to our interpolate model SSH
 
 try:
-    stats = altimetry_nwes.basic_stats("sla_filtered", "interp_ssh")
-    altimetry_nwes.basic_stats("sla_filtered", "interp_ssh", create_new_object=False)
+    stats = altimetry_nwes.basic_stats("ocean_tide_standard_name", "interp_ssh")
+    altimetry_nwes.basic_stats("ocean_tide_standard_name", "interp_ssh", create_new_object=False)
 
     # TEST: Check new object resembles internal object
     check1 = all(stats.dataset.error == altimetry_nwes.dataset.error)
     # TEST: Check lengths and values
-    check2 = stats.dataset.absolute_error.shape[0] == altimetry_nwes.dataset.sla_filtered.shape[0]
+    check2 = stats.dataset.absolute_error.shape[0] == altimetry_nwes.dataset.ocean_tide_standard_name.shape[0]
     if check1 and check2:
         print(str(sec) + chr(subsec) + " OK - Basic Stats for Altimetry")
     else:
@@ -1105,14 +1107,14 @@
     date_end = np.datetime64("now") - np.timedelta64(10, "D")
     eg = coast.TideGauge()
     # Extract the data between explicit dates
-    eg.dataset = eg.read_EA_API_to_xarray(date_start=date_start, date_end=date_end)
+    eg.dataset = eg.read_ea_api_to_xarray(date_start=date_start, date_end=date_end)
     check1 = eg.dataset.site_name == "Liverpool"
     check2 = len(eg.dataset.sea_level) > 0
     # eg.plot_timeseries()
 
     # Alternatively extract the data for the last ndays, here for a specific
     # (the default) station.
-    eg.dataset = eg.read_EA_API_to_xarray(ndays=1, stationId="E70124")
+    eg.dataset = eg.read_ea_api_to_xarray(n_days=1, station_id="E70124")
     check3 = eg.dataset.site_name == "Liverpool"
     # eg.plot_timeseries()
 
@@ -1248,7 +1250,7 @@
 try:
     date0 = datetime.datetime(2007, 1, 10)
     date1 = datetime.datetime(2007, 1, 12)
-    tidegauge_list = coast.TideGauge.create_multiple_tidegauge("./example_files/tide_gauges/l*", date0, date1)
+    tidegauge_list = coast.TideGauge.create_multiple("./example_files/tide_gauges/l*", date0, date1)
 
     # TEST: Check length of list
     check1 = len(tidegauge_list) == 2
@@ -1321,18 +1323,18 @@
 
     # Initiate a TideGauge object, if a filename is passed it assumes it is a GESLA type object
     tg = coast.TideGauge()
-    tg.dataset = tg.read_HLW_to_xarray(filnam, date_start, date_end)
+    tg.dataset = tg.read_hlw_to_xarray(filnam, date_start, date_end)
 
     check1 = len(tg.dataset.sea_level) == 37
-    check2 = tg.get_tidetabletimes(np.datetime64("2020-10-13 12:48"), method="nearest_HW").values == 8.01
-    check3 = tg.get_tidetabletimes(np.datetime64("2020-10-13 12:48"), method="nearest_1").time.values == np.datetime64(
-        "2020-10-13 14:36"
-    )
+    check2 = tg.get_tide_table_times(np.datetime64("2020-10-13 12:48"), method="nearest_HW").values == 8.01
+    check3 = tg.get_tide_table_times(
+        np.datetime64("2020-10-13 12:48"), method="nearest_1"
+    ).time.values == np.datetime64("2020-10-13 14:36")
     check4 = np.array_equal(
-        tg.get_tidetabletimes(np.datetime64("2020-10-13 12:48"), method="nearest_2").values, [2.83, 8.01]
+        tg.get_tide_table_times(np.datetime64("2020-10-13 12:48"), method="nearest_2").values, [2.83, 8.01]
     )
     check5 = np.array_equal(
-        tg.get_tidetabletimes(np.datetime64("2020-10-13 12:48"), method="window", winsize=24).values,
+        tg.get_tide_table_times(np.datetime64("2020-10-13 12:48"), method="window", winsize=24).values,
         [3.47, 7.78, 2.8, 8.01, 2.83, 8.45, 2.08, 8.71],
     )
 
@@ -1374,9 +1376,9 @@
     f.savefig(dn_fig + "tidegauge_optima.png")
 
     if check1 and check2 and check3 and check4:
-        print(str(sec) + chr(subsec) + " OK - Tidegauge local extrema found")
+        print(str(sec) + chr(subsec) + " OK - TideGauge local extrema found")
     else:
-        print(str(sec) + chr(subsec) + " X - Tidegauge local extrema")
+        print(str(sec) + chr(subsec) + " X - TideGauge local extrema")
 except:
     print(str(sec) + chr(subsec) + " FAILED.")
 
@@ -1416,9 +1418,9 @@
     f.savefig(dn_fig + "tidegauge_optima.png")
 
     if check1.all() and check2.all():
-        print(str(sec) + chr(subsec) + " OK - Tidegauge cubic extrema found")
+        print(str(sec) + chr(subsec) + " OK - TideGauge cubic extrema found")
     else:
-        print(str(sec) + chr(subsec) + " X - Tidegauge cubic extrema")
+        print(str(sec) + chr(subsec) + " X - TideGauge cubic extrema")
 except:
     print(str(sec) + chr(subsec) + " FAILED.")
 
@@ -1587,13 +1589,17 @@
 # Create Profile object and read EN4 example data file
 
 try:
-    profiles = coast.Profile()
-    profiles.read_EN4(fn_EN4)
+    # Create object without config file
+    profiles = coast.Profile(file_path=fn_profile)
+    check0 = profiles is not None
+
+    # Create object with config file
+    profiles = coast.Profile(file_path=fn_profile, config=fn_profile_config)
 
     # TEST: Check some data
     check1 = profiles.dataset.dims["z_dim"] == 400
     check2 = profiles.dataset.longitude[11].values == 9.89777
-    if check1 and check2:
+    if check0 and 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")
@@ -1928,6 +1934,7 @@
     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 internal_tide_pycnocline_diagnostics  # 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")
@@ -1953,10 +1960,10 @@
         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
+        # subsec = subsec + 1
+        # from example_scripts import internal_tide_pycnocline_diagnostics
 
-        print(str(sec) + chr(subsec) + " OK - tutorial on internal tides")
+        # print(str(sec) + chr(subsec) + " OK - tutorial on internal tides")
     else:
         print("Don't forget to test on a LIVLJOBS machine")
 
diff --git a/unit_testing/unit_test2.py b/unit_testing/unit_test2.py
index 02322aa9..ecd81f34 100644
--- a/unit_testing/unit_test2.py
+++ b/unit_testing/unit_test2.py
@@ -73,24 +73,26 @@
         print(f"location f{dn_files} cannot be found")
 
 dn_fig = "unit_testing/figures/"
-json_config_file_t_grid = path.join(dn_files, "example_t_nemo_config.json")
-json_config_file_f_grid = path.join(dn_files, "example_f_nemo_config.json")
-json_config_file_u_grid = path.join(dn_files, "example_u_nemo_config.json")
-json_config_file_v_grid = path.join(dn_files, "example_v_nemo_config.json")
-json_config_file_w_grid = path.join(dn_files, "example_w_nemo_config.json")
+json_config_file_t_grid = path.join("./config", "example_nemo_grid_t.json")
+json_config_file_f_grid = path.join("./config", "example_nemo_grid_f.json")
+json_config_file_u_grid = path.join("./config", "example_nemo_grid_u.json")
+json_config_file_v_grid = path.join("./config", "example_nemo_grid_v.json")
+json_config_file_w_grid = path.join("./config", "example_nemo_grid_w.json")
 fn_nemo_grid_t_dat_summer = "nemo_data_T_grid_Aug2015.nc"
 fn_nemo_grid_t_dat = "nemo_data_T_grid.nc"
 fn_nemo_grid_u_dat = "nemo_data_U_grid.nc"
 fn_nemo_grid_v_dat = "nemo_data_V_grid.nc"
-fn_nemo_dat = "Coast_example_NEMO_data.nc"
-fn_nemo_dat_subset = "Coast_example_NEMO_subset_data.nc"
-fn_nemo_dom = "Coast_example_NEMO_domain.nc"
-fn_altimetry = "Coast_example_altimetry_data.nc"
+fn_nemo_dat = "coast_example_nemo_data.nc"
+fn_nemo_dat_subset = "Coast_example_nemo_subset_data.nc"
+fn_nemo_dom = "coast_example_nemo_domain.nc"
+fn_altimetry = "coast_example_altimetry_data.nc"
 fn_tidegauge = dn_files + "tide_gauges/lowestoft-p024-uk-bodc"
 fn_tidegauge2 = dn_files + "tide_gauges/LIV2010.txt"
-fn_EN4 = dn_files + "EN4_example.nc"
 fn_nemo_harmonics = "coast_nemo_harmonics.nc"
 fn_nemo_harmonics_dom = "coast_nemo_harmonics_dom.nc"
+# EN4 profile data (NetCDF)
+fn_profile = dn_files + "EN4_example.nc"
+fn_profile_config = "config/example_en4_profiles.json"
 
 sec = 1
 subsec = 96  # Code for '`' (1 below 'a')
@@ -874,18 +876,19 @@
 subsec = 96
 # This section is for testing and demonstrating the use of the Altimetry
 # object. First begin by reloading Nemo t-grid test data:
-sci = coast.Gridded(dn_files + fn_nemo_dat, dn_files + fn_nemo_dom, config=json_config_file_t_grid)
+sci = coast.Nemo(dn_files + fn_nemo_dat, dn_files + fn_nemo_dom, grid_ref="t-grid")
+
 
 # -----------------------------------------------------------------------------#
-# %% ( 6a ) Load example altimetry data                                          #
+#%% ( 6a ) Load example altimetry data                                          #
 #                                                                             #
 
 subsec = subsec + 1
 # We can load altimetry data straight from a CMEMS netcdf file on initialisation
 try:
-    altimetry = coast.Altimetry(dn_files + fn_altimetry)
+    altimetry = coast.Altimetry(dn_files + fn_altimetry, config="./config/example_altimetry.json")
 
-    # Test the data has loaded using attribute comparison, as for NEMO_data
+    # Test the data has loaded using attribute comparison, as for nemo_data
     alt_attrs_ref = dict(
         [
             ("source", "Jason-1 measurements"),
@@ -903,8 +906,9 @@
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
+
 # -----------------------------------------------------------------------------#
-# %% ( 6b ) Altimetry subsetting                                                 #
+#%% ( 6b ) Altimetry subsetting                                                 #
 #                                                                             #
 
 subsec = subsec + 1
@@ -922,8 +926,9 @@
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
+
 # -----------------------------------------------------------------------------#
-# %% ( 6c ) Interpolate model to altimetry                                       #
+#%% ( 6c ) Interpolate model to altimetry                                       #
 #                                                                             #
 
 subsec = subsec + 1
@@ -944,8 +949,9 @@
 except:
     print(str(sec) + chr(subsec) + " FAILED")
 
+
 # -----------------------------------------------------------------------------#
-# %% ( 6d ) Altimetry CRPS                                                       #
+#%% ( 6d ) Altimetry CRPS                                                       #
 #                                                                             #
 
 
@@ -955,10 +961,10 @@
 # to the existing object. Here we look at the first option.
 
 try:
-    crps = altimetry_nwes.crps(sci, "ssh", "sla_filtered")
+    crps = altimetry_nwes.crps(sci, "ssh", "ocean_tide_standard_name")
 
     # TEST: Check length of crps and that it contains values
-    check1 = crps.dataset.crps.shape[0] == altimetry_nwes.dataset.sla_filtered.shape[0]
+    check1 = crps.dataset.crps.shape[0] == altimetry_nwes.dataset.ocean_tide_standard_name.shape[0]
     check2 = False in np.isnan(crps.dataset.crps)
     if check1 and check2:
         print(str(sec) + chr(subsec) + " OK - Altimetry CRPS")
@@ -968,8 +974,9 @@
 except:
     print(str(sec) + chr(subsec) + " FAILED.")
 
+
 # -----------------------------------------------------------------------------#
-# %% ( 6e ) Altimetry Stats methods                                              #
+#%% ( 6e ) Altimetry Stats methods                                              #
 #                                                                             #
 
 subsec = subsec + 1
@@ -978,13 +985,13 @@
 # Here we compare an altimetry variable to our interpolate model SSH
 
 try:
-    stats = altimetry_nwes.basic_stats("sla_filtered", "interp_ssh")
-    altimetry_nwes.basic_stats("sla_filtered", "interp_ssh", create_new_object=False)
+    stats = altimetry_nwes.basic_stats("ocean_tide_standard_name", "interp_ssh")
+    altimetry_nwes.basic_stats("ocean_tide_standard_name", "interp_ssh", create_new_object=False)
 
     # TEST: Check new object resembles internal object
     check1 = all(stats.dataset.error == altimetry_nwes.dataset.error)
     # TEST: Check lengths and values
-    check2 = stats.dataset.absolute_error.shape[0] == altimetry_nwes.dataset.sla_filtered.shape[0]
+    check2 = stats.dataset.absolute_error.shape[0] == altimetry_nwes.dataset.ocean_tide_standard_name.shape[0]
     if check1 and check2:
         print(str(sec) + chr(subsec) + " OK - Basic Stats for Altimetry")
     else:
@@ -993,8 +1000,9 @@
 except:
     print(str(sec) + chr(subsec) + " FAILED.")
 
+
 # -----------------------------------------------------------------------------#
-# %% ( 6f ) Altimetry quick_plot()                                               #
+#%% ( 6f ) Altimetry quick_plot()                                               #
 #                                                                             #
 
 subsec = subsec + 1
@@ -1011,6 +1019,7 @@
 
 plt.close("all")
 
+
 """
 #################################################
 ## ( 7 ) TideGauge Methods                     ##
@@ -1612,13 +1621,17 @@
 # Create Profile object and read EN4 example data file
 
 try:
-    profiles = coast.Profile()
-    profiles.read_EN4(fn_EN4)
+    # Create object without config file
+    profiles = coast.Profile(file_path=fn_profile)
+    check0 = profiles is not None
+
+    # Create object with config file
+    profiles = coast.Profile(file_path=fn_profile, config=fn_profile_config)
 
     # TEST: Check some data
     check1 = profiles.dataset.dims["z_dim"] == 400
     check2 = profiles.dataset.longitude[11].values == 9.89777
-    if check1 and check2:
+    if check0 and 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")