Merge pull request #41 from GeoscienceAustralia/rbt

Add latest changes to ensemble tide modelling code and upper intertidal cleanup

intertidal/exposure.py

@@ -26,7 +26,6 @@ def exposure(
     time_range,
     tide_model="FES2014",
     tide_model_dir="/var/share/tide_models",
-    ancillary_points="data/raw/corr_points.geojson"
 ):
     """
     Calculate exposure percentage for each pixel based on tide-height
@@ -42,11 +41,15 @@ def exposure(
         Tuple containing start and end time of time range to be used for
         tide model in the format of "YYYY-MM-DD".
     tide_model : str, optional
-        The tide model used to model tides, as supported by the `pyTMD`
-        Python package. Options include:
+        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)
-        - "TPXO8-atlas"
         - "TPXO9-atlas-v5"
+        - "TPXO8-atlas"
+        - "EOT20"
+        - "HAMTIDE11"
+        - "GOT4.10"
+        - "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
@@ -79,19 +82,32 @@ def exposure(
     pc_range = np.linspace(0, 1, 101)
 
     # Run the pixel_tides function with the calculate_quantiles option.
-    # For each pixel, an array of tideheights is returned, corresponding
+    # For each pixel, an array of tide heights is returned, corresponding
     # to the percentiles from pc_range of the timerange-tide model that
     # each tideheight appears in the model.
-    tide_cq, _ = pixel_tides_ensemble(
-        dem,
-        resample=True,
-        ancillary_points=ancillary_points,
-        calculate_quantiles=pc_range,
-        times=time_range,
-        model=tide_model,
-        directory=tide_model_dir,
-        cutoff=np.inf,
-    )
+    if (tide_model[0] == "ensemble") or (tide_model == "ensemble"):
+        # Use ensemble model combining multiple input ocean tide models
+        tide_cq, _ = pixel_tides_ensemble(
+            dem,
+            calculate_quantiles=pc_range,
+            times=time_range,
+            directory=tide_model_dir,
+            ancillary_points="data/raw/tide_correlations_2017-2019.geojson",
+            top_n=3,
+            reduce_method='mean',
+            resolution=3000,
+        )
+
+    else:
+        # Use single input ocean tide model
+        tide_cq, _ = pixel_tides(
+            dem,
+            calculate_quantiles=pc_range,
+            times=time_range,
+            resample=True,
+            model=tide_model,
+            directory=tide_model_dir,
+        )
 
     # Calculate the tide-height difference between the elevation value and
     # each percentile value per pixel

intertidal/extents.py

@@ -7,26 +7,25 @@
 
 from odc.algo import mask_cleanup
 from odc.geo.geom import Geometry
-import rioxarray
 import odc.geo.xr
 
 
 def load_reproject(
-    path, gbox, name=None, chunks={"x": 2048, "y": 2048}, **reproj_kwargs
+    path, gbox, chunks={"x": 2048, "y": 2048}, masked=True, **reproj_kwargs
 ):
     """
     Load and reproject part of a raster dataset into a given GeoBox.
     """
     ds = (
-        rioxarray.open_rasterio(
+        xr.open_dataset(
             path,
-            masked=True,
+            engine="rasterio",
+            masked=masked,
             chunks=chunks,
         )
         .squeeze("band")
         .odc.reproject(how=gbox, **reproj_kwargs)
     )
-    ds.name = name
 
     return ds
 
@@ -152,18 +151,49 @@ def extents(
     dc = datacube.Datacube(app="ocean_masking")
 
     # Load the land use dataset to mask out misclassified extents classes caused by urban land class
-    landuse_ds = load_reproject(
+    landuse_da = load_reproject(
         path=land_use_mask,
         gbox=dem.odc.geobox,
         resampling="nearest",
-    ).compute()
+    ).band_data.compute()
 
     # Separate out the 'intensive urban' land use summary class and set
     # all other pixels to False
-    reclassified = landuse_ds.isin(
-        [500, 530, 531, 532, 533, 534, 535, 536, 537, 538, 540, 541, 
-        550, 551, 552, 553, 554, 555, 560, 561, 562, 563, 564, 565, 
-        566, 567, 570, 571, 572, 573, 574, 575]
+    reclassified = landuse_da.isin(
+        [
+            500,
+            530,
+            531,
+            532,
+            533,
+            534,
+            535,
+            536,
+            537,
+            538,
+            540,
+            541,
+            550,
+            551,
+            552,
+            553,
+            554,
+            555,
+            560,
+            561,
+            562,
+            563,
+            564,
+            565,
+            566,
+            567,
+            570,
+            571,
+            572,
+            573,
+            574,
+            575,
+        ]
     )
 
     """--------------------------------------------------------------------"""

intertidal/tidal_bias_offset.py

@@ -4,91 +4,87 @@
 
 from dea_tools.spatial import subpixel_contours, points_on_line
 
-def bias_offset(
-                tide_m,
-                tide_cq,
-                extents,
-                lat_hat=True,
-                lot_hot=None
-                ):
+
+def bias_offset(tide_m, tide_cq, extents, lat_hat=True, lot_hot=None):
     """
-    Calculate the pixel-based sensor-observed spread and high/low 
-    offsets in tide heights compared to the full modelled tide range. 
-    Optionally, also return the highest and lowest astronomical and 
+    Calculate the pixel-based sensor-observed spread and high/low
+    offsets in tide heights compared to the full modelled tide range.
+    Optionally, also return the highest and lowest astronomical and
     sensor-observed tides for each pixel.
 
     Parameters
     ----------
     tide_m : xr.DataArray
