Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extents updates #91

Merged
merged 14 commits into from
Sep 24, 2024
Merged

Extents updates #91

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
adb9048
Integrate new tide modelling, connectivity
robbibt Sep 20, 2024
5c348a5
Update all, bugfinding
robbibt Sep 20, 2024
57b443c
Update all, bugfinding
robbibt Sep 20, 2024
203a9fe
Add latest versions
robbibt Sep 22, 2024
848bdba
Fix import
robbibt Sep 23, 2024
cdca3dd
Update packaging
robbibt Sep 23, 2024
3826470
Update packaging
robbibt Sep 23, 2024
07930bf
Commit to repo to enable PR
erialC-P Sep 23, 2024
169c402
manually creating new requirements.txt file
erialC-P Sep 23, 2024
4c98007
Update nodata value
robbibt Sep 23, 2024
aef8bfd
Fix tests
robbibt Sep 23, 2024
3738e36
Automatically update integration test validation results
robbibt Sep 24, 2024
b957436
Update notebooks and docstrings
robbibt Sep 24, 2024
c7f4ca3
Automatically update integration test validation results
robbibt Sep 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 34 additions & 80 deletions intertidal/elevation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@
load_data,
load_topobathy_mask,
load_aclum_mask,
load_ocean_mask,
prepare_for_export,
tidal_metadata,
export_dataset_metadata,
Expand All @@ -31,46 +30,10 @@
configure_logging,
round_date_strings,
)
from intertidal.tide_modelling import pixel_tides_ensemble
from intertidal.extents import extents#, ocean_connection
from intertidal.extents import extents, load_connectivity_mask
from intertidal.exposure import exposure
from intertidal.tidal_bias_offset import bias_offset

def ocean_connection(water, ocean_da, connectivity=2):
"""
Identifies areas of water pixels that are adjacent to or directly
connected to intertidal pixels.

Parameters:
-----------
water : xarray.DataArray
An array containing True for water pixels.
ocean_da : xarray.DataArray
An array containing True for ocean pixels.
connectivity : integer, optional
An integer passed to the 'connectivity' parameter of the
`skimage.measure.label` function.

Returns:
--------
ocean_connection : xarray.DataArray
An array containing the a mask consisting of identified
ocean-connected pixels as True.
"""

# First, break `water` array into unique, discrete
# regions/blobs.
blobs = xr.apply_ufunc(label, water, 0, False, connectivity)

# For each unique region/blob, use region properties to determine
# whether it overlaps with a feature from `intertidal`. If
# it does, then it is considered to be adjacent or directly connected
# to intertidal pixels
ocean_connection = blobs.isin(
[i.label for i in regionprops(blobs.values, ocean_da.values) if i.max_intensity]
)

return ocean_connection

