Add support for fetching source HLS scenes to OPERA metadata export

next_pass.py

@@ -134,6 +134,11 @@ def create_parser() -> argparse.ArgumentParser:
             "in format YYYY-MM-DDTHH:MM"
         ),
     )
+    parser.add_argument(
+        "--include-hls",
+        action="store_true",
+        help="Fetch and append corresponding source HLS scenes for OPERA HLS products.",
+    )
     return parser
 
 
@@ -254,6 +259,7 @@ def run_next_pass(
     functionality: str = "both",
     compute_cloudiness: bool = False,
     compute_tide: bool = False,
+    include_hls: bool = False,
     products: List[str] | str | None = None,
     satellites: List[str] | str | None = None,
 ):
@@ -286,6 +292,8 @@ def run_next_pass(
         cli_args.append("-c")
     if compute_tide:
         cli_args.append("-t")
+    if include_hls:
+        cli_args.append("--include-hls")
 
     if date:
         cli_args += ["-d", date]
@@ -389,7 +397,10 @@ def main(cli_args: Any = None):
             timestamp_dir,
         )
         export_opera_products(
-            results_opera, timestamp_dir, compute_cloudiness=args.cloudiness
+            results_opera,
+            timestamp_dir,
+            compute_cloudiness=args.cloudiness,
+            include_hls=getattr(args, "include_hls", False),
         )
         make_opera_granule_map(results_opera, args.bbox, timestamp_dir)
 

utils/opera_products.py

@@ -3,6 +3,7 @@
 from datetime import datetime, timezone
 from pathlib import Path
 
+import earthaccess
 import leafmap
 import pandas as pd
 from dateutil.relativedelta import relativedelta
