Merge pull request #99 from GeoscienceAustralia/eo_tides

Update DEA Intertidal to use `eo-tides`

Author: vnewey

intertidal/elevation.py

@@ -15,7 +15,7 @@
 from odc.geo.geom import BoundingBox
 from odc.algo import xr_quantile
 from datacube.utils.aws import configure_s3_access
-from dea_tools.coastal import pixel_tides
+from eo_tides.eo import pixel_tides
 from dea_tools.dask import create_local_dask_cluster
 
 from intertidal.io import (
@@ -29,6 +29,7 @@
 from intertidal.utils import (
     configure_logging,
     round_date_strings,
+    spearman_correlation,
 )
 from intertidal.extents import extents, load_connectivity_mask
 from intertidal.exposure import exposure
@@ -43,6 +44,7 @@ def ds_to_flat(
     max_freq=0.99,
     min_correlation=0.15,
     corr_method="pearson",
+    apply_threshold=True,
     correct_seasonality=False,
     valid_mask=None,
 ):
@@ -77,6 +79,11 @@ def ds_to_flat(
     corr_method : str, optional
         Correlation method to use. Defaults to "pearson", also supports
         "spearman".
+    apply_threshold : bool, optional
+        Whether to threshold the water index timeseries before calculating
+        correlations, to ensure small changes in index values beneath the
+        water surface are not included when calcualting correlations.
+        Defaults to True.
     correct_seasonality : bool, optional
         If True, remove any seasonal signal from the tide height data
         by subtracting monthly mean tide height from each value. This
@@ -132,12 +139,17 @@ def ds_to_flat(
     clear = clear.stack(z=("y", "x"))
 
     # Calculate correlations between NDWI water observations and tide
-    # height. Because we are only interested in pixels with inundation
-    # patterns (e.g.transitions from dry to wet) are driven by tide, we
-    # first convert NDWI into a boolean dry/wet layer before running the
-    # correlation. This prevents small changes in NDWI beneath the water
-    # surface from producing correlations with tide height.
-    wet_dry = flat_ds[index] > ndwi_thresh
+    # height. By default, because we are only interested in pixels with
+    # inundation patterns (e.g.transitions from dry to wet) that are driven
+    # by the tide, we first convert NDWI into a boolean dry/wet layer before
+    # running the correlation. This prevents small changes in NDWI beneath
+    # the water surface from producing correlations with tide height.
+    # This can be turned off by passing `apply_threshold=False`.
+    if apply_threshold:
+        wet_dry = flat_ds[index] > ndwi_thresh
+    else:
+        print("Using raw water index values for correlation")
+        wet_dry = flat_ds[index]
 
     # Use either tides directly or correct to remove seasonal signal
     if correct_seasonality:
@@ -149,19 +161,13 @@ def ds_to_flat(
 
     # Calculate correlation
     if corr_method == "pearson":
-        corr = xr.corr(wet_dry, tide_array, dim="time").rename("qa_ndwi_corr")
+        corr = xr.corr(wet_dry, tide_array, dim="time")
     elif corr_method == "spearman":
-        import xskillscore
-
-        corr = xskillscore.spearman_r(
-            flat_ds[index], tide_array, dim="time", skipna=True, keep_attrs=True
-        ).rename("qa_ndwi_corr")
-
-    # TODO: investigate alternative function from DEA Tools
-    # (doesn't currently handle multiple tide models)
-    # corr = lag_linregress_3D(x=flat_ds.tide_m, y=wet_dry).cor.rename("qa_ndwi_corr")
+        print("Applying Spearman correlation")
+        corr = spearman_correlation(x=wet_dry, y=tide_array, dim="time")
 
     # Keep only pixels with correlations that meet min threshold
+    corr = corr.rename("qa_ndwi_corr")
     corr_mask = corr >= min_correlation
     flat_ds = flat_ds.where(corr_mask, drop=True)
 
@@ -785,14 +791,14 @@ def elevation(
     window_prop_tide=0.15,
     correct_seasonality=False,
     max_workers=None,
-    tide_model="FES2014",
+    tide_model="EOT20",
     tide_model_dir="/var/share/tide_models",
     run_id=None,
     log=None,
 ):
     """
-    Calculates DEA Intertidal Elevation using satellite imagery and
-    tidal modeling.
+    Generates DEA Intertidal Elevation outputs using satellite imagery
+    and tidal modeling.
 
     Parameters
     ----------
@@ -835,18 +841,19 @@ def elevation(
         determine workers.
     tide_model : str, optional
         The tide model or a list of models used to model tides, as
-        supported by the `pyTMD` Python package. Options include:
-        - "FES2014" (default; pre-configured on DEA Sandbox)
-        - "TPXO9-atlas-v5"
-        - "TPXO8-atlas"
-        - "EOT20"
-        - "HAMTIDE11"
-        - "GOT4.10"
+        supported by the `eo-tides` Python package. Options include:
+        - "EOT20" (default)
+        - "TPXO10-atlas-v2-nc"
+        - "FES2022"
+        - "FES2022_extrapolated"
+        - "FES2014"
+        - "FES2014_extrapolated"
+        - "GOT5.6"
         - "ensemble" (experimental: combine all above into single ensemble)
     tide_model_dir : str, optional
         The directory containing tide model data files. Defaults to
         "/var/share/tide_models"; for more information about the
-        directory structure, refer to `dea_tools.coastal.model_tides`.
+        directory structure, refer to `eo-tides.utils.list_models`.
     run_id : string, optional
         An optional string giving the name of the analysis; used to
         prefix log entries.
@@ -881,8 +888,8 @@ def elevation(
     # dataset (x by y by time). If `model` is "ensemble" this will model
     # tides by combining the best local tide models.
     log.info(f"{run_id}: Modelling tide heights for each pixel")
-    tide_m, _ = pixel_tides(
-        ds=satellite_ds,
+    tide_m = pixel_tides(
+        data=satellite_ds,
         model=tide_model,
         directory=tide_model_dir,
     )
@@ -1097,18 +1104,18 @@ def elevation(
     "--tide_model",
     type=str,
     multiple=True,
-    default=["FES2014"],
+    default=["EOT20"],
     help="The model used for tide modelling, as supported by the "
-    "`pyTMD` Python package. Options include 'FES2014' (default), "
-    "'TPXO9-atlas-v5', 'TPXO8-atlas-v1', 'EOT20', 'HAMTIDE11', 'GOT4.10'. ",
+    "`eo-tides` Python package. Options include 'EOT20' (default), "
+    "'TPXO10-atlas-v2-nc', 'FES2022', 'FES2014', 'GOT5.6', 'ensemble'."
 )
 @click.option(
     "--tide_model_dir",
     type=str,
     default="/var/share/tide_models",
     help="The directory containing tide model data files. Defaults to "
     "'/var/share/tide_models'; for more information about the required "
-    "directory structure, refer to `dea_tools.coastal.model_tides`.",
+    "directory structure, refer to `eo-tides.utils.list_models`.",
 )
 @click.option(
     "--modelled_freq",

intertidal/exposure.py

@@ -9,7 +9,7 @@
 import pandas as pd
 
 from math import ceil
-from dea_tools.coastal import _pixel_tides_resample, pixel_tides
+from eo_tides.eo import _pixel_tides_resample, pixel_tides
 from intertidal.utils import configure_logging, round_date_strings
 
 
@@ -235,7 +235,7 @@ def exposure(
     start_date,
     end_date,
     modelled_freq="30min",
-    tide_model="FES2014",
+    tide_model="EOT20",
     tide_model_dir="/var/share/tide_models",
     filters=None,
     filters_combined=None,
@@ -281,19 +281,19 @@ def exposure(
         cadence. Defaults to '30min'.
     tide_model : str, optional
         The tide model or a list of models used to model tides, as
-        supported by the `pyTMD` Python package. Options include:
-        - "FES2014" (default; pre-configured on DEA Sandbox)
+        supported by the `eo-tides` Python package. Options include:
+        - "EOT20" (default)
+        - "TPXO10-atlas-v2-nc"
         - "FES2022"
-        - "TPXO9-atlas-v5"
-        - "TPXO8-atlas"
-        - "EOT20"
-        - "HAMTIDE11"
-        - "GOT4.10"
-        - "ensemble" (combine above into single ensemble)
+        - "FES2022_extrapolated"
+        - "FES2014"
+        - "FES2014_extrapolated"
+        - "GOT5.6"
+        - "ensemble" (experimental: combine all above into single ensemble)
     tide_model_dir : str, optional
         The directory containing tide model data files. Defaults to
         "/var/share/tide_models"; for more information about the
-        directory structure, refer to `dea_tools.coastal.model_tides`.
+        directory structure, refer to `eo-tides.utils.list_models`.
     filters : list of strings, optional
         An optional list of customisation options to input into the tidal
         modelling to calculate exposure. Filters include:
@@ -441,9 +441,9 @@ def exposure(
 
     # Run tide model at low resolution
     modelledtides_lowres = pixel_tides(
-        ds=dem,
+        data=dem,
+        time=time_range,
         model=tide_model,
-        times=time_range,
         directory=tide_model_dir,
         resample=False,
     )
@@ -455,17 +455,18 @@ def exposure(
     # pixel resolution if any filter is "unfiltered"
     if "unfiltered" in filters:
 
-        # Convert to quantiles
-        modelledtides_lowres_quantiles = modelledtides_lowres.quantile(
-            q=calculate_quantiles, dim="time"
-        ).astype(modelledtides_lowres.dtype)
-
-        # Reproject into pixel resolution, after making sure CRS is present
-        modelledtides_highres, _ = _pixel_tides_resample(
-            tides_lowres=modelledtides_lowres_quantiles.odc.assign_crs(
-                dem.odc.geobox.crs
-            ),
-            ds=dem,
+        # Convert to quantiles, and make sure CRS is present
+        modelledtides_lowres_quantiles = (
+            modelledtides_lowres.quantile(q=calculate_quantiles, dim="time")
+            .astype(modelledtides_lowres.dtype)
+            .odc.assign_crs(dem.odc.geobox.crs)
+        )
+
+        # Reproject into pixel resolution
+        modelledtides_highres = _pixel_tides_resample(
+            tides_lowres=modelledtides_lowres_quantiles,
+            dask_chunks=dem.shape,
+            gbox=dem.odc.geobox,
         )
 
         # Add pixel resolution tides into to output dataset

intertidal/tidal_bias_offset.py

@@ -11,6 +11,8 @@ def bias_offset(tide_m, tide_cq, lat_hat=True, lot_hot=None):
     Optionally, also return the highest and lowest astronomical and
     sensor-observed tides for each pixel.
 
+    TODO: update to use `eo-tides.stats` functionality
+
     Parameters
     ----------
     tide_m : xr.DataArray
@@ -84,69 +86,3 @@ def bias_offset(tide_m, tide_cq, lat_hat=True, lot_hot=None):
     else:
         return spread, offset_lowtide, offset_hightide
 
-
-# def tidal_offset_tidelines(extents, offset_hightide, offset_lowtide, distance=500):
-#     """
-#     This function extracts high and low tidelines from a rasterised
-#     'sometimes wet' layer in the extents input xr.DataArray,
-#     calculates the tidal offsets at each point on the lines, and returns
-#     the offset values in separate `geopandas.GeoDataFrame` objects.
-
-#     Parameters
-#     ----------
-#     extents : xarray.DataArray
-#         An xarray.DataArray containing binary shoreline information,
-#         depicting always, sometimes and never wet pixels.
-#     offset_hightide: xarray.DataArray
-#         An xarray.DataArray containing the percentage high-tide offset of the
-#         satellite observed tide heights from the modelled heights.
-#     offset_lowtide: xarray.DataArray
-#         An xarray.DataArray containing the percentage low-tide offset of the
-#         satellite observed tide heights from the modelled heights.
-#     distance : integer or float, optional
-#         A number giving the interval at which to generate points along
-#         the line feature. Defaults to 500, which will generate a point
-#         at every 500 metres along the line.
-
-#     Returns
-#     -------
-#     tuple
-#         A tuple of two `geopandas.GeoDataFrame` objects containing the
-#         high and low tidelines with their respective tidal offsets and
-#         a `geopandas.GeoDataFrame` containing the multilinestring tidelines.
-#     """
-#     ## Create a three class extents dataset: tidal wet/intertidal/dry
-#     extents = extents.where((extents == 0) | (extents == 1) | (extents == 2), np.nan)
-
-#     # Extract the high/low tide boundaries
-#     tidelines_gdf = subpixel_contours(da=extents, z_values=[1.5, 0.5])
-
-#     # Translate the high/Low tidelines into point data at regular intervals
-#     lowtideline = points_on_line(tidelines_gdf, 0, distance=distance)
-#     hightideline = points_on_line(tidelines_gdf, 1, distance=distance)
-
-#     # Extract the point coordinates into xarray for the hightideline dataset
-#     x_indexer_high = xr.DataArray(hightideline.centroid.x, dims=["point"])
-#     y_indexer_high = xr.DataArray(hightideline.centroid.y, dims=["point"])
-
-#     # Extract the point coordinates into xarray for the lowtideline dataset
-#     x_indexer_low = xr.DataArray(lowtideline.centroid.x, dims=["point"])
-#     y_indexer_low = xr.DataArray(lowtideline.centroid.y, dims=["point"])
-
-#     # Extract the high or low tide offset at each point in the high and low tidelines respectively
-#     # From https://stackoverflow.com/questions/67425567/extract-values-from-xarray-dataset-using-geopandas-multilinestring
-#     highlineoffset = offset_hightide.sel(
-#         x=x_indexer_high, y=y_indexer_high, method="nearest"
-#     )
-#     lowlineoffset = offset_lowtide.sel(
-#         x=x_indexer_low, y=y_indexer_low, method="nearest"
-#     )
-
-#     # Replace the offset values per point into the master dataframes
-#     hightideline["offset_hightide"] = highlineoffset
-#     lowtideline["offset_lowtide"] = lowlineoffset
-
-#     ## Consider adding the points to the tidelines
-#     ## https://gis.stackexchange.com/questions/448788/merging-points-to-linestrings-using-geopandas
-
-#     return hightideline, lowtideline, tidelines_gdf