Classify dep staging

Makefile

@@ -58,7 +58,6 @@ geomad-test-ausp:
 		--year 2000 \
 		--version 0-3-0-test \
 		--decimated \
-		--no-single-region \
 		--bucket data.ldn.auspatious.com \
 		--overwrite;
 geomad-test-dep-staging:
@@ -69,7 +68,6 @@ geomad-test-dep-staging:
 		--version 0-3-0-test \
 		--collection-url-root="https://stac.staging.digitalearthpacific.io/collections" \
 		--decimated \
-		--single-region \
 		--bucket dep-public-staging \
 		--overwrite;
 
@@ -115,73 +113,76 @@ collection-geomad-test-dep-staging:
 
 
 
-# #### Training Data
-# training-data-generate:
-# 	for site in $(PACIFIC_TRAINING_TILES); do \
-# 		tile_id=$$(echo $$site | cut -d: -f1); \
-# 		region=$$(echo $$site | cut -d: -f2); \
-# 		country_name=$$(echo $$site | cut -d: -f3 | tr '_' ' '); \
-# 		country_code=$$(echo $$site | cut -d: -f4); \
-# 		ldn training generate-training-data \
-# 			--tile-id $$tile_id \
-# 			--region $$region \
-# 			--country-name "$$country_name" \
-# 			--country-code "$$country_code"; \
-# 	done;
-
-# # poetry run ldn training generate-training-data \
-# # 	--tile-id 028_030 \
-# # 	--region pacific \
-# # 	--country-name "Papua New Guinea" \
-# # 	--country-code "PNG" \
-# # 	--geomad-version 0-2-1;
+#### Training Data
+# Geomad version: 0-2-1 in DEP staging, 0-3-0 in DEP public.
+training-data-generate:
+	for site in $(PACIFIC_TRAINING_TILES) do \
+		tile_id=$$(echo $$site | cut -d: -f1); \
+		region=$$(echo $$site | cut -d: -f2); \
+		country_name=$$(echo $$site | cut -d: -f3 | tr '_' ' '); \
+		country_code=$$(echo $$site | cut -d: -f4); \
+		ldn training generate-training-data \
+			--tile-id $$tile_id \
+			--region $$region \
+			--country-name "$$country_name" \
+			--country-code "$$country_code" \
+			--geomad-version 0-2-1 \
+			--geomad-bucket dep-public-staging \
+			--output-bucket dep-public-staging \
+			--single-region \
+			--product-owner dep \
+			--no-overwrite; \
+	done;
 
+#### Make the model using ldn-lulc/notebooks/1_Train_Model.ipynb
 
 
 
 # ###### LULC Classification/Prediction
 
 # # Predict LULC for the test tiles and one year (2025).
 
-# # 1. Print tasks
+# # Print tasks
 # print-tasks-lulc-2020:
 # 	ldn print-tasks \
 # 	--years="2020" \
 # 	--region="pacific" \
 # 	--dataset="lulc";
 
 
