Merge pull request #84 from GeoscienceAustralia/upgrade/dem-logic

DEM Logic and Backscatter Convention Specification

Author: abradley60

.dockerignore

@@ -0,0 +1 @@
+.pixi

Docker/Dockerfile

@@ -18,16 +18,15 @@ chown -R rtc_user:rtc_user /home/rtc_user/OPERA
 
 USER rtc_user
 
-ENV CONDA_PREFIX=/home/rtc_user/miniconda3
-
-# Install Miniconda
+# Install Miniforge
 WORKDIR /home/rtc_user
-RUN curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -o miniconda.sh &&\
-bash miniconda.sh -b -p ${CONDA_PREFIX} &&\
-rm $HOME/miniconda.sh
+ENV CONDA_PREFIX=/home/rtc_user/miniforge3
+ENV PATH=$CONDA_PREFIX/bin:$PATH
 
-ENV PATH=${CONDA_PREFIX}/bin:${PATH}
-RUN ${CONDA_PREFIX}/bin/conda init bash
+RUN curl -sSL https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh -o miniforge.sh &&\
+    bash miniforge.sh -b -p $CONDA_PREFIX &&\
+    rm miniforge.sh &&\
+    $CONDA_PREFIX/bin/conda init bash
 
 # Download and install the Geoscience Australia fork of the opera/RTC project
 # UPDATE VERSION AND TAG HERE IF NEEDED

docs/workflows/aws.md

