global_geomorphometric_layers.py

# ---
# title:   Global Geomorphometric Layers (Geomorpho90m)
# author:  Brendan Casey
# created: 2026-07-10
# inputs:
#   - Geomorpho90m ImageCollections
#     (projects/sat-io/open-datasets/Geomorpho90m)
#   - AB2020 provincial boundary (EE asset)
# outputs:
#   - Multiband Geomorpho90m GeoTIFF for Alberta, at native
#     (~90 m) or on the ABMI 1 km reference grid, per
#     EXPORT_TARGET (exported to Google Drive). The 1 km
#     product drops the raw 'aspect' band (see notes).
# notes:
#   This script loads multiple geomorphometric variables
#   from the Geomorpho90m dataset, mosaics and clips them
#   to the AOI, and combines them into a single multiband
#   image. Map visualization layers from the original GEE
#   JavaScript are dropped.
#
#   Citation:
#   Amatulli, G., McInerney, D., Sethi, T., Strobl, P.,
#   Domisch, S. (2020). Geomorpho90m, empirical evaluation
#   and accuracy assessment of global high-resolution
#   geomorphometric layers. Scientific Data 7(1), 1-18.
#
#   Setup (once):
#     pip install earthengine-api
#     earthengine authenticate
#   Then set EE_PROJECT in _gee_config.py to your
#   registered Earth Engine cloud project and run the
#   script.
# ---

import os
import sys

import ee

# Make utils importable regardless of the working
# directory VS Code runs the script from
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

from _gee_config import DRIVE_FOLDER, PROVINCIAL_BOUNDARY_ASSET
from utils.compute_report import ComputeReport
from utils.gee_utils import (
    export_image_to_drive,
    export_to_reference_grid,
    initialize_ee,
)

# 1. Setup ----

# 1.1 User parameters ----
EXPORT_SCALE = 90  # meters (Geomorpho90m native ~90 m)
EXPORT_CRS = "EPSG:3400"  # AB 10-TM (Forest)
# Raster export target. "native" exports the full ~90 m stack;
# "reference_grid" aggregates (area mean) onto the ABMI 1 km
# grid so it stacks with the other 1 km covariates. Geomorpho90m
# is a stored (pyramided) dataset, so 90 m -> 1 km is well under
# Earth Engine's per-tile reprojection limit.
EXPORT_TARGET = "reference_grid"  # "native" or "reference_grid"

# Compute ring grown around the aoi before the source is
# clipped, sized at 2x the output scale. Every output pixel -
# a 1 km grid cell or a native pixel - is then built from a
# full neighbourhood rather than one truncated at the aoi
# edge; a 1 km cell can touch the aoi at a corner and still
# reach a full diagonal (1414 m) beyond it. The exported
# image is clipped back to the plain aoi, so the ring never
# widens the output.
COARSE_SCALE = 1000  # ABMI reference grid cell (m)
AGG_BUFFER_M = 2 * (
    COARSE_SCALE
    if EXPORT_TARGET == "reference_grid"
    else EXPORT_SCALE
)
BUFFER_MAX_ERROR_M = 100
PRINT_STATS = True  # min/max check (slow for large AOIs)
USE_TEST_AOI = True  # True: small test AOI; False: Alberta
COMPUTE_REPORT = True  # write EECU usage report (txt)
# Block until every export task finishes so its batch
# EECU-seconds land in the compute report. Costs the full
# export runtime (hours for a province-wide run), so keep it
# False for production runs and turn it on when profiling a
# test AOI.
WAIT_FOR_EXPORTS = False

# Base path and Geomorpho90m collections to combine, in the
# order they are stacked into the multiband export image.
BASE_PATH = "projects/sat-io/open-datasets/Geomorpho90m/"
COLLECTION_NAMES = [
    "aspect",           # Aspect
    "aspect-cosine",    # Aspect-Cosine
    "aspect-sine",      # Aspect-Sine
    "convergence",      # Convergence Index
    "cti",              # Compound Topographic Index (CTI)
    "dev-magnitude",    # Deviation Magnitude
    "dev-scale",        # Deviation Scale
    "eastness",         # Eastness
    "elev-stdev",       # Elevation Standard Deviation
    "northness",        # Northness
    "rough-magnitude",  # Multiscale Roughness Magnitude
    "rough-scale",      # Multiscale Roughness Scale
    "roughness",        # Roughness
    "slope",            # Slope
    "spi",              # Stream Power Index
    "tpi",              # Topographic Position Index (TPI)
    "tri",              # Terrain Ruggedness Index (TRI)
    "vrm",              # Vector Ruggedness Measure (VRM)
]

# 1.2 Initialize Earth Engine ----
# Project ID is read from _gee_config.py
initialize_ee()

