Merge pull request #150 from GeoscienceAustralia/external_tides

Support external tide data, upgrade DEA Tools

Author: erialC-P

intertidal/elevation.py

@@ -5,6 +5,7 @@
 import click
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import seaborn as sns
 import xarray as xr
 from dea_tools.dask import create_local_dask_cluster
@@ -744,6 +745,7 @@ def clean_edge_pixels(ds):
 
 def elevation(
     satellite_ds,
+    tide_data=None,
     valid_mask=None,
     ndwi_thresh=0.1,
     min_freq=0.01,
@@ -759,14 +761,18 @@ def elevation(
     log=None,
     **model_tides_kwargs,
 ):
-    """Generates DEA Intertidal Elevation outputs using satellite imagery
-    and tidal modeling.
+    """Generate DEA Intertidal Elevation outputs using satellite imagery and tide data.
 
     Parameters
     ----------
     satellite_ds : xarray.Dataset
         A satellite data time series containing an "ndwi" water index
         variable.
+    tide_data : array, optional
+        An optional array of tide heights or sea level measurements
+        matching the timesteps in `satellite_ds`. If this is provided,
+        these values will be used instead of modelling tides for each
+        satellite observation timestep.
     valid_mask : xr.DataArray, optional
         A boolean mask used to optionally constrain the analysis area,
         with the same spatial dimensions as `satellite_ds`. For example,
@@ -803,7 +809,8 @@ def elevation(
         determine workers.
     tide_model : str, optional
         The tide model or a list of models used to model tides, as
-        supported by the `eo-tides` Python package. Options include:
+        supported by the `eo-tides` Python package; ignored if `tide_data`
+        is provided. Options include:
         - "EOT20" (default)
         - "TPXO10-atlas-v2-nc"
         - "FES2022"
@@ -816,6 +823,7 @@ def elevation(
         The directory containing tide model data files. Defaults to
         "/var/share/tide_models"; for more information about the
         directory structure, refer to `eo-tides.utils.list_models`.
+        Ignored if `tide_data` is provided.
     run_id : string, optional
         An optional string giving the name of the analysis; used to
         prefix log entries.
@@ -852,22 +860,49 @@ def elevation(
     # Use run ID name for logs if it exists
     run_id = "Processing" if run_id is None else run_id
 
-    # Model tides into every pixel in the three-dimensional satellite
-    # 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(
-        data=satellite_ds,
-        model=tide_model,
-        directory=tide_model_dir,
-        **model_tides_kwargs,
-    )
+    # If tide heights are provided, use these directly after validating
+    # they have the current number of timesteps
+    if tide_data is not None:
+        log.info(f"{run_id}: Using provided tide heights")
+        # Convert to xarray if required
+        if isinstance(tide_data, pd.Series):
+            tide_data = tide_data.rename_axis("time").to_xarray()
+        elif isinstance(tide_data, xr.DataArray):
+            assert "time" in tide_data.dims, "Provided tide data must include a 'time' dimension."
+        else:
+            raise ValueError("Tide data must be provided in `pd.Series` or `xr.DataArray` format.")
+
+        # Verify that data includes the expected number of timesteps
+        if len(tide_data.time) != len(satellite_ds.time):
+            err_msg = (
+                "`tide_data` has a different number of timesteps "
+                f"({len(tide_data.time)}) than `satellite_ds` "
+                f"({len(satellite_ds.time)}). Ensure `tide_data` includes only a "
+                "single measurement for each satellite observation in `satellite_ds`."
+            )
+            raise ValueError(err_msg)
+
+        # Add to dataset
+        satellite_ds["tide_m"] = tide_data
+        tide_m = tide_data
+
+    # Otherwise, model tides into every pixel in the 3D (i.e. `x` by `y`
+    # by `time`) satellite dataset . If `model` is "ensemble" this will
+    # model tides by combining the best local tide models.
+    else:
+        log.info(f"{run_id}: Modelling tide heights for each pixel")
+        tide_m = pixel_tides(
+            data=satellite_ds,
+            model=tide_model,
+            directory=tide_model_dir,
+            **model_tides_kwargs,
+        )
 
-    # Set tide array pixels to nodata if the satellite data array pixels
-    # contain nodata. This ensures that we ignore any tide observations
-    # where we don't have matching satellite imagery
-    log.info(f"{run_id}: Masking nodata and adding tide heights to satellite data array")
-    satellite_ds["tide_m"] = tide_m.where(~satellite_ds.to_array().isel(variable=0).isnull().drop_vars("variable"))
+        # Set tide array pixels to nodata if the satellite data array pixels
+        # contain nodata. This ensures that we ignore any tide observations
+        # where we don't have matching satellite imagery
+        log.info(f"{run_id}: Masking nodata and adding tide heights to satellite data array")
+        satellite_ds["tide_m"] = tide_m.where(~satellite_ds.to_array().isel(variable=0).isnull().drop_vars("variable"))
 
     # Flatten array from 3D (time, y, x) to 2D (time, z) and drop pixels
     # with no correlation with tide. This greatly improves processing

notebooks/Getting_started_with_DEA_Intertidal.ipynb

@@ -14,7 +14,7 @@
     "\n",
     "<div class=\"alert alert-info\">\n",
     "    \n",
-    "**Note:** For more information on setting up tide models for DEA Intertidal using the `eo-tides` Python pacakge, refer to the [guide located here](https://geoscienceaustralia.github.io/eo-tides/setup/).\n",
+    "**Note:** For more information on setting up tide models for DEA Intertidal using the `eo-tides` Python package, refer to the [guide located here](https://geoscienceaustralia.github.io/eo-tides/setup/).\n",
     "\n",
     "</div>"
    ]
@@ -53,7 +53,7 @@
    "id": "c586a481-013f-4884-8e53-f3ddb9c438d7",
    "metadata": {},
    "source": [
-    "Install latest version of DEA Intertidal:"
+    "Install latest version of DEA Intertidal if required:"
    ]
   },
   {
@@ -73,7 +73,7 @@
     }
    ],
    "source": [
-    "%pip install -e . --quiet"
+    "# %pip install -e . --quiet"
    ]
   },
   {
@@ -1204,21 +1204,13 @@
     "fix, axes = plt.subplots(2, 2, figsize=(10, 8))\n",
     "axes = axes.flatten()\n",
     "\n",
-    "ds.elevation.plot.imshow(\n",
-    "    cmap=\"viridis\", robust=True, ax=axes[0], label=False, add_labels=False\n",
-    ")\n",
+    "ds.elevation.plot.imshow(cmap=\"viridis\", robust=True, ax=axes[0], label=False, add_labels=False)\n",
     "axes[0].set_title(\"Elevation (m)\")\n",
-    "ds.elevation_uncertainty.plot.imshow(\n",
-    "    cmap=\"inferno\", robust=True, ax=axes[1], label=False, add_labels=False\n",
-    ")\n",
+    "ds.elevation_uncertainty.plot.imshow(cmap=\"inferno\", robust=True, ax=axes[1], label=False, add_labels=False)\n",
     "axes[1].set_title(\"Elevation uncertainty (m)\")\n",
-    "ds.qa_ndwi_corr.plot.imshow(\n",
-    "    cmap=\"RdBu\", vmin=-0.7, vmax=0.7, ax=axes[2], label=False, add_labels=False\n",
-    ")\n",
+    "ds.qa_ndwi_corr.plot.imshow(cmap=\"RdBu\", vmin=-0.7, vmax=0.7, ax=axes[2], label=False, add_labels=False)\n",
     "axes[2].set_title(\"Tide correlation\")\n",
-    "ds.qa_ndwi_freq.plot.imshow(\n",
-    "    cmap=\"Blues\", vmin=0, vmax=1, ax=axes[3], label=False, add_labels=False\n",
-    ")\n",
+    "ds.qa_ndwi_freq.plot.imshow(cmap=\"Blues\", vmin=0, vmax=1, ax=axes[3], label=False, add_labels=False)\n",
     "axes[3].set_title(\"Inundation frequency (%)\");"
    ]
   },
