Task/tpi summarize hfi

.vscode/settings.json

@@ -54,19 +54,25 @@
     ],
     "typescript.preferences.importModuleSpecifier": "non-relative",
     "cSpell.words": [
+        "actuals",
         "aiobotocore",
         "Albers",
         "allclose",
         "anyio",
         "APCP",
         "botocore",
         "cffdrs",
+        "cutline",
         "excinfo",
+        "FBAN",
         "firezone",
         "GDPS",
+        "geoalchemy",
         "GEOGCS",
+        "geopackage",
         "geospatial",
         "geotiff",
+        "gpkg",
         "grib",
         "gribs",
         "hourlies",
@@ -75,7 +81,9 @@
         "maxx",
         "maxy",
         "miny",
+        "morecast",
         "ndarray",
+        "Neighbour",
         "numba",
         "osgeo",
         "PMTILES",
@@ -86,6 +94,7 @@
         "PRIMEM",
         "PROJCS",
         "RDPS",
+        "reproject",
         "rocketchat",
         "sfms",
         "sqlalchemy",

api/alembic/versions/6910d017b626_add_advisory_tpi_table.py

@@ -0,0 +1,54 @@
+"""add advisory tpi table
+
+Revision ID: 6910d017b626
+Revises: be128a7bb4fd
+Create Date: 2024-07-31 16:27:31.642156
+
+"""
+
+from alembic import op
+import sqlalchemy as sa
+from sqlalchemy.dialects import postgresql
+
+# revision identifiers, used by Alembic.
+revision = "6910d017b626"
+down_revision = "be128a7bb4fd"
+branch_labels = None
+depends_on = None
+
+
+def upgrade():
+    # ### commands auto generated by Alembic ###
+    op.create_table(
+        "advisory_tpi_stats",
+        sa.Column("id", sa.Integer(), nullable=False),
+        sa.Column("advisory_shape_id", sa.Integer(), nullable=False),
+        sa.Column("run_parameters", sa.Integer(), nullable=False),
+        sa.Column("valley_bottom", sa.Integer(), nullable=False),
+        sa.Column("mid_slope", sa.Integer(), nullable=False),
+        sa.Column("upper_slope", sa.Integer(), nullable=False),
+        sa.Column("pixel_size_metres", sa.Integer(), nullable=False),
+        sa.ForeignKeyConstraint(
+            ["advisory_shape_id"],
+            ["advisory_shapes.id"],
+        ),
+        sa.ForeignKeyConstraint(
+            ["run_parameters"],
+            ["run_parameters.id"],
+        ),
+        sa.PrimaryKeyConstraint("id"),
+        comment="Elevation TPI stats per fire shape",
+    )
+    op.create_index(op.f("ix_advisory_tpi_stats_advisory_shape_id"), "advisory_tpi_stats", ["advisory_shape_id"], unique=False)
+    op.create_index(op.f("ix_advisory_tpi_stats_id"), "advisory_tpi_stats", ["id"], unique=False)
+    op.create_index(op.f("ix_advisory_tpi_stats_run_parameters"), "advisory_tpi_stats", ["run_parameters"], unique=False)
+    # ### end Alembic commands ###
+
+
+def downgrade():
+    # ### commands auto generated by Alembic ###
+    op.drop_index(op.f("ix_advisory_tpi_stats_run_parameters"), table_name="advisory_tpi_stats")
+    op.drop_index(op.f("ix_advisory_tpi_stats_id"), table_name="advisory_tpi_stats")
+    op.drop_index(op.f("ix_advisory_tpi_stats_advisory_shape_id"), table_name="advisory_tpi_stats")
+    op.drop_table("advisory_tpi_stats")
+    # ### end Alembic commands ###

api/app/.env.example

@@ -60,4 +60,5 @@ OBJECT_STORE_SECRET=object_store_secret
 OBJECT_STORE_BUCKET=object_store_bucket
 DEM_NAME=dem_mosaic_250_max.tif
 TPI_DEM_NAME=bc_dem_50m_tpi.tif
+CLASSIFIED_TPI_DEM_NAME=bc_dem_50m_tpi_win100_classified.tif
 SENTRY_DSN=some_dsn

api/app/auto_spatial_advisory/elevation.py

@@ -1,10 +1,12 @@
 """Takes a classified HFI image and calculates basic elevation statistics associated with advisory areas per fire zone."""
 
+from dataclasses import dataclass
 from datetime import date, datetime
 from time import perf_counter
 import logging
 import os
 import tempfile
+from typing import Dict
 import numpy as np
 from osgeo import gdal
 from sqlalchemy.ext.asyncio import AsyncSession
@@ -13,16 +15,47 @@
 from app import config
 from app.auto_spatial_advisory.classify_hfi import classify_hfi
 from app.auto_spatial_advisory.run_type import RunType
-from app.db.crud.auto_spatial_advisory import get_run_parameters_id, save_advisory_elevation_stats
+from app.db.crud.auto_spatial_advisory import get_run_parameters_id, save_advisory_elevation_stats, save_advisory_elevation_tpi_stats
 from app.db.database import get_async_read_session_scope, get_async_write_session_scope, DB_READ_STRING
-from app.db.models.auto_spatial_advisory import AdvisoryElevationStats
+from app.db.models.auto_spatial_advisory import AdvisoryElevationStats, AdvisoryTPIStats
+from app.auto_spatial_advisory.hfi_filepath import get_raster_filepath, get_raster_tif_filename
 from app.utils.s3 import get_client
+from app.utils.geospatial import GeospatialOptions, WarpOptions, cut_raster_by_shape_id, data_array_to_raster, read_raster_into_memory
 
 
 logger = logging.getLogger(__name__)
 DEM_GDAL_SOURCE = None
 
 
+async def process_elevation_tpi(run_type: RunType, run_datetime: datetime, for_date: date):
+    """
+    Create new elevation statistics records for the given parameters.
+
+    :param hfi_s3_key: the object store key pointing to the hfi tif to intersect with tpi layer
+    :param run_type: The type of run to process. (is it a forecast or actual run?)
+    :param run_datetime: The date and time of the run to process. (when was the hfi file created?)
+    :param for_date: The date of the hfi to process. (when is the hfi for?)
+    """
+    logger.info("Processing elevation stats %s for run date: %s, for date: %s", run_type, run_datetime, for_date)
+    perf_start = perf_counter()
+    # Get the id from run_parameters associated with the provided run_type, for_date and for_datetime
+    async with get_async_read_session_scope() as session:
+        run_parameters_id = await get_run_parameters_id(session, run_type, run_datetime, for_date)
+
+        stmt = select(AdvisoryTPIStats).where(AdvisoryTPIStats.run_parameters == run_parameters_id)
+
+        exists = (await session.execute(stmt)).scalars().first() is not None
+        if not exists:
+            fire_zone_stats = await process_tpi_by_firezone(run_type, run_datetime.date(), for_date)
+            await store_elevation_tpi_stats(session, run_parameters_id, fire_zone_stats)
+        else:
+            logger.info("Elevation stats already computed")
+
+    perf_end = perf_counter()
+    delta = perf_end - perf_start
+    logger.info("%f delta count before and after processing elevation stats", delta)
+
+
 async def process_elevation(source_path: str, run_type: RunType, run_datetime: datetime, for_date: date):
     """Create new elevation statistics records for the given parameters.
 
@@ -178,6 +211,75 @@ def apply_threshold_mask_to_dem(threshold: int, mask_path: str, temp_dir: str):
     return masked_dem_path
 
 
+@dataclass(frozen=True)
+class FireZoneTPIStats:
+    """
+    Captures fire zone stats of TPI pixels hitting >4K HFI threshold via
+    a dictionary, fire_zone_stats, of {source_identifier: {1: X, 2: Y, 3: Z}}, where 1 = valley bottom, 2 = mid slope, 3 = upper slope
+    and X, Y, Z are pixel counts at each of those elevation classes respectively.
+
+    Also includes the TPI raster's pixel size in metres.
+    """
+
+    fire_zone_stats: Dict[int, Dict[int, int]]
+    pixel_size_metres: int
+
+
+async def process_tpi_by_firezone(run_type: RunType, run_date: date, for_date: date) -> FireZoneTPIStats:
+    """
+    Given run parameters, lookup associated snow-masked HFI and static classified TPI geospatial data.
+    Cut out each fire zone shape from the above and intersect the TPI and HFI pixels, counting each pixel contributing to the TPI class.
+    Capture all fire zone stats keyed by its source_identifier.
+
+    :param run_parameters_id: The RunParameter object id associated with this run_type, for_date and run_datetime
+    :return: fire zone TPI status
+    """
+    raster_options = GeospatialOptions(include_geotransform=True, include_projection=True, include_x_size=True, include_y_size=True)
+    fire_zone_stats: Dict[int, Dict[int, int]] = {}
+
+    hfi_raster_filename = get_raster_tif_filename(for_date)
+    hfi_raster_key = get_raster_filepath(run_date, run_type, hfi_raster_filename)
+
+    async with get_client() as (client, bucket):
+        tpi_result = await read_raster_into_memory(
+            client,
+            bucket,
+            f'dem/tpi/{config.get("CLASSIFIED_TPI_DEM_NAME")}',
+            raster_options,
+        )
+        output_options = GeospatialOptions(
+            warp_options=WarpOptions(
+                source_geotransform=tpi_result.data_geotransform,
+                source_projection=tpi_result.data_projection,
+                source_x_size=tpi_result.data_x_size,
+                source_y_size=tpi_result.data_y_size,
+            )
+        )
+        hfi_result = await read_raster_into_memory(client, bucket, hfi_raster_key, output_options)
+
+        # We only want pixels with HFI >4k as values to intersect with TPI raster
+        hfi_result.data_array = np.where(hfi_result.data_array >= 1, 1, 0)
+        hfi_masked_tpi = np.multiply(tpi_result.data_array, hfi_result.data_array)
+
+        hfi_masked_tpi_ds = data_array_to_raster(hfi_masked_tpi, hfi_raster_key, output_options)
+        async with get_async_write_session_scope() as session:
+            stmt = text("SELECT id, source_identifier FROM advisory_shapes;")
+            result = await session.execute(stmt)
+            rows = result.all()
+
+            for row in rows:
+                cut_hfi_masked_tpi = cut_raster_by_shape_id(row[0], row[1], hfi_masked_tpi_ds, raster_options)
+                # Get unique values and their counts
+                tpi_classes, counts = np.unique(cut_hfi_masked_tpi.data_array, return_counts=True)
+                tpi_class_freq_dist = dict(zip(tpi_classes, counts))
+
+                # Drop TPI class 4, this is the no data value from the TPI raster
+                tpi_class_freq_dist.pop(4, None)
+                fire_zone_stats[row[1]] = tpi_class_freq_dist
+
+            return FireZoneTPIStats(fire_zone_stats=fire_zone_stats, pixel_size_metres=tpi_result.data_geotransform[1])
+
+
 async def process_elevation_by_firezone(threshold: int, masked_dem_path: str, run_parameters_id: int):
     """
     Given a tif that only contains elevations values at pixels where HFI exceeds the threshold, calculate statistics
@@ -205,7 +307,7 @@ def intersect_raster_by_firezone(threshold: int, advisory_shape_id: int, source_
 
     :param threshold: The current threshold being processed, 1 = 4k-10k, 2 = > 10k
     :param advisory_shape_id: The id of the fire zone (aka advisory_shape object) to clip with
-    :param source_identifier: The source identifer of the fire zone.
+    :param source_identifier: The source identifier of the fire zone.
     :param raster_path: The path to the raster to be clipped.
     :param temp_dir: A temporary location for storing intermediate files
     """
@@ -261,3 +363,26 @@ async def store_elevation_stats(session: AsyncSession, threshold: int, shape_id:
         threshold=threshold,
     )
     await save_advisory_elevation_stats(session, advisory_elevation_stats)
+
+
+async def store_elevation_tpi_stats(session: AsyncSession, run_parameters_id: int, fire_zone_tpi_stats: FireZoneTPIStats):
+    """
+    Writes elevation TPI statistics to the database.
+
+    :param shape_id: The advisory shape id.
+    :param run_parameters_id: The RunParameter object id associated with this run_type, for_date and run_datetime
+    :param fire_zone_stats: Dictionary keying shape id to a dictionary of classified tpi hfi pixel counts
+    """
+    advisory_tpi_stats_list = []
+    for shape_id, tpi_freq_count in fire_zone_tpi_stats.fire_zone_stats:
+        advisory_tpi_stats = AdvisoryTPIStats(
+            advisory_shape_id=shape_id,
+            run_parameters=run_parameters_id,
+            valley_bottom=tpi_freq_count[1],
+            mid_slope=tpi_freq_count[2],
+            upper_slope=tpi_freq_count[3],
+            pixel_size_metres=fire_zone_tpi_stats.pixel_size_metres,
+        )
+        advisory_tpi_stats_list.append(advisory_tpi_stats)
+
+    await save_advisory_elevation_tpi_stats(session, advisory_tpi_stats_list)

api/app/auto_spatial_advisory/process_elevation_hfi.py

@@ -1,31 +1,27 @@
-""" Code relating to processing HFI data related to elevation
-"""
+"""Code relating to processing HFI data related to elevation"""
+
 import logging
 from datetime import date, datetime
 from time import perf_counter
-from app.auto_spatial_advisory.common import get_s3_key
-from app.auto_spatial_advisory.elevation import process_elevation
+from app.auto_spatial_advisory.elevation import process_elevation_tpi
 from app.auto_spatial_advisory.run_type import RunType
 
 logger = logging.getLogger(__name__)
 
 
 async def process_hfi_elevation(run_type: RunType, run_date: date, run_datetime: datetime, for_date: date):
-    """ Create a new elevation based hfi analysis records for the given date.
+    """Create a new elevation based hfi analysis records for the given date.
 
     :param run_type: The type of run to process. (is it a forecast or actual run?)
     :param run_date: The date of the run to process. (when was the hfi file created?)
     :param for_date: The date of the hfi to process. (when is the hfi for?)
     """
 
-    logger.info('Processing HFI elevation %s for run date: %s, for date: %s', run_type, run_date, for_date)
+    logger.info("Processing HFI elevation %s for run date: %s, for date: %s", run_type, run_date, for_date)
     perf_start = perf_counter()
 
-    key = get_s3_key(run_type, run_date, for_date)
-    logger.info(f'Key to HFI in object storage: {key}')
-
-    await process_elevation(key, run_type, run_datetime, for_date)
+    await process_elevation_tpi(run_type, run_datetime, for_date)
 
     perf_end = perf_counter()
     delta = perf_end - perf_start
-    logger.info('%f delta count before and after processing HFI elevation', delta)
+    logger.info("%f delta count before and after processing HFI elevation", delta)