def ds_to_flat(
satellite_ds,
Expand Down Expand Up @@ -138,27 +101,31 @@ def ds_to_flat(
corr : xr.DataArray
Correlation of NDWI pixel wetness with tide height.
"""

# If an overall valid data mask is provided, apply to the data first
if valid_mask is not None:
satellite_ds = satellite_ds.where(valid_mask)

# Flatten satellite dataset by stacking "y" and "x" dimensions, then
# drop any pixels that are empty across all-of-time
flat_ds = satellite_ds.stack(z=("y", "x")).dropna(dim="time", how="all")

# Calculate frequency of wet per pixel, then threshold
# to exclude always wet and always dry
freq = (
(flat_ds[index] > ndwi_thresh)
.where(~flat_ds[index].isnull())
(satellite_ds[index] > ndwi_thresh)
.where(~satellite_ds[index].isnull())
.mean(dim="time")
.rename("qa_ndwi_freq")
)

# Mask out pixels outside of frequency bounds
freq_mask = (freq >= min_freq) & (freq <= max_freq)
satellite_ds = satellite_ds.where(freq_mask)

# Flatten to 1D, dropping any pixels that are not in frequency mask
flat_ds = flat_ds.where(freq_mask, drop=True)
# If an overall valid data mask is provided, apply to the data
if valid_mask is not None:
satellite_ds = satellite_ds.where(valid_mask)

# Flatten satellite and freq data by stacking "y" and "x" dims.
# Drop any pixels that are always empty, or empty timesteps
flat_ds = (
satellite_ds.stack(z=("y", "x"))
.dropna(dim="time", how="all")
.dropna(dim="z", how="all")
)
freq = freq.stack(z=("y", "x"))

# Calculate correlations between NDWI water observations and tide
# height. Because we are only interested in pixels with inundation
Expand All @@ -176,6 +143,7 @@ def ds_to_flat(
else:
tide_array = flat_ds.tide_m

# Calculate correlation
if corr_method == "pearson":
corr = xr.corr(wet_dry, tide_array, dim="time").rename("qa_ndwi_corr")
elif corr_method == "spearman":
Expand Down Expand Up @@ -805,7 +773,6 @@ def clean_edge_pixels(ds):
def elevation(
satellite_ds,
valid_mask=None,
ocean_mask=None,
ndwi_thresh=0.1,
min_freq=0.01,
max_freq=0.99,
Expand Down Expand Up @@ -834,12 +801,6 @@ def elevation(
this could be a mask generated from a topo-bathy DEM, used to
limit the analysis to likely intertidal pixels. Default is None,
which will not apply a mask.
ocean_mask : xr.DataArray, optional
An optional mask identifying ocean pixels within the analysis
area, with the same spatial dimensions as `satellite_ds`.
If provided, this will be used to restrict the analysis to pixels
that are directly connected to ocean waters. Defaults is None,
which will not apply a mask.
ndwi_thresh : float, optional
A threshold value for the normalized difference water index
(NDWI) above which pixels are considered water, by default 0.1.
Expand Down Expand Up @@ -916,13 +877,14 @@ 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_ensemble(
tide_m, _ = pixel_tides(
ds=satellite_ds,
ancillary_points="data/raw/tide_correlations_2017-2019.geojson",
model=tide_model,
directory=tide_model_dir,
ranking_points="data/raw/tide_correlations_2017-2019.geojson",
ensemble_top_n=3,
)

# 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
Expand Down Expand Up @@ -999,16 +961,6 @@ def elevation(
elevation_bands = [d for d in ds.data_vars if "elevation" in d]
ds[elevation_bands] = clean_edge_pixels(ds[elevation_bands])

# Mask out any non-ocean connected elevation pixels.
# `~(ds.qa_ndwi_freq < min_freq)` ensures that nodata pixels are
# treated as wet
if ocean_mask is not None:
log.info(f"{run_id}: Restricting outputs to ocean-connected waters")
ocean_connected_mask = ocean_connection(
~(ds.qa_ndwi_freq < min_freq), ocean_mask
)
ds[elevation_bands] = ds[elevation_bands].where(ocean_connected_mask)

# Return output data and tide height array
log.info(f"{run_id}: Successfully completed intertidal elevation modelling")
return ds, tide_m
Expand Down Expand Up @@ -1252,11 +1204,11 @@ def intertidal_cli(
satellite_ds.load()

# Load topobathy mask from GA's AusBathyTopo 250m 2023 Grid,
# urban land use class mask from ABARES CLUM, and ocean mask
# from geodata_coast_100k
topobathy_mask = load_topobathy_mask(dc, satellite_ds.odc.geobox.compat)
reclassified_aclum = load_aclum_mask(dc, satellite_ds.odc.geobox.compat)
ocean_mask = load_ocean_mask(dc, satellite_ds.odc.geobox.compat)
# urban land use class mask from ABARES CLUM, and coastal mask
# from least-cost connectivity analysis
topobathy_mask = load_topobathy_mask(dc, satellite_ds.odc.geobox)
urban_mask = load_aclum_mask(dc, satellite_ds.odc.geobox)
coastal_mask, _ = load_connectivity_mask(dc, satellite_ds.odc.geobox)

# Also load ancillary dataset IDs to use in metadata
# (both layers are continental continental products with only
Expand All @@ -1269,8 +1221,7 @@ def intertidal_cli(
log.info(f"{run_id}: Calculating Intertidal Elevation")
ds, tide_m = elevation(
satellite_ds,
valid_mask=topobathy_mask,
# ocean_mask=ocean_mask,
valid_mask=topobathy_mask & coastal_mask,
ndwi_thresh=ndwi_thresh,
min_freq=min_freq,
max_freq=max_freq,
Expand All @@ -1287,8 +1238,11 @@ def intertidal_cli(
# Calculate extents (to be included in next version)
log.info(f"{run_id}: Calculating Intertidal Extents")
ds["extents"] = extents(
dc=dc,
ds=ds
dem=ds.elevation,
freq=ds.qa_ndwi_freq,
corr=ds.qa_ndwi_corr,
urban_mask=urban_mask,
coastal_mask=coastal_mask,
)

if exposure_offsets:
Expand Down
10 changes: 5 additions & 5 deletions intertidal/exposure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,7 @@
import pandas as pd

from math import ceil
from dea_tools.coastal import _pixel_tides_resample
from intertidal.tide_modelling import pixel_tides_ensemble
from dea_tools.coastal import _pixel_tides_resample, pixel_tides
from intertidal.utils import configure_logging, round_date_strings


Expand Down Expand Up @@ -441,12 +440,13 @@ def exposure(
), f'Nominated filter "{x}" is not in {all_filters}. Check spelling and retry'

# Run tide model at low resolution
modelledtides_lowres = pixel_tides_ensemble(
dem,
modelledtides_lowres = pixel_tides(
ds=dem,
model=tide_model,
times=time_range,
directory=tide_model_dir,
ancillary_points="data/raw/tide_correlations_2017-2019.geojson",
ranking_points="data/raw/tide_correlations_2017-2019.geojson",
ensemble_top_n=3,
resample=False,
)

Expand Down
93 changes: 48 additions & 45 deletions intertidal/extents.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,7 @@
from dea_tools.spatial import xr_rasterize
from rasterio.features import sieve
from intertidal.io import (
load_ocean_mask,
load_aclum_mask,
load_topobathy_mask,
extract_geobox,
)

Expand Down Expand Up @@ -251,43 +249,56 @@ def load_gmw_mask(ds, gmw_path="/gdata1/data/mangroves/gmw_v3_2007_vec_aus.geojs
return gmw_da


def extents(dc, ds, buffer=20000, min_correlation=0.15, sieve_size=5):





# TODO: pass in ACLUM and connectivity mask directly


def extents(
dem, freq, corr, urban_mask, coastal_mask, buffer=20000, min_correlation=0.15, sieve_size=5, **connectivity_kwargs
):
"""
Classify coastal ecosystems into broad classes based
on their respective patterns of wetting frequency,
proximity to the ocean, and relationships to tidal
inundation, elevation and urban land use (to mask misclassifications).
inundation, elevation and urban land use (to mask
misclassifications).

Parameters:
-----------
dc : Datacube
A Datacube instance for loading data.
ds : xarray.Dataset
ds : xarray.Dataset
An xarray.Dataset that must include xarray.DataArrays for
at least 'elevation', 'qa_ndwi_freq' and 'qa_ndwi_corr'.
These arrays are generated as an output from the
intertidal.elevation algorithm.
`intertidal.elevation` algorithm.
buffer : int, optional
The distance by which to buffer the ds GeoBox to reduce edge
The distance by which to buffer the `ds` GeoBox to reduce edge
effects. This buffer will eventually be removed and clipped back
to the original GeoBox extent. Defaults to 20,000 metres.
min_correlation : int, optional
min_correlation : int, optional
Minimum correlation between water index and tide height
required for a pixel to be included in the intertidal pixel
analysis, 0.15 by default, aligning with the default value
used in the intertidal.elevation algorithm.
sieve_size : int, optional
used in the `intertidal.elevation` algorithm.
sieve_size : int, optional
Maximum number of grouped pixels belonging to any single
class that are sieved out to remove small noisy dataset
features. Class values are replaced with the dominant
neighbouring class.
**connectivity_kwargs :
Optional keyword arguments to pass to `load_connectivity_mask`.

Returns:
--------
extents: xarray.DataArray
A categorical xarray.DataArray depicting the following pixel classes:
- Nodata (0),
- Ocean and coastal water (1),
- Nodata (255),
- Ocean and coastal waters (1),
- Exposed intertidal (low confidence) (2),
- Exposed intertidal (high confidence) (3),
- Inland waters (4),
Expand All @@ -297,57 +308,49 @@ class that are sieved out to remove small noisy dataset
------
Classes are defined as follows:

0 : Nodata
1 : Ocean and coastal water
Pixels that are wet in 50 % or more observations
2 : Exposed intertidal (low confidence)
Pixels with water index and tidal correlations higher than 0.15.
Pixels must be located within the costdist_mask connectivity mask
and not be included in the intertidal elevation (high confidence)
dataset
3 : Exposed intertidal (high confidence)
Pixels that are included in the intertidal elevation dataset
4 : Inland waters
Pixels that are wet in more than 50 % of observations and fall
outside of the costdist_mask connectivity mask.
5 : Land
Pixels that are wet in less than 50 % of observations
255: Nodata
1: Ocean and coastal waters
Pixels that are wet in 50% or more observations and located within
the coastal connectivity mask
2: Exposed intertidal (low confidence)
Pixels with water index/tide correlations higher than 0.15.
Pixels must be located within the coastal connectivity mask
and not be intertidal elevation (high confidence) pixels
3: Exposed intertidal (high confidence)
Pixels that are included in the intertidal elevation dataset
4: Inland waters
Pixels that are wet in more than 50% of observations and fall
outside of the coastal connectivity mask.
5: Land
Pixels that are wet in less than 50% of observations
"""

# Identify dataset geobox
geobox = ds.odc.geobox

# Load other inputs
ocean_mask = load_ocean_mask(dc, geobox)
urban_mask = load_aclum_mask(dc, geobox)
bathy_mask = load_topobathy_mask(dc, geobox)
costdist_mask, _ = load_connectivity_mask(dc, geobox, buffer=buffer)
geobox = dem.odc.geobox

# Identify any pixels that are nodata in both frequency and bathy mask
# (this is a temporary hack due to us not having any other way of telling
# ocean nodata pixels apart from inland nodata pixels
is_nan = (ds.qa_ndwi_freq.isnull()) & bathy_mask
# Identify any pixels that are nodata in frequency
is_nan = freq.isnull()

# Spilt pixels into those that were mostly wet vs mostly dry.
# Identify subset of mostly wet pixels that were inland
mostly_dry = (ds.qa_ndwi_freq < 50) & ~is_nan
mostly_wet = (ds.qa_ndwi_freq >= 50) & ~is_nan
mostly_wet_inland = mostly_wet & ~costdist_mask
mostly_dry = (freq < 0.50) & ~is_nan
mostly_wet = (freq >= 0.50) & ~is_nan
mostly_wet_inland = mostly_wet & ~coastal_mask

# Identify low-confidence pixels as those with greater than 0.15
# correlation. Use connectivity mask to mask out any that are "inland"
intertidal_lc = (ds.qa_ndwi_corr >= min_correlation) & costdist_mask
intertidal_lc = (corr >= min_correlation) & coastal_mask

# Identify high confidence intertidal as those in our elevation data
intertidal_hc = ds.elevation.notnull()
intertidal_hc = dem.notnull()

# Identify likely misclassified urban pixels as those that overlap with
# the urban mask
urban_misclass = mostly_wet_inland & urban_mask

# Combine all classifications - this is done one-by-one, pasting each
# new layer over the top of the existing data
extents = xr_zeros(geobox=geobox, dtype="int16") # start with 0
extents = xr_zeros(geobox=geobox, dtype="int16") + 255 # start with 255
extents.values[mostly_wet] = 1 # Add in mostly wet pixels
extents.values[mostly_wet_inland] = 4 # Add in mostly wet inland pixels on top
extents.values[urban_misclass] = (
Expand All @@ -364,6 +367,6 @@ class that are sieved out to remove small noisy dataset
extents.values[intertidal_hc] = 3

# Export to file
extents.attrs["nodata"] = 0
extents.attrs["nodata"] = 255

return extents
Loading