@@ -172,6 +173,19 @@ def find_print_available_opera_products(
     return results_dict
 
 
+def fetch_hls_granule_links(granule_id: str) -> list:
+    """Fetch the CMR metadata for a specific HLS granule ID to get download links."""
+    collection = "HLSS30" if "S30" in granule_id else "HLSL30"
+    try:
+        results = earthaccess.search_data(short_name=collection, granule_ur=granule_id)
+        if results:
+            return results[0].data_links()
+    except Exception as e:
+        LOGGER.error("Failed to fetch HLS granule %s: %s", granule_id, e)
+
+    return []
+
+
 def describe_cloud_cover(cover_percent: float) -> str:
     """Return a short description string for a given cloud cover %."""
     if cover_percent > 75:
@@ -188,7 +202,11 @@ def describe_cloud_cover(cover_percent: float) -> str:
 
 
 def export_opera_products(
-    results_dict: dict, timestamp_dir, result_s1=None, compute_cloudiness: bool = True
+    results_dict: dict,
+    timestamp_dir,
+    result_s1=None,
+    compute_cloudiness: bool = True,
+    include_hls: bool = False,
 ) -> None:
     """
     Export OPERA products to an Excel file and log cloudiness summary.
@@ -203,6 +221,8 @@ def export_opera_products(
         Currently unused, kept for API compatibility.
     compute_cloudiness : bool
         Whether to compute cloudiness from CLOUD layers. Set to False to skip and save time.
+    include_hls : bool
+        Whether to include HLS products in the export. Set to False to skip HLS products in xls output.
     """
     output_file = timestamp_dir / "opera_products_metadata.xlsx"
     wb = Workbook()
@@ -229,6 +249,18 @@ def export_opera_products(
         "Download URL CSLC-VV",
         "Geometry (WKT)",
     ]
+
+    if include_hls:
+        headers.extend(
+            [
+                "Source HLS Granule ID",
+                "HLS Download URL (B04/Red)",
+                "HLS Download URL (B03/Green)",
+                "HLS Download URL (B02/Blue)",
+                "HLS Download URL (B8A/B05/NIR)",
+                "HLS Download URL (Fmask)",
+            ]
+        )
     ws.append(headers)
 
     # Apply bold to header cells
@@ -323,27 +355,114 @@ def export_opera_products(
                     overall_cloudy_area += area * cloud_cover_percent / 100.0
                     overall_area += area
 
-            # Write data row
-            ws.append(
-                [
-                    dataset,
-                    granule_id,
-                    start_time,
-                    end_time,
-                    cloud_cover_percent,
-                    urls["water"],
-                    urls["bwater"],
-                    urls["water_conf"],
-                    urls["veg_anom_max"],
-                    urls["veg_dist_status"],
-                    urls["veg_dist_date"],
-                    urls["veg_dist_conf"],
-                    urls["rtc-vv"],
-                    urls["rtc-vh"],
-                    urls["cslc-vv"],
-                    geom_wkt,
-                ]
-            )
+            # Write base data row
+            row_data = [
+                dataset,
+                granule_id,
+                start_time,
+                end_time,
+                cloud_cover_percent,
+                urls["water"],
+                urls["bwater"],
+                urls["water_conf"],
+                urls["veg_anom_max"],
+                urls["veg_dist_status"],
+                urls["veg_dist_date"],
+                urls["veg_dist_conf"],
+                urls["rtc-vv"],
+                urls["rtc-vh"],
+                urls["cslc-vv"],
+                geom_wkt,
+            ]
+
+            # Only check for HLS granules if it is an HLS-derived OPERA product
+            if include_hls:
+                # Initialize variables here to guarantee they exist
+                hls_granule_id = "N/A"
+                hls_red = "N/A"
+                hls_green = "N/A"
+                hls_blue = "N/A"
+                hls_nir = "N/A"
+                hls_fmask = "N/A"
+
+                if "HLS" in dataset:
+                    hls_links = []
+
+                    # Try to extract directly from InputGranules (Standard for DSWx)
+                    input_granules = umm.get("InputGranules", [])
+                    raw_hls = next(
+                        (g for g in input_granules if g.startswith("HLS.")), "N/A"
+                    )
+
+                    if raw_hls != "N/A":
+                        parts = raw_hls.split(".")
+                        hls_granule_id = (
+                            ".".join(parts[:6]) if len(parts) >= 6 else raw_hls
+                        )
+                        hls_links = fetch_hls_granule_links(hls_granule_id)
+
+                    # Fallback: Search CMR dynamically via Tile ID and Date (Required for DIST)
+                    else:
+                        try:
+                            # Extract Tile ID (e.g., T11SLT) from the OPERA GranuleUR
+                            parts = granule_id.split("_")
+                            tile_id = next(
+                                (p for p in parts if p.startswith("T") and len(p) == 6),
+                                None,
+                            )
+
+                            if tile_id and start_time != "N/A":
+                                date_only = start_time.split("T")[0]
+                                search_bounds = geom.bounds if geom else None
+
+                                # Query CMR for all HLS granules on that day over the bounding box
+                                hls_results = earthaccess.search_data(
+                                    short_name=["HLSS30", "HLSL30"],
+                                    temporal=(
+                                        f"{date_only}T00:00:00",
+                                        f"{date_only}T23:59:59",
+                                    ),
+                                    bounding_box=search_bounds,
+                                )
+
+                                # Filter the results to find the one matching the exact Tile ID
+                                for r in hls_results:
+                                    g_name = r.get("umm", {}).get("GranuleUR", "")
+                                    if f".{tile_id}." in g_name:
+                                        hls_links = r.data_links()
+                                        hls_granule_id = ".".join(g_name.split(".")[:6])
+                                        break
+                        except Exception as e:
+                            LOGGER.warning(
+                                f"Failed to dynamically locate HLS granule for {granule_id}: {e}"
+                            )
+
+                    # Map bands based on HLSS30 or HLSL30 naming conventions
+                    for href in hls_links:
+                        if href.endswith(".tif"):
+                            if "B04" in href or "band04" in href.lower():
+                                hls_red = href
+                            elif "B03" in href or "band03" in href.lower():
+                                hls_green = href
+                            elif "B02" in href or "band02" in href.lower():
+                                hls_blue = href
+                            elif "S30" in hls_granule_id and (
+                                "B8A" in href or "band8a" in href.lower()
+                            ):
+                                hls_nir = href
+                            elif "L30" in hls_granule_id and (
+                                "B05" in href or "band05" in href.lower()
+                            ):
+                                hls_nir = href
+                            elif "Fmask" in href or "fmask" in href.lower():
+                                hls_fmask = href
+
+                # Appends all 6 elements to keep data rows and headers perfectly 1-to-1
+                row_data.extend(
+                    [hls_granule_id, hls_red, hls_green, hls_blue, hls_nir, hls_fmask]
+                )
+
+            ws.append(row_data)
 
         if compute_cloudiness and overall_area > 0:
             overall_cloud_cover_percent = 100.0 * (overall_cloudy_area / overall_area)