11import os
22import sys
3+
34import click
4- import numpy as np
5- import xarray as xr
65import datacube
6+ import numpy as np
77import odc .geo .xr
8- from odc .geo .geom import BoundingBox
8+ import xarray as xr
9+ from datacube .utils .aws import configure_s3_access
10+ from dea_tools .dask import create_local_dask_cluster
11+ from eo_tides .eo import pixel_tides
912from odc .algo import (
1013 int_geomedian ,
1114 keep_good_only ,
12- xr_quantile ,
1315)
14- from datacube .utils .aws import configure_s3_access
15- from eo_tides .eo import pixel_tides
16- from dea_tools .dask import create_local_dask_cluster
16+ from odc .geo .geom import BoundingBox
1717
18- from intertidal .utils import configure_logging
1918from intertidal .io import (
19+ export_dataset_metadata ,
2020 load_data ,
2121 prepare_for_export ,
2222 tidal_metadata ,
23- export_dataset_metadata ,
2423)
24+ from intertidal .utils import configure_logging
2525
2626
2727# Function to rename the bands
2828def rename_bands (ds , old_string , new_string ):
2929 # Create a new dataset with renamed bands
30- ds_renamed = ds .rename (
31- {band : band .replace (old_string , new_string ) for band in ds .data_vars }
32- )
30+ ds_renamed = ds .rename ({band : band .replace (old_string , new_string ) for band in ds .data_vars })
3331 return ds_renamed
3432
3533
@@ -66,16 +64,38 @@ def tidal_thresholds(
6664
6765
6866def filter_granules (dataset ):
69- """
70- Return False for any Sentinel-2 dataset with a MGRS
67+ """Return False for any Sentinel-2 dataset with a MGRS
7168 granule region code in the list of bad region codes.
7269 """
73- drop_list = ["50HKG" , "50HNF" , "51LWD" , "51LXE" , "51LZF" ,
74- "52LBL" , "52LCL" , "52LDK" , "53HNA" , "53LRC" ,
75- "54GYU" , "54LWR" , "54LXR" , "54LYR" , "55GBP" ,
76- "55KEA" , "55KFV" , "55KGV" , "55KHT" , "55KHU" ,
77- "56KKC" , "56KLC" , "56KMC" , "56KMV" , "56KNU" ,
78- "54LWQ" , "54LWP" ]
70+ drop_list = [
71+ "50HKG" ,
72+ "50HNF" ,
73+ "51LWD" ,
74+ "51LXE" ,
75+ "51LZF" ,
76+ "52LBL" ,
77+ "52LCL" ,
78+ "52LDK" ,
79+ "53HNA" ,
80+ "53LRC" ,
81+ "54GYU" ,
82+ "54LWR" ,
83+ "54LXR" ,
84+ "54LYR" ,
85+ "55GBP" ,
86+ "55KEA" ,
87+ "55KFV" ,
88+ "55KGV" ,
89+ "55KHT" ,
90+ "55KHU" ,
91+ "56KKC" ,
92+ "56KLC" ,
93+ "56KMC" ,
94+ "56KMV" ,
95+ "56KNU" ,
96+ "54LWQ" ,
97+ "54LWP" ,
98+ ]
7999 return dataset .metadata .region_code not in drop_list
80100
81101
@@ -92,8 +112,7 @@ def tidal_composites(
92112 run_id = None ,
93113 log = None ,
94114):
95- """
96- Calculates Geometric Median composites of the coastal zone at low
115+ """Calculates Geometric Median composites of the coastal zone at low
97116 and high tide using satellite imagery and tidal modeling.
98117
99118 This function uses tools from `odc.algo` to keep data in its
@@ -157,8 +176,8 @@ def tidal_composites(
157176 ds_hightide : xarray.Dataset
158177 xarray.Dataset object containing a geomedian of the observations
159178 with the highest X quantile tide values for each pixel.
160- """
161179
180+ """
162181 # Set up logs if no log is passed in
163182 if log is None :
164183 log = configure_logging ()
@@ -184,9 +203,7 @@ def tidal_composites(
184203 tides_highres = tides_highres .where (nodata_array )
185204
186205 # Calculate low and high tide thresholds from masked tide data
187- log .info (
188- f"{ run_id } : Calculating low and high tide thresholds with minimum { min_obs } observations"
189- )
206+ log .info (f"{ run_id } : Calculating low and high tide thresholds with minimum { min_obs } observations" )
190207 low_threshold , high_threshold = tidal_thresholds (
191208 tides_highres = tides_highres ,
192209 threshold_lowtide = threshold_lowtide ,
@@ -206,13 +223,9 @@ def tidal_composites(
206223 ds_high = satellite_ds .sel (time = high_keep )
207224
208225 # Load low and high subsets of data into memory
209- log .info (
210- f"{ run_id } : Loading { len (ds_low .time )} low tide satellite images into memory"
211- )
226+ log .info (f"{ run_id } : Loading { len (ds_low .time )} low tide satellite images into memory" )
212227 ds_low .load ()
213- log .info (
214- f"{ run_id } : Loading { len (ds_high .time )} high tide satellite images into memory"
215- )
228+ log .info (f"{ run_id } : Loading { len (ds_high .time )} high tide satellite images into memory" )
216229 ds_high .load ()
217230
218231 # Use `keep_good_only` to set any pixels outside of the tide masks to nodata
@@ -239,9 +252,7 @@ def tidal_composites(
239252 # Calculate clear count (both low and high tide clear counts
240253 # are identical, so we can just use one)
241254 log .info (f"{ run_id } : Calculating clear counts" )
242- ds_lowtide ["qa_count_clear" ] = (
243- (ds_low_masked .nbart_red != nodata ).sum (dim = "time" ).astype ("int16" )
244- )
255+ ds_lowtide ["qa_count_clear" ] = (ds_low_masked .nbart_red != nodata ).sum (dim = "time" ).astype ("int16" )
245256
246257 # Add low and high tide thresholds to the output datasets
247258 ds_lowtide ["qa_low_threshold" ] = low_threshold
@@ -255,8 +266,7 @@ def tidal_composites(
255266 "--study_area" ,
256267 type = str ,
257268 required = True ,
258- help = "A string providing a GridSpec tile ID (e.g. in the form "
259- "'x123y123') to run the analysis on." ,
269+ help = "A string providing a GridSpec tile ID (e.g. in the form 'x123y123') to run the analysis on." ,
260270)
261271@click .option (
262272 "--start_date" ,
@@ -287,7 +297,7 @@ def tidal_composites(
287297 "--output_version" ,
288298 type = str ,
289299 required = True ,
290- help = "The version number to use for output files and metadata (e.g. " " '0.0.1')." ,
300+ help = "The version number to use for output files and metadata (e.g. '0.0.1')." ,
291301)
292302@click .option (
293303 "--output_dir" ,
@@ -300,15 +310,13 @@ def tidal_composites(
300310 "--product_maturity" ,
301311 type = str ,
302312 default = "provisional" ,
303- help = "Product maturity metadata to use for the output dataset. "
304- "Defaults to 'provisional', can also be 'stable'." ,
313+ help = "Product maturity metadata to use for the output dataset. Defaults to 'provisional', can also be 'stable'." ,
305314)
306315@click .option (
307316 "--dataset_maturity" ,
308317 type = str ,
309318 default = "final" ,
310- help = "Dataset maturity metadata to use for the output dataset. "
311- "Defaults to 'final', can also be 'interim'." ,
319+ help = "Dataset maturity metadata to use for the output dataset. Defaults to 'final', can also be 'interim'." ,
312320)
313321@click .option (
314322 "--resolution" ,
@@ -350,8 +358,7 @@ def tidal_composites(
350358 "--gqa_filter/--no-gqa_filter" ,
351359 type = bool ,
352360 default = True ,
353- help = "Whether to filter scenes when loading data based on gqa values. "
354- "Defaults to True" ,
361+ help = "Whether to filter scenes when loading data based on gqa values. Defaults to True" ,
355362)
356363@click .option (
357364 "--include_coastal_aerosol/--no-include_coastal_aerosol" ,
@@ -433,14 +440,13 @@ def tidal_composites_cli(
433440 overwrite ,
434441):
435442 # Create sample filename to test if data exists on file system
436- filename = f"{ output_dir } ga_s2_tidal_composites_cyear_3/{ output_version .replace ('.' ,'-' )} /{ study_area [:4 ]} /{ study_area [4 :]} /{ label_date } --P1Y/ga_s2_tidal_composites_cyear_3_{ study_area } _{ label_date } --P1Y_final.stac-item.json"
443+ filename = f"{ output_dir } ga_s2_tidal_composites_cyear_3/{ output_version .replace ('.' , '-' )} /{ study_area [:4 ]} /{ study_area [4 :]} /{ label_date } --P1Y/ga_s2_tidal_composites_cyear_3_{ study_area } _{ label_date } --P1Y_final.stac-item.json"
437444
438445 process_tile = True
439446 if overwrite :
440447 process_tile = True
441- else :
442- if os .path .exists (filename ):
443- process_tile = False
448+ elif os .path .exists (filename ):
449+ process_tile = False
444450
445451 # Create a unique run ID based on input params and use for logs
446452 input_params = locals ()
@@ -454,9 +460,7 @@ def tidal_composites_cli(
454460 configure_s3_access (cloud_defaults = True , aws_unsigned = aws_unsigned )
455461
456462 if process_tile :
457-
458463 try :
459-
460464 # Create local dask cluster to improve data load time
461465 client = create_local_dask_cluster (return_client = True )
462466
@@ -466,9 +470,7 @@ def tidal_composites_cli(
466470 # Use a custom polygon if in testing mode
467471 if study_area == "testing" :
468472 log .info (f"{ run_id } : Running in testing mode using custom study area" )
469- geom = BoundingBox (
470- 467510 , - 1665790 , 468260 , - 1664840 , crs = "EPSG:3577"
471- ).polygon
473+ geom = BoundingBox (467510 , - 1665790 , 468260 , - 1664840 , crs = "EPSG:3577" ).polygon
472474 else :
473475 geom = None
474476
@@ -494,15 +496,11 @@ def tidal_composites_cli(
494496 dtype = "int16" ,
495497 dataset_predicate = filter_granules ,
496498 )
497- log .info (
498- f"{ run_id } : Found { len (satellite_ds .time )} satellite data timesteps"
499- )
499+ log .info (f"{ run_id } : Found { len (satellite_ds .time )} satellite data timesteps" )
500500
501501 # Fail early if not enough observations
502502 if len (satellite_ds .time ) < 50 :
503- raise Exception (
504- "Insufficient satellite data available to process composites; skipping."
505- )
503+ raise Exception ("Insufficient satellite data available to process composites; skipping." )
506504
507505 # Calculate high and low tide geomedian composites
508506 log .info (f"{ run_id } : Running DEA Tidal Composites workflow" )
@@ -528,9 +526,7 @@ def tidal_composites_cli(
528526 ds_tidalcomposites = xr .merge ([ds_lowtide , ds_hightide ])
529527
530528 # Ensure spatial information is still attached
531- ds_tidalcomposites = odc .geo .xr .assign_crs (
532- ds_tidalcomposites , satellite_ds .odc .crs
533- )
529+ ds_tidalcomposites = odc .geo .xr .assign_crs (ds_tidalcomposites , satellite_ds .odc .crs )
534530
535531 custom_dtypes = {
536532 "low_coastal_aerosol" : (np .int16 , - 999 ),
0 commit comments