-# # 2. Classify
-# predict-lulc-test-tiles-2020:
-# 	for site in $(TEST_TILES); do \
-# 		tile_id=$${site%%:*}; \
-# 		region=$${site#*:}; region=$${region%%:*}; \
-# 		ldn lulc run \
-# 			--tile-id $$tile_id \
-# 			--year 2020 \
-# 			--version $(LULC_VERSION) \
-# 			--geomad-version $(GEOMAD_VERSION) \
-# 			--region $$region \
-# 			$(DECIMATED) \
-# 			--overwrite; \
-# 	done;
-
-# lulc-2-regions-decimated:
-# 	for site in $(TEST_TILES_2_REGIONS); do \
-# 		tile_id=$${site%%:*}; \
-# 		region=$${site#*:}; region=$${region%%:*}; \
-# 		ldn lulc run \
-# 			--tile-id $$tile_id \
-# 			--year 2010 \
-# 			--version $(LULC_VERSION) \
-# 			--geomad-version $(GEOMAD_VERSION) \
-# 			--region $$region \
-# 			--decimated \
-# 			--overwrite; \
-# 	done;
-
-
-# # 3. Update the STAC-Geoparquet index after all tiles/years have run.
+# # Classify
+lulc-predict-test:
+	ldn lulc run \
+		--tile-id 028_030 \
+		--year 2000 \
+		--region pacific \
+		--version 0-0-9 \
+		--geomad-version 0-2-1 \
+		--bucket dep-public-staging \
+		--product-owner dep \
+		--model-path="/Users/wj/Projects/ldn-lulc/ldn-lulc/ldn/models/0-0-9/pacific/2020/lulc_random_forest_model_pacific_2020.joblib" \
+		--no-overwrite;
+
+# 		--model-path="https://dep-public-staging.s3.us-west-2.amazonaws.com/dep_ls_lulc/models/0-0-9/pacific/2020/lulc_random_forest_model_pacific_2020.joblib" \
+
+lulc-predict-test-2:
+	for site in $(PACIFIC_TRAINING_TILES) do \
+		tile_id=$$(echo $$site | cut -d: -f1); \
+		region=$$(echo $$site | cut -d: -f2); \
+		for year in 2000 2025; do \
+			ldn lulc run \
+				--tile-id $$tile_id \
+				--year $$year \
+				--region $$region \
+				--version 0-0-9 \
+				--geomad-version 0-2-1 \
+				--bucket dep-public-staging \
+				--product-owner dep \
+				--model-path="/Users/wj/Projects/ldn-lulc/ldn-lulc/ldn/models/0-0-9/pacific/2020/lulc_random_forest_model_pacific_2020.joblib" \
+				--no-overwrite;
+		done;
+	done;
+# # Update the STAC-Geoparquet index after all tiles/years have run.
 # index-lulc:
 # 	ldn index-to-stac-geoparquet \
 # 	--dataset "lulc" \

ldn/cli_geomad.py

