@@ -33,10 +33,39 @@ def rename_bands(ds, old_string, new_string):
3333 return ds_renamed
3434
3535
36+ def tidal_thresholds (
37+ tides_highres ,
38+ threshold_lowtide = 0.15 ,
39+ threshold_hightide = 0.85 ,
40+ min_obs = 0 ,
41+ ):
42+ # Calculate per-pixel integer rankings for each tide height
43+ rank_n = tides_highres .rank (dim = "time" )
44+
45+ # Calculate low and high ranking thresholds from total rankings.
46+ # Low threshold needs to be rounded up ("ceil"), and high tide
47+ # rounded down ("floor") to ensure we capture all matching values.
48+ rank_max = rank_n .max (dim = "time" )
49+ rank_thresh_low = np .ceil (rank_max * threshold_lowtide )
50+ rank_thresh_high = np .floor (rank_max * threshold_hightide )
51+
52+ # Update thresholds to ensure minimum number of valid observations
53+ if min_obs > 0 :
54+ rank_thresh_low = np .maximum (rank_thresh_low , min_obs )
55+ rank_thresh_high = np .minimum (rank_thresh_high , rank_max - min_obs )
56+
57+ # Calculate tide thresholds by masking tides by ranking threshold
58+ tide_thresh_low = tides_highres .where (rank_n <= rank_thresh_low ).max (dim = "time" )
59+ tide_thresh_high = tides_highres .where (rank_n >= rank_thresh_high ).min (dim = "time" )
60+
61+ return tide_thresh_low , tide_thresh_high
62+
63+
3664def tidal_composites (
3765 satellite_ds ,
38- threshold_lowtide = 0.2 ,
39- threshold_hightide = 0.8 ,
66+ threshold_lowtide = 0.15 ,
67+ threshold_hightide = 0.85 ,
68+ min_obs = 0 ,
4069 eps = 1e-4 ,
4170 cpus = None ,
4271 max_iters = 10000 ,
@@ -55,7 +84,7 @@ def tidal_composites(
5584 to filter satellite data to low and high tide images prior to
5685 loading it into memory, allowing more efficient processing.
5786
58- Based on the method described in:
87+ Pixel-based implementation of the method originally published in:
5988
6089 Sagar, S., Phillips, C., Bala, B., Roberts, D., & Lymburner, L.
6190 (2018). Generating Continental Scale Pixel-Based Surface Reflectance
@@ -67,9 +96,12 @@ def tidal_composites(
6796 satellite_ds : xarray.Dataset
6897 A satellite data time series containing spectral bands.
6998 threshold_lowtide : float, optional
70- Quantile used to identify low tide observations, by default 0.2 .
99+ Quantile used to identify low tide observations, by default 0.15 .
71100 threshold_hightide : float, optional
72- Quantile used to identify high tide observations, by default 0.8.
101+ Quantile used to identify high tide observations, by default 0.85.
102+ min_obs : int, optional
103+ Minimum number of clear observations to enforce when calculating tide
104+ height thresholds. Defaults to 0, which will not apply any minimum.
73105 eps: float, optional
74106 Termination criteria passed on to the geomedian algorithm.
75107 cpus: int, optional
@@ -134,14 +166,22 @@ def tidal_composites(
134166 tides_highres = tides_highres .where (nodata_array )
135167
136168 # Calculate low and high tide thresholds from masked tide data
137- log .info (f"{ run_id } )
138- threshold_ds = xr_quantile (
139- src = tides_highres .to_dataset (),
140- quantiles = [threshold_lowtide , threshold_hightide ],
141- nodata = np .nan ,
169+ log .info (
170+ f"{ run_id } { min_obs }
171+ )
172+ # threshold_ds = xr_quantile(
173+ # src=tides_highres.to_dataset(),
174+ # quantiles=[threshold_lowtide, threshold_hightide],
175+ # nodata=np.nan,
176+ # )
177+ # low_threshold = threshold_ds.isel(quantile=0).tide_height.drop("quantile")
178+ # high_threshold = threshold_ds.isel(quantile=-1).tide_height.drop("quantile")
179+ low_threshold , high_threshold = tidal_thresholds (
180+ tides_highres = tides_highres ,
181+ threshold_lowtide = threshold_lowtide ,
182+ threshold_hightide = threshold_hightide ,
183+ min_obs = min_obs ,
142184 )
143- low_threshold = threshold_ds .isel (quantile = 0 ).tide_height .drop ("quantile" )
144- high_threshold = threshold_ds .isel (quantile = - 1 ).tide_height .drop ("quantile" )
145185
146186 # Create masks for selecting satellite observations below and above the
147187 # low and high tide thresholds
@@ -270,14 +310,21 @@ def tidal_composites(
270310@click .option (
271311 "--threshold_lowtide" ,
272312 type = float ,
273- default = 0.2 ,
274- help = "The quantile used to identify low tide observations. " " Defaults to 0.2 ."
313+ default = 0.15 ,
314+ help = "The quantile used to identify low tide observations. Defaults to 0.15 ." ,
275315)
276316@click .option (
277317 "--threshold_hightide" ,
278318 type = float ,
279- default = 0.8 ,
280- help = "The quantile used to identify high tide observations. " "Defaults to 0.8." ,
319+ default = 0.85 ,
320+ help = "The quantile used to identify high tide observations. Defaults to 0.85." ,
321+ )
322+ @click .option (
323+ "--min_obs" ,
324+ type = int ,
325+ default = 0 ,
326+ help = "Minimum number of clear observations to enforce when calculating tide "
327+ "height thresholds. Defaults to 0, which will not apply any minimum." ,
281328)
282329@click .option (
283330 "--mask_sunglint" ,
@@ -355,6 +402,7 @@ def tidal_composites_cli(
355402 resolution ,
356403 threshold_lowtide ,
357404 threshold_hightide ,
405+ min_obs ,
358406 mask_sunglint ,
359407 include_coastal_aerosol ,
360408 eps ,
@@ -429,7 +477,7 @@ def tidal_composites_cli(
429477 )
430478
431479 # Fail early if not enough observations
432- if len (satellite_ds .time ) < 20 :
480+ if len (satellite_ds .time ) < 50 :
433481 raise Exception (
434482 "Insufficient satellite data available to process composites; skipping."
435483 )
@@ -440,6 +488,7 @@ def tidal_composites_cli(
440488 satellite_ds = satellite_ds ,
441489 threshold_lowtide = threshold_lowtide ,
442490 threshold_hightide = threshold_hightide ,
491+ min_obs = min_obs ,
443492 eps = eps ,
444493 cpus = cpus ,
445494 max_iters = max_iters ,
0 commit comments