@@ -1481,9 +1473,7 @@
     }
    ],
    "source": [
-    "ds.exposure.plot.imshow(\n",
-    "    cmap=cmocean.cm.matter_r, vmin=0, vmax=100, label=False, add_labels=False\n",
-    ")\n",
+    "ds.exposure.plot.imshow(cmap=cmocean.cm.matter_r, vmin=0, vmax=100, label=False, add_labels=False)\n",
     "plt.gca().set_title(\"Exposure\");"
    ]
   },
@@ -1561,105 +1551,7 @@
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {
-    "state": {
-     "04ff1dc786bd43ae9b40ac9e79054b76": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "ProgressStyleModel",
-      "state": {
-       "description_width": ""
-      }
-     },
-     "14cb61b551c5421d9c6787bb1059c213": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "HTMLModel",
-      "state": {
-       "layout": "IPY_MODEL_f7eb721ebfd84ab793c9ea035b513eee",
-       "style": "IPY_MODEL_e0202762aaa64a15bac48e8846e4efcc",
-       "value": "100%"
-      }
-     },
-     "76c1f85381964e93b4aa38fde6175eb8": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "HTMLModel",
-      "state": {
-       "layout": "IPY_MODEL_f9754946f72d4baf81112a1835d76a3d",
-       "style": "IPY_MODEL_c2532225aa694c5e9efdcbe107336c19",
-       "value": " 55/55 [01:04&lt;00:00,  1.16s/it]"
-      }
-     },
-     "84ab042d9b364ae8b79d2aae6f3d036e": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "FloatProgressModel",
-      "state": {
-       "bar_style": "success",
-       "layout": "IPY_MODEL_a059d536fe8d4aaf9fc3ce37fb0b3d77",
-       "max": 55,
-       "style": "IPY_MODEL_04ff1dc786bd43ae9b40ac9e79054b76",
-       "value": 55
-      }
-     },
-     "a059d536fe8d4aaf9fc3ce37fb0b3d77": {
-      "model_module": "@jupyter-widgets/base",
-      "model_module_version": "2.0.0",
-      "model_name": "LayoutModel",
-      "state": {}
-     },
-     "bdcdf2be0d65434486dc96907baca417": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "HBoxModel",
-      "state": {
-       "children": [
-        "IPY_MODEL_14cb61b551c5421d9c6787bb1059c213",
-        "IPY_MODEL_84ab042d9b364ae8b79d2aae6f3d036e",
-        "IPY_MODEL_76c1f85381964e93b4aa38fde6175eb8"
-       ],
-       "layout": "IPY_MODEL_f5bf1a0aa1344fafb10502093a0d2e95"
-      }
-     },
-     "c2532225aa694c5e9efdcbe107336c19": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "HTMLStyleModel",
-      "state": {
-       "description_width": "",
-       "font_size": null,
-       "text_color": null
-      }
-     },
-     "e0202762aaa64a15bac48e8846e4efcc": {
-      "model_module": "@jupyter-widgets/controls",
-      "model_module_version": "2.0.0",
-      "model_name": "HTMLStyleModel",
-      "state": {
-       "description_width": "",
-       "font_size": null,
-       "text_color": null
-      }
-     },
-     "f5bf1a0aa1344fafb10502093a0d2e95": {
-      "model_module": "@jupyter-widgets/base",
-      "model_module_version": "2.0.0",
-      "model_name": "LayoutModel",
-      "state": {}
-     },
-     "f7eb721ebfd84ab793c9ea035b513eee": {
-      "model_module": "@jupyter-widgets/base",
-      "model_module_version": "2.0.0",
-      "model_name": "LayoutModel",
-      "state": {}
-     },
-     "f9754946f72d4baf81112a1835d76a3d": {
-      "model_module": "@jupyter-widgets/base",
-      "model_module_version": "2.0.0",
-      "model_name": "LayoutModel",
-      "state": {}
-     }
-    },
+    "state": {},
     "version_major": 2,
     "version_minor": 0
    }