@@ -61,6 +61,7 @@ The AWS pipeline runs using a docker container. At runtime, the script [run_aws_
 --output_crs=""
 --dem_type="cop_glo30"
 --product="RTC_S1"
+--backscatter_convention="gamma0" 
 --s3_bucket="deant-data-public-dev"
 --s3_project_folder="experimental"
 --collection="s1_rtc_c1"
@@ -82,9 +83,10 @@ The AWS pipeline runs using a docker container. At runtime, the script [run_aws_
 - `output_crs` -> The target crs of the products. If not specified, the UTM of the scene center will be used. Expects integer values (e.g. 3031)
 - `dem_type` -> The type of digital elevation model (DEM) to download and use for processing. Supported values: `cop_glo30`, `REMA_32`, `REMA_10`, `REMA_2`.
 - `product` -> The product being created with the workflow. Must be `RTC_S1` or `RTC_S1_STATIC`.
+- `backscatter_convention` -> the output backscatter convention from the workflow. Note sigma0 data is referenced to the DEM. To create sigma0 ellipsoid referenced data, the beta0 layer and static incidence_angle layer is required; sigma0_ellipsoid = beta0*sin(incidence_angle).
 - `s3_bucket` -> the bucket to upload the products
 - `s3_project_folder` -> The project folder to upload to.
-- `collection` -> The collection which the set of products belongs. Must end with 'cX' where X is a valid integer reffering to the collection number of the product. e.g. rtc_s1_c1.
+- `collection` -> The collection which the set of products belongs. Must end with 'cX' where X is a valid integer referring to the collection number of the product. e.g. rtc_s1_c1.
 - `make_existing_products` -> Whether to generate products even if they already exist in AWS S3 under the specified product folder path `s3_bucket/s3_project_folder/collection/...`. 
   - **WARNING** - Passing this flag will create duplicate files and overwrite existing metadata, which may affect downstream workflows.
 - `skip_upload_to_s3` -> Make the products, but skip uploading them to S3.

environment.yaml

@@ -32,6 +32,6 @@ dependencies:
   - ruamel-yaml>=0.18.10
   - s1-orbits>=0.1.3
   - cdsetool>=0.2.13
-  - dem-handler @ git+https://github.com/GeoscienceAustralia/dem-handler.git@v0.2.2
+  - dem-handler @ git+https://github.com/GeoscienceAustralia/dem-handler.git@v0.2.3
   - python-dotenv>=1.1.1, <2
   - pystac[validation]==1.13.0

pixi.lock

@@ -24,7 +24,7 @@ dependencies = [
     "ruamel-yaml>=0.18.10", 
     "s1-orbits>=0.1.3", 
     "cdsetool>=0.2.13",
-    "dem-handler @ git+https://github.com/GeoscienceAustralia/dem-handler.git@v0.2.2", 
+    "dem-handler @ git+https://github.com/GeoscienceAustralia/dem-handler.git@v0.2.3", 
     "python-dotenv>=1.1.1,<2", 
     "pystac[validation]==1.13.0", 
 ]

pyproject.toml

@@ -20,7 +20,6 @@
 )
 
 from sar_pipeline.aws.preparation.config import RTCConfigManager
-from sar_pipeline.aws.metadata.h5 import update_h5_file_metadata
 from sar_pipeline.aws.metadata.stac import BurstH5toStacManager
 from sar_pipeline.aws.metadata.odc import (
     make_static_layer_base_url,
@@ -35,11 +34,14 @@
 from dem_handler.utils.spatial import (
     check_s1_bounds_cross_antimeridian,
     get_correct_bounds_from_shape_at_antimeridian,
+    check_dem_type_in_bounds,
 )
 
 logging.basicConfig(level=logging.INFO)
 logger = logging.getLogger(__name__)
 
+VALID_DEMS = ["cop_glo30", "REMA_32", "REMA_10", "REMA_2"]
+
 
 @click.command()
 @click.option(
@@ -69,15 +71,27 @@
 @click.option(
     "--dem-type",
     required=True,
-    type=click.Choice(["cop_glo30", "REMA_32", "REMA_10", "REMA_2"]),
-    help="The type of DEM that should be downloaded for processing the scene.",
+    type=str,
+    help="The type of DEM that should be downloaded for processing the scene. "
+    "Can be passed as a string or list of preferences separated by a space. "
+    "If DEM data does not exist in the area of the first preference, the next will be used. "
+    "E.g. `--dem-type REMA_32 cop_glo30` will first look for the Antarctic specific REMA DEM @32m before settling on the cop_glo30. "
+    f"Values must be one of {VALID_DEMS}",
 )
 @click.option(
     "--product",
     required=True,
     type=click.Choice(["RTC_S1", "RTC_S1_STATIC"]),
     help="The product to be made",
 )
+@click.option(
+    "--backscatter-convention",
+    required=False,
+    default="gamma0",
+    type=click.Choice(["gamma0", "sigma0", "beta0"]),
+    help="Backscatter convention of the product to be made (gamma0, sigma0 or beta0). "
+    "Note, sigma0 data is referenced to the DEM. See docs for information on creating ellipsoid referenced sigma0 data.",
+)
 @click.option(
     "--s3-bucket",
     required=True,
@@ -187,6 +201,7 @@ def get_data_for_scene_and_make_run_config(
     output_crs,
     dem_type,
     product,
+    backscatter_convention,
     s3_bucket,
     s3_project_folder,
     collection,
@@ -212,6 +227,10 @@ def get_data_for_scene_and_make_run_config(
     logger.info(f"Data source for scene download : {scene_data_source}")
     logger.info(f"Data source for orbit download : {orbit_data_source}")
 
+    # get the preference list of dem types
+    dem_types = dem_type.split(" ")
+    logger.info(f"The order of DEM preference is : {dem_types}")
+
     # make the base .yaml for RTC processing
     if product == "RTC_S1":
         RTC_RUN_CONFIG = RTCConfigManager(base_config="S1_RTC.yaml")
@@ -220,21 +239,22 @@ def get_data_for_scene_and_make_run_config(
     else:
         raise ValueError("product must be S1_RTC or S1_RTC_STATIC")
 
+    if backscatter_convention not in ["gamma0", "sigma0", "beta0"]:
+        raise ValueError("backscatter_convention must be one of gamma0, sigma0, beta0")
+
     # ensure the collection ends with cX, where X is a positive integer.
     # Raise error for invalid naming
     _ = get_collection_number(collection)
 
     # sub-folders for downloads
     orbit_folder = download_folder / "orbits"
     scene_folder = download_folder / "scenes"
-    dem_folder = download_folder / "dem" / dem_type
 
     if make_folders:
         logger.info(f"Making output folders if not existing")
         download_folder.mkdir(parents=True, exist_ok=True)
         orbit_folder.mkdir(parents=True, exist_ok=True)
         scene_folder.mkdir(parents=True, exist_ok=True)
-        dem_folder.mkdir(parents=True, exist_ok=True)
         out_folder.mkdir(parents=True, exist_ok=True)
         scratch_folder.mkdir(parents=True, exist_ok=True)
         run_config_save_path.parent.mkdir(parents=True, exist_ok=True)
@@ -331,9 +351,6 @@ def get_data_for_scene_and_make_run_config(
     )
     logger.info(f"File downloaded to : {ORBIT_PATHS[0]}")
 
-    # download the dem
-    DEM_PATH = dem_folder / f"{scene}_dem.tif"
-
     # get the shape of the area covering the bursts to be processed
     burst_geoms_to_process = [
         all_scene_burst_info[id_]["geometry"] for id_ in burst_id_list_to_process
@@ -363,13 +380,36 @@ def get_data_for_scene_and_make_run_config(
         bounds = get_correct_bounds_from_shape_at_antimeridian(scene_bounds)
         logger.info(f"The corrected scene bounds are : {bounds}")
         # increase the buffer to ensure DEM sufficiently covers area
-        # shape is a little odd due to merging either side of AM
+        # cop30 shape is a little odd due to merging either side of AM
         cop30_buffer_degrees = 0.8
     else:
         logger.info(f"Downloading DEM data for the combined burst bounds")
         bounds = combined_burst_bounds
         cop30_buffer_degrees = 0.3
 
+    logger.info(f"Finding the best DEM in order of preference: {dem_types}")
+
+    for dem_type in dem_types:
+        if dem_type == "cop_glo30":
+            dem_resolution = 30
+        elif dem_type in ["REMA_32", "REMA_10", "REMA_2"]:
+            dem_resolution = int(dem_type.split("_")[1])
+        dem_in_bounds = check_dem_type_in_bounds(dem_type, dem_resolution, bounds)
+        if dem_in_bounds:
+            # dem has data in the bounds we want, exit with dem_type set
+            break
+        elif dem_type == dem_types[-1] and not dem_in_bounds:
+            logger.warning(
+                "No DEM data could not be found in the requested bounds. "
+                "The bursts will be processed with ellipsoidal heights. "
+                "Change `dem_type` if this is not desired."
+            )
+
+    dem_folder = download_folder / "dem" / dem_type
+    DEM_PATH = dem_folder / f"{scene}_dem.tif"
+    if make_folders:
+        dem_folder.mkdir(parents=True, exist_ok=True)
+
     logger.info(f"Downloading DEM type `{dem_type}` to path : {DEM_PATH}")
     if dem_type == "cop_glo30":
         get_cop30_dem_for_bounds(
@@ -385,7 +425,6 @@ def get_data_for_scene_and_make_run_config(
             download_geoid=True,
         )
     elif dem_type in ["REMA_32", "REMA_10", "REMA_2"]:
-        dem_resolution = int(dem_type.split("_")[1])
         get_rema_dem_for_bounds(
             bounds=bounds,
             bounds_src_crs=4326,
@@ -413,6 +452,14 @@ def get_data_for_scene_and_make_run_config(
     RTC_RUN_CONFIG.set(f"{gk}.input_file_group.orbit_file_path", [str(ORBIT_PATHS[0])])
     RTC_RUN_CONFIG.set(f"{gk}.dynamic_ancillary_file_group.dem_file", str(DEM_PATH))
 
+    # set the backscatter convention and if rtc_apply is true
+    if backscatter_convention == "beta0":
+        RTC_RUN_CONFIG.set(f"{gk}.processing.apply_rtc", False)
+        # note output_type is still set to gamma0, but it is ignored as apply_rtc=False
+    else:
+        RTC_RUN_CONFIG.set(f"{gk}.processing.apply_rtc", True)
+        RTC_RUN_CONFIG.set(f"{gk}.processing.rtc.output_type", backscatter_convention)
+
     # set the dem input source
     if dem_type == "cop_glo30":
         demSource = "https://registry.opendata.aws/copernicus-dem/"
@@ -495,6 +542,13 @@ def get_data_for_scene_and_make_run_config(
     type=click.Choice(["RTC_S1", "RTC_S1_STATIC"]),
     help="The product being made. Determines bucket structure",
 )
+@click.option(
+    "--backscatter-convention",
+    required=False,
+    default="gamma0",
+    type=click.Choice(["gamma0", "sigma0", "beta0"]),
+    help="Backscatter convention of the product to be made (gamma0, sigma0 or beta0)",
+)
 @click.option(
     "--collection",
     required=True,
@@ -550,6 +604,7 @@ def make_rtc_opera_stac_and_upload_bursts(
     results_folder,
     run_config_path,
     product,
+    backscatter_convention,
     collection,
     s3_bucket,
     s3_project_folder,
@@ -579,17 +634,13 @@ def make_rtc_opera_stac_and_upload_bursts(
                 f"This error might be caused by repeat runs. Delete duplicate files or change run settings."
             )
         burst_h5_filepath = burst_folder / burst_h5_files[0]
-        burst_tif_files = list(burst_folder.glob("*.tif"))
-        # update GeoTIFFs and h5 metadata to reflect GA processing
-        logging.info(f"Updating .h5 file with GA parameters")
-        update_h5_file_metadata(burst_h5_filepath)
-        # TODO Update .tif files with GA parameters
         # make the stac metadata from the .h5 metadata
         logging.info(f"Making stac metadata from .h5 file")
         # initialise the class to convert data from the .h5 to a stac doc
         burst_stac_manager = BurstH5toStacManager(
             h5_filepath=burst_h5_filepath,
             product=product,
+            backscatter_convention=backscatter_convention,
             collection=collection,
             s3_bucket=s3_bucket,
             s3_project_folder=s3_project_folder,