# 1.3 Set up compute usage report ----
# Profiles EECU usage per section and per export task.
# Best used with USE_TEST_AOI = True to find choke
# points cheaply before a full-province run.
report = ComputeReport(
    "global_geomorphometric_layers",
    enabled=COMPUTE_REPORT,
)

# 2. Define study area ----
# This section defines the export geometry. It uses a
# small test polygon when USE_TEST_AOI is True; otherwise
# it uses the AB2020 provincial boundary asset.

if USE_TEST_AOI:
    # Small aoi for testing purposes
    aoi = ee.Geometry.Polygon([
        [-113.5, 55.5],  # Top-left corner
        [-113.5, 55.0],  # Bottom-left corner
        [-112.8, 55.0],  # Bottom-right corner
        [-112.8, 55.5],  # Top-right corner
    ])
else:
    aoi = ee.FeatureCollection(
        PROVINCIAL_BOUNDARY_ASSET
    ).geometry()

# 3. Geomorpho90m processing ----
# This section loads, mosaics, clips, and renames each
# Geomorpho90m collection, then combines them into a single
# multiband image.


def load_and_process(collection_name, aoi):
    """Load, mosaic, clip, and rename a collection.

    Parameters
    ----------
    collection_name : str
        Geomorpho90m collection short name.
    aoi : ee.Geometry
        Area of interest to clip to.

    Returns
    -------
    ee.Image
        Single-band image renamed to ``collection_name``.
    """
    return (
        ee.ImageCollection(BASE_PATH + collection_name)
        .mosaic()
        .clip(aoi)
        .rename(collection_name)
    )


# Aggregation reads from the ring (AGG_BUFFER_M); the 1 km
# result is clipped back to the plain aoi downstream.
clip_geom = (
    aoi.buffer(AGG_BUFFER_M, BUFFER_MAX_ERROR_M)
    if AGG_BUFFER_M
    else aoi
)

geomorpho90m = load_and_process(COLLECTION_NAMES[0], clip_geom)
for name in COLLECTION_NAMES[1:]:
    geomorpho90m = geomorpho90m.addBands(
        load_and_process(name, clip_geom)
    )

# 3.1 Check min and max values (optional) ----
# Also runs when COMPUTE_REPORT is on: Earth Engine is
# lazy, so the profiler needs an evaluated computation
# (getInfo) to measure per-algorithm EECU usage.
if PRINT_STATS or COMPUTE_REPORT:
    with report.section("Geomorpho90m min/max (reduceRegion)"):
        stats = geomorpho90m.reduceRegion(
            reducer=ee.Reducer.minMax(),
            geometry=aoi,
            scale=EXPORT_SCALE,
            maxPixels=1e13,
            bestEffort=True,
        ).getInfo()
    print("Geomorpho90m min and max values:", stats)

# 4. Export data ----
# Export at the target chosen by EXPORT_TARGET. "native" writes
# the full ~90 m stack; "reference_grid" aggregates (area mean)
# onto the ABMI 1 km grid. The raw 'aspect' band is dropped from
# the 1 km product because averaging a circular angle (0-360
# deg) is meaningless; 'aspect-cosine' and 'aspect-sine' carry
# aspect correctly and mean cleanly. setDefaultProjection pins
# the native base for reduceResolution. Set wait=True to block;
# otherwise monitor at https://code.earthengine.google.com/tasks

if EXPORT_TARGET == "reference_grid":
    grid_bands = [n for n in COLLECTION_NAMES if n != "aspect"]
    geomorpho90m_grid = geomorpho90m.select(
        grid_bands
    ).setDefaultProjection(crs=EXPORT_CRS, scale=EXPORT_SCALE)
    task = export_to_reference_grid(
        image=geomorpho90m_grid,
        aoi=aoi,
        description="Geomorpho90m_AB_abmi1km",
        folder=DRIVE_FOLDER,
        file_name_prefix="global_geomorphometric_layers_abmi1km",
        aggregate=True,
        wait=False,
    )
elif EXPORT_TARGET == "native":
    task = export_image_to_drive(
        image=geomorpho90m.clip(aoi),
        description="Geomorpho90m_AB_native",
        region=aoi,
        folder=DRIVE_FOLDER,
        file_name_prefix="global_geomorphometric_layers_native",
        scale=EXPORT_SCALE,
        crs=EXPORT_CRS,
        max_pixels=1e13,
        wait=False,
    )
else:
    raise ValueError(
        "Unknown EXPORT_TARGET: "
        f"{EXPORT_TARGET!r} (use 'native' or 'reference_grid')"
    )

# 5. Compute usage report ----
# This section waits for the export to finish, records
# its total EECU-seconds, and writes the txt report to
# gee_compute_reports/. Note: a full-province export can
# take hours; for a quick profile use the test AOI.

if WAIT_FOR_EXPORTS:
    report.log_task(task)
report.write()

# End of script ----