-        An xarray.DataArray representing sensor observed tide heights 
-        for each pixel. Should have 'time', 'x' and 'y' in its 
+        An xarray.DataArray representing sensor observed tide heights
+        for each pixel. Should have 'time', 'x' and 'y' in its
         dimensions.
     tide_cq : xr.DataArray
-        An xarray.DataArray representing modelled tidal heights for 
-        each pixel. Should have 'quantile', 'x' and 'y' in its 
+        An xarray.DataArray representing modelled tidal heights for
+        each pixel. Should have 'quantile', 'x' and 'y' in its
         dimensions.
     extents : xr.DataArray
-        An xarray.DataArray representing 5 ecosystem class extents 
-        of the intertidal zone. Should have the same 
+        An xarray.DataArray representing 5 ecosystem class extents
+        of the intertidal zone. Should have the same
         dimensions as `tide_m` and `tide_cq`.
     lat_hat : bool, optional
-        Lowest/highest astronomical tides. This work considers the 
-        modelled tides to be equivalent to the astronomical tides. 
-        Default is True. 
+        Lowest/highest astronomical tides. This work considers the
+        modelled tides to be equivalent to the astronomical tides.
+        Default is True.
     lot_hot : bool, optional
         Lowest/highest sensor-observed tides. Default is None.
 
     Returns
     -------
-    Depending on the values of `lat_hat` and `lot_hot`, returns a tuple 
+    Depending on the values of `lat_hat` and `lot_hot`, returns a tuple
     with some or all of the following as xarray.DataArrays:
         * `lat`: The lowest astronomical tide.
         * `hat`: The highest astronomical tide.
         * `lot`: The lowest sensor-observed tide.
         * `hot`: The highest sensor-observed tide.
-        * `spread`: The spread of the observed tide heights as a 
+        * `spread`: The spread of the observed tide heights as a
         percentage of the modelled tide heights.
-        * `offset_lowtide`: The low tide offset measures the offset of the 
+        * `offset_lowtide`: The low tide offset measures the offset of the
         sensor-observed lowest tide from the minimum modelled tide.
-        * `offset_hightide`: The high tide measures the offset of the 
-        sensor-observed highest tide from the maximum modelled tide.         
+        * `offset_hightide`: The high tide measures the offset of the
+        sensor-observed highest tide from the maximum modelled tide.
     """
-    
-    # Set the maximum and minimum values per pixel for the observed and 
+
+    # Set the maximum and minimum values per pixel for the observed and
     # modelled datasets
-    max_obs = tide_m.max(dim='time')
-    min_obs = tide_m.min(dim='time')
-    max_mod = tide_cq.max(dim='quantile')
-    min_mod = tide_cq.min(dim='quantile')
-    
+    max_obs = tide_m.max(dim="time")
+    min_obs = tide_m.min(dim="time")
+    max_mod = tide_cq.max(dim="quantile")
+    min_mod = tide_cq.min(dim="quantile")
+
     # Set the maximum range in the modelled and observed tide heights
     mod_range = max_mod - min_mod
     obs_range = max_obs - min_obs
 
     # Calculate the spread of the observed tide heights as a percentage
     # of the modelled tide heights
-    spread = obs_range/mod_range * 100
-    
-    # Calculate the high and low tide offset of the observed tide 
+    spread = obs_range / mod_range * 100
+
+    # Calculate the high and low tide offset of the observed tide
     # heights as a percentage of the modelled highest and lowest tides.
-    offset_hightide = (abs(max_mod - max_obs))/mod_range * 100
-    offset_lowtide = (abs(min_mod - min_obs))/mod_range * 100
-        
+    offset_hightide = (abs(max_mod - max_obs)) / mod_range * 100
+    offset_lowtide = (abs(min_mod - min_obs)) / mod_range * 100
+
     # Add the lowest and highest astronomical tides
+    valid_mask = extents.isin([5, 4, 3])
     if lat_hat:
-        lat = min_mod.where((extents == 1)|(extents==2))
-        hat = max_mod.where((extents == 1)|(extents==2))
-    
+        lat = min_mod.where(valid_mask)
+        hat = max_mod.where(valid_mask)
+
     # Add the lowest and highest sensor-observed tides
     if lot_hot:
-        lot = min_obs.where((extents == 2)|(extents==1))
-        hot = max_obs.where((extents == 2)|(extents==1))
-        
+        lot = min_obs.where(valid_mask)
+        hot = max_obs.where(valid_mask)
+
     # Mask out non-intertidal pixels using ds extents
-    spread = spread.where((extents == 2)|(extents==1))
-    offset_hightide = offset_hightide.where((extents == 2)|(extents==1))
-    offset_lowtide = offset_lowtide.where((extents == 2)|(extents==1))
-    
+    spread = spread.where(valid_mask)
+    offset_hightide = offset_hightide.where(valid_mask)
+    offset_lowtide = offset_lowtide.where(valid_mask)
+
     if lat_hat:
         if lot_hot:
             return lat, hat, lot, hot, spread, offset_lowtide, offset_hightide
@@ -98,14 +94,11 @@ def bias_offset(
         return lot, hot, spread, offset_lowtide, offset_hightide
     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 
+
+
+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.
@@ -116,10 +109,10 @@ def tidal_offset_tidelines(extents,
         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 
+        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 
+        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
@@ -131,37 +124,40 @@ def tidal_offset_tidelines(extents,
     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.   
-    '''
+        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)
-    
+    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'])
-    
+    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'])
+    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')
+    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
-    
+    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
-