@@ -85,10 +85,6 @@ def run(
     year: Annotated[str, typer.Option()],
     version: Annotated[str, typer.Option()],
     region: Annotated[Literal["pacific", "non-pacific"], typer.Option()],
-    single_region: Annotated[
-        bool,
-        typer.Option(help="Whether to use the single region prefix for the collection_url_root (e.g. 'dep_ls_geomad')"),
-    ],
     bucket: Annotated[str | None, typer.Option(help="S3 bucket for data.")] = None,
     product_owner: Annotated[str | None, typer.Option(help="Override the region-derived owner prefix.")] = None,
     overwrite: Annotated[bool, typer.Option()] = False,
@@ -112,7 +108,7 @@ def run(
         str | None,
         typer.Option(
             help="Override the default collection URL root"
-            " e.g for a STAC API like 'https://stac.digitalearthpacific.org/collections/dep_ls_geomad'"
+            " e.g for a STAC API like 'https://stac.digitalearthpacific.org/collections'"
         ),
     ] = None,
     sensor: Annotated[str, typer.Option(help="Sensor name, e.g. 'ls'.")] = SENSOR,
@@ -221,6 +217,7 @@ def run(
         overwrite,
         collection_url_root=collection_url_root,
         s3_client=s3_client,
+        sensor=sensor,
     )
     if components is None:
         return  # Skip due to no overwrite.

ldn/cli_lulc.py

@@ -7,10 +7,9 @@
 from ldn.utils import (
     GEOMAD_VERSION,
     LULC_VERSION,
-    MODEL_VERSION,
+    SENSOR,
     LdnError,
     get_env_var,
-    owner_for_region,
 )
 
 classify_app = typer.Typer()
@@ -36,7 +35,7 @@ def run(
     product_owner: str | None = typer.Option(None, help="Override the region-derived owner prefix."),
     model_path: str = typer.Option(
         # TODO: defaults to pacific. Later have per region/time period models.
-        f"https://s3.us-west-2.amazonaws.com/data.ldn.auspatious.com/models/{MODEL_VERSION}/pacific/2020/lulc_random_forest_model_pacific_2020.joblib",
+        "https://dep-public-staging.s3.us-west-2.amazonaws.com/dep_ls_lulc/models/0-0-9/pacific/2020/lulc_random_forest_model_pacific_2020.joblib",
         help="Model to use for LULC classification.",
     ),
     decimated: bool = typer.Option(
@@ -66,18 +65,19 @@ def run(
             help="Chunk size in pixels for x and y dimensions. Larger chunk sizes may be faster but use more memory."
         ),
     ] = 1024,
-    single_region: bool = typer.Option(
-        ...,
-        help="Whether to use the single region prefix (e.g. 'dep_ls_geomad') "
-        "or the generic prefix (e.g. 'ls_geomad') when accessing GeoMAD data.",
-    ),
+    sensor: str = typer.Option(SENSOR, help="Sensor to use for LULC classification. Defaults to 'ls'."),
+    collection_url_root: Annotated[
+        str | None,
+        typer.Option(
+            help="Override the default collection URL root"
+            " e.g for a STAC API like 'https://stac.digitalearthpacific.org/collections'"
+        ),
+    ] = None,
 ) -> None:
     if int(year) < 2000 or int(year) > 2025:
         raise LdnError("Year must be between 2000 and 2025.")
 
     bucket = bucket or get_env_var("BUCKET")  # Default
-    owner = owner_for_region(region, product_owner)
-    # TODO: Use build_prefix() here?
 
     run_classify_task(
         tile_id,
@@ -86,7 +86,7 @@ def run(
         geomad_version=geomad_version,
         region=region,
         bucket=bucket,
-        owner=owner,
+        product_owner=product_owner,
         model_path=model_path,
         xy_chunk_size=xy_chunk_size,
         decimated=decimated,
@@ -97,5 +97,6 @@ def run(
         memory_limit=memory_limit,
         n_workers=n_workers,
         threads_per_worker=threads_per_worker,
-        single_region=single_region,
+        sensor=sensor,
+        collection_url_root=collection_url_root,
     )

ldn/lulc.py

@@ -22,7 +22,7 @@
 from sklearn.ensemble import RandomForestClassifier
 from typing_extensions import Annotated
 
-from ldn.aws import configure_s3_access_profile
+from ldn.aws import configure_s3_access_profile, s3_client
 from ldn.geomad import AwsStacTask as Task
 from ldn.grids import get_gridspec
 from ldn.raster import (
@@ -33,15 +33,19 @@
     scale_offset_landsat,
 )
 from ldn.utils import (
+    GEOMAD_DATASET_ID,
     GEOMAD_VERSION,
     LULC_DATASET_ID,
     LULC_VERSION,
     SOURCE_COOP_PREFIX_LULC,
     WGS84,
     LdnError,
     get_analysis_epsg,
+    get_public_url_base,
+    get_stac_geoparquet_key,
     get_stac_geoparquet_url,
     is_bucket_source_coop,
+    owner_for_region,
     parse_tile_id,
 )
 
@@ -206,7 +210,13 @@ def do_prediction(
     nodata_mask = stacked.isnull().any(dim="variable")
 
     # Build observation table: fill NaN with nodata_value (masked pixels are excluded below).
-    obs = stacked.squeeze().fillna(nodata_value).transpose().to_dataframe()
+    obs = (
+        stacked.squeeze()
+        .fillna(nodata_value)
+        .transpose()
+        .to_dataset(dim="variable")  # pivot "variable" into named columns
+        .to_dataframe()
+    )
 
     # Validate that all model features are present before reindexing.
     missing = set(model.feature_names_in_) - set(obs.columns)
@@ -382,7 +392,7 @@ def run_classify_task(
     geomad_version: Annotated[str, typer.Option()],
     region: Literal["pacific", "non-pacific"],
     bucket: str,
-    owner: str,
+    product_owner: str | None,
     model_path: str,
     xy_chunk_size: int,
     decimated: bool,
@@ -393,7 +403,8 @@ def run_classify_task(
     memory_limit: str,
     n_workers: int,
     threads_per_worker: int,
-    single_region: bool,
+    sensor: str,
+    collection_url_root: str | None,
 ) -> None:
     """Run LULC prediction for a single tile and year, writing results to S3.
 
@@ -407,7 +418,7 @@ def run_classify_task(
         geomad_version: Version of the GeoMAD data to use (e.g. "0-0-1").
         region: Grid region, either "pacific" or "non-pacific".
         bucket: S3 bucket for output COGs, STAC metadata, and input GeoMAD source data.
-        owner: Output prefix for paths (e.g. "dep" or "ci" or owner override).
+        product_owner: Override the region-derived owner prefix.
         model_path: Path or URL to the trained joblib model.
         xy_chunk_size: Chunk size in pixels for lazy loading.
         decimated: If True, use 10x lower resolution (for testing).
@@ -418,6 +429,9 @@ def run_classify_task(
         memory_limit: Per-worker Dask memory limit.
         n_workers: Number of Dask workers.
         threads_per_worker: Number of threads per Dask worker.
+        single_region: If True, use the single region prefix (e.g. 'dep_ls_geomad') for GeoMAD data.
+        sensor: Sensor to use for LULC classification e.g. 'ls'.
+        collection_url_root: Override the default collection URL root for STAC metadata.
     """
     logger.info(f"Starting processing. Tile ID: {tile_id}, Year: {year}, Region: {region}, Version: {version}.")
     logger.info(
@@ -433,11 +447,11 @@ def run_classify_task(
         logger.info(
             f"Overriding the latest LULC prediction version ({LULC_VERSION}) with the specified version ({version})."
         )
+    owner = owner_for_region(region, product_owner)
 
-    # geomad_stac_geoparquet_key = get_stac_geoparquet_key(
-    #     bucket, single_region, product_owner, sensor, "geomad", geomad_version
-    # )
-    geomad_stac_geoparquet_key = ""  # TODO: Fix this.
+    geomad_stac_geoparquet_key = get_stac_geoparquet_key(
+        bucket, product_owner, sensor, GEOMAD_DATASET_ID, geomad_version
+    )
     geomad_stac_geoparquet_url = get_stac_geoparquet_url(bucket, geomad_stac_geoparquet_key)
 
     tile_id_tuple = parse_tile_id(tile_id)
@@ -464,18 +478,30 @@ def run_classify_task(
     logger.info("Loading model")
     loaded_model = _load_joblib_model(model_path)
 
+    # TODO: Could this block be a function because it will be done in cli_collection.py?
+    _is_source_coop = is_bucket_source_coop(bucket)
+    public_url = get_public_url_base(bucket)
+    if _is_source_coop:
+        public_url = f"{public_url}/{SOURCE_COOP_PREFIX_LULC}"
+    collection_url_root = collection_url_root or f"{public_url}/collections"
+
+    owner = owner_for_region(region, product_owner)
+
     components = build_pipeline_components(
         tile_id_tuple,
         year,
         version,
         bucket,
-        owner,
+        owner,  # This owner respects single_region because geomad always writes to an owner.
         LULC_DATASET_ID,
-        SOURCE_COOP_PREFIX_LULC if is_bucket_source_coop(bucket) else None,
+        SOURCE_COOP_PREFIX_LULC if _is_source_coop else None,
         overwrite,
+        collection_url_root=collection_url_root,
+        s3_client=s3_client,
+        sensor=sensor,
     )
     if components is None:
-        return  # Task exists and overwrite is False, so skipping processing.
+        return  # Skip due to no overwrite.
     itempath, stac_creator, writer = components
 
     searcher = StacGeoparquetSearcher(

ldn/raster.py

@@ -19,7 +19,6 @@
 from shapely.geometry import box
 
 from ldn.utils import (
-    SENSOR,
     WGS84,
     LdnError,
     get_public_url_base,
@@ -308,6 +307,7 @@ def build_pipeline_components(
     overwrite: bool,
     collection_url_root: str,
     s3_client: BaseClient,
+    sensor: str,
 ) -> tuple[PrefixedS3ItemPath, StacCreator, AwsDsCogWriter] | None:
     """Build shared pipeline components for GeoMAD and classify tasks.
 
@@ -319,7 +319,7 @@ def build_pipeline_components(
         key_prefix=source_coop_prefix,
         prefix=owner,
         bucket=bucket,
-        sensor=SENSOR,
+        sensor=sensor,
         dataset_id=dataset_id,
         version=version,
         time=year,