TPI proportions

api/alembic/versions/1e6932096921_fuel_area_per_tpi_class.py

@@ -0,0 +1,41 @@
+"""Fuel area per TPI class
+
+Revision ID: 1e6932096921
+Revises: 4014ddf1f874
+Create Date: 2025-01-21 09:48:41.621501
+
+"""
+
+from alembic import op
+import sqlalchemy as sa
+from sqlalchemy.dialects import postgresql
+
+# revision identifiers, used by Alembic.
+revision = "1e6932096921"
+down_revision = "4014ddf1f874"
+branch_labels = None
+depends_on = None
+
+
+def upgrade():
+    op.create_table(
+        "tpi_fuel_area",
+        sa.Column("id", sa.Integer(), nullable=False),
+        sa.Column("advisory_shape_id", sa.Integer(), nullable=False),
+        sa.Column("tpi_class", sa.Enum("valley_bottom", "mid_slope", "upper_slope", name="tpiclassenum"), nullable=False),
+        sa.Column("fuel_area", sa.Float(), nullable=False),
+        sa.ForeignKeyConstraint(
+            ["advisory_shape_id"],
+            ["advisory_shapes.id"],
+        ),
+        sa.PrimaryKeyConstraint("id"),
+        comment="Combustible area in each TPI class per fire zone unit.",
+    )
+    op.create_index(op.f("ix_tpi_fuel_area_advisory_shape_id"), "tpi_fuel_area", ["advisory_shape_id"], unique=False)
+    op.create_index(op.f("ix_tpi_fuel_area_id"), "tpi_fuel_area", ["id"], unique=False)
+
+
+def downgrade():
+    op.drop_index(op.f("ix_tpi_fuel_area_id"), table_name="tpi_fuel_area")
+    op.drop_index(op.f("ix_tpi_fuel_area_advisory_shape_id"), table_name="tpi_fuel_area")
+    op.drop_table("tpi_fuel_area")

api/alembic/versions/fa4b3ecb57fe_populate_tpi_fuel_area_table.py

@@ -0,0 +1,25 @@
+"""Populate tpi_fuel_area table
+
+Revision ID: fa4b3ecb57fe
+Revises: 1e6932096921
+Create Date: 2025-01-21 10:04:48.838953
+
+"""
+
+from app.auto_spatial_advisory.process_tpi_fuel_area import process_tpi_fuel_area
+
+
+# revision identifiers, used by Alembic.
+revision = "fa4b3ecb57fe"
+down_revision = "1e6932096921"
+branch_labels = None
+depends_on = None
+
+
+def upgrade():
+    process_tpi_fuel_area()
+    pass
+
+
+def downgrade():
+    pass

api/app/auto_spatial_advisory/elevation.py

@@ -32,7 +32,6 @@ async def process_elevation_tpi(run_type: RunType, run_datetime: datetime, for_d
     """
     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?)

api/app/auto_spatial_advisory/process_tpi_fuel_area.py

@@ -0,0 +1,118 @@
+import numpy as np
+from osgeo import gdal, ogr, osr
+
+from sqlalchemy import func, select
+from sqlalchemy.orm import Session
+from app import config
+from app.auto_spatial_advisory.process_fuel_type_area import get_fuel_type_s3_key
+from app.db.database import get_write_session_scope
+from app.db.models.auto_spatial_advisory import Shape, TPIClassEnum, TPIFuelArea
+from app.utils.geospatial import GDALResamplingMethod, warp_to_match_raster
+from app.utils.s3 import set_s3_gdal_config
+
+
+REPROJECTED_FUEL_LAYER_NAME = "warped_fuel.tif"
+MASKED_TPI_NAME = "masked_tpi.tif"
+
+
+def store_tpi_area_data(session: Session, advisory_shape_id: int, data: np.ndarray, pixel_size: int):
+    """
+    Save TPIFuelArea records to the API database.
+
+        :param session: A synchronous sqlalchemy session object.
+        :param advisory_shape_id: A fire zone id used as a foreign key to the advisory_shapes table.
+        :param data: A numpy ndarray containing classified TPI values that have been masked by the fuel layer.
+        :param pixel_size: The size of the cells in the TPI layer.
+    """
+    unique_values, counts = np.unique(data, return_counts=True)
+    for value, count in zip(unique_values, counts):
+        if value in TPIClassEnum:
+            tpi_enum = TPIClassEnum(value)
+            fuel_area = count * pixel_size * pixel_size
+            tpi_fuel_area = TPIFuelArea(advisory_shape_id=advisory_shape_id, tpi_class=tpi_enum, fuel_area=fuel_area)
+            session.add(tpi_fuel_area)
+
+
+def prepare_masked_tif():
+    """
+    Creates a masked TPI raster using a classified TPI raster from S3 storage and masking it using the fuel layer
+    also from S3 storage
+    """
+    # Open up our rasters with gdal
+    set_s3_gdal_config()
+    bucket = config.get("OBJECT_STORE_BUCKET")
+    tpi_raster_name = config.get("CLASSIFIED_TPI_DEM_NAME")
+    fuel_raster_key = get_fuel_type_s3_key(bucket)
+    tpi_raster_key = f"/vsis3/{bucket}/dem/tpi/{tpi_raster_name}"
+    fuel_ds: gdal.Dataset = gdal.Open(fuel_raster_key, gdal.GA_ReadOnly)  # LCC projection
+    tpi_ds: gdal.Dataset = gdal.Open(tpi_raster_key, gdal.GA_ReadOnly)  # BC Albers 3005 projection
+
+    # Warp the fuel raster to match extent, spatial reference and cell size of the TPI raster
+    warped_fuel_path = f"/vsimem/{REPROJECTED_FUEL_LAYER_NAME}"
+    warped_fuel_ds: gdal.Dataset = warp_to_match_raster(fuel_ds, tpi_ds, warped_fuel_path, GDALResamplingMethod.NEAREST_NEIGHBOUR)
+
+    # Convert the warped fuel dataset into a binary mask by classifying fuel cells as 1 and non-fuel cells as 0.
+    warped_fuel_band: gdal.Band = warped_fuel_ds.GetRasterBand(1)
+    warped_fuel_data: np.ndarray = warped_fuel_band.ReadAsArray()
+    mask = np.where((warped_fuel_data > 0) & (warped_fuel_data < 99), 1, 0)
+
+    # Some helpful things for creating the final masked TPI raster
+    geo_transform = tpi_ds.GetGeoTransform()
+    tpi_ds_srs = tpi_ds.GetProjection()
+    tpi_band: gdal.Band = tpi_ds.GetRasterBand(1)
+    tpi_data = tpi_band.ReadAsArray()
+
+    # Apply the fuel layer mask to the classified TPI raster and store the result in an in-memory gdal dataset
+    masked_tpi_data = np.multiply(mask, tpi_data)
+    output_driver: gdal.Driver = gdal.GetDriverByName("MEM")
+    masked_tpi_dataset: gdal.Dataset = output_driver.Create(f"/vsimem/{MASKED_TPI_NAME}", xsize=tpi_band.XSize, ysize=tpi_band.YSize, bands=1, eType=gdal.GDT_Byte)
+    masked_tpi_dataset.SetGeoTransform(geo_transform)
+    masked_tpi_dataset.SetProjection(tpi_ds_srs)
+    masked_fuel_type_band: gdal.Band = masked_tpi_dataset.GetRasterBand(1)
+    masked_fuel_type_band.SetNoDataValue(0)
+    masked_fuel_type_band.WriteArray(masked_tpi_data)
+    fuel_ds = None
+    tpi_ds = None
+    return masked_tpi_dataset
+
+
+def prepare_wkt_geom_for_gdal(wkt_geom: str):
+    """
+    Given a wkt geometry as a string, convert it to an ogr.Geometry that can be used by gdal.
+
+    :param wkt_geom: The wky geometry string.
+    :return: An osr.Geometry.
+    """
+    geometry: ogr.Geometry = ogr.CreateGeometryFromWkt(wkt_geom)
+    source_srs = osr.SpatialReference()
+    source_srs.ImportFromEPSG(3005)
+    geometry.AssignSpatialReference(source_srs)
+    transform = osr.CoordinateTransformation(geometry.GetSpatialReference(), source_srs)
+    geometry.Transform(transform)
+    return geometry
+
+
+def process_tpi_fuel_area():
+    """
+    Entry point for calculating the fuel layer masked area of each TPI class (valley bottom, mid slope and upper slope) per fire zone unit.
+    """
+    masked_tpi_ds: gdal.Dataset = prepare_masked_tif()
+    masked_tpi_pixel_size = masked_tpi_ds.GetGeoTransform()[1]
+
+    with get_write_session_scope() as session:
+        # Iterate through the fire zone units from the advisory_shapes table
+        stmt = select(Shape.id, func.ST_AsText(Shape.geom))
+        result = session.execute(stmt)
+        for row in result:
+            # Convert the WKT geometry of a fire zone unit into a form that can be used by gdal
+            geometry = prepare_wkt_geom_for_gdal(row[1])
+
+            # Us gdal.Warp to clip out our fire zone unit from the masked tpi raster and then store areas in the tpi_fuel_area table
+            warp_options = gdal.WarpOptions(cutlineWKT=geometry, cutlineSRS=geometry.GetSpatialReference(), cropToCutline=True)
+            intersected_path = "/vsimem/intersected.tif"
+            intersected_ds: gdal.Dataset = gdal.Warp(intersected_path, masked_tpi_ds, options=warp_options)
+            intersected_band: gdal.Band = intersected_ds.GetRasterBand(1)
+            intersected_data: np.ndarray = intersected_band.ReadAsArray()
+            store_tpi_area_data(session, advisory_shape_id=row[0], data=intersected_data, pixel_size=masked_tpi_pixel_size)
+            intersected_ds = None
+    masked_tpi_ds = None

api/app/db/crud/auto_spatial_advisory.py

@@ -23,6 +23,7 @@
     AdvisoryElevationStats,
     AdvisoryTPIStats,
     ShapeType,
+    TPIFuelArea,
 )
 from app.db.models.hfi_calc import FireCentre
 
@@ -454,6 +455,22 @@ async def get_centre_tpi_stats(session: AsyncSession, fire_centre_name: str, run
     result = await session.execute(stmt)
     return result.all()
 
+async def get_fire_zone_tpi_fuel_areas(session: AsyncSession, fire_zone_id):
+    stmt = select(TPIFuelArea).join(Shape, Shape.id == TPIFuelArea.advisory_shape_id).where(Shape.source_identifier == fire_zone_id)
+    result = await session.execute(stmt)
+    return result.all()
+
+
+async def get_fire_centre_tpi_fuel_areas(session: AsyncSession, fire_centre_name: str):
+    stmt = (
+        select(TPIFuelArea.tpi_class, TPIFuelArea.fuel_area, Shape.source_identifier)
+        .join(Shape, Shape.id == TPIFuelArea.advisory_shape_id)
+        .join(FireCentre, FireCentre.id == Shape.fire_centre)
+        .where(FireCentre.name == fire_centre_name)
+    )
+    result = await session.execute(stmt)
+    return result.all()
+
 
 async def get_provincial_rollup(session: AsyncSession, run_type: RunTypeEnum, run_datetime: datetime, for_date: date) -> List[Row]:
     logger.info("gathering provincial rollup")

api/app/db/models/auto_spatial_advisory.py

@@ -30,6 +30,11 @@ class RunTypeEnum(enum.Enum):
     forecast = "forecast"
     actual = "actual"
 
+class TPIClassEnum(enum.Enum):
+    valley_bottom = 1
+    mid_slope = 2
+    upper_slope = 3
+
 
 class ShapeType(Base):
     """Identify some kind of area type, e.g. "Zone", or "Fire" """
@@ -227,3 +232,15 @@ class CriticalHours(Base):
     fuel_type = Column(Integer, ForeignKey(SFMSFuelType.id), nullable=False, index=True)
     start_hour = Column(Integer, nullable=False)
     end_hour = Column(Integer, nullable=False)
+
+class TPIFuelArea(Base):
+    """
+    Combustible area in each TPI class per fire zone unit.
+    """
+
+    __tablename__ = "tpi_fuel_area"
+    __table_args__ = {"comment": "Combustible area in each TPI class per fire zone unit."}
+    id = Column(Integer, primary_key=True, index=True)
+    advisory_shape_id = Column(Integer, ForeignKey(Shape.id), nullable=False, index=True)
+    tpi_class = Column(Enum(TPIClassEnum), nullable=False)
+    fuel_area = Column(Float, nullable=False)

api/app/routers/fba.py

@@ -10,6 +10,8 @@
 from app.db.crud.auto_spatial_advisory import (
     get_all_sfms_fuel_types,
     get_all_hfi_thresholds,
+    get_fire_centre_tpi_fuel_areas,
+    get_fire_zone_tpi_fuel_areas,
     get_hfi_area,
     get_precomputed_stats_for_shape,
     get_provincial_rollup,
@@ -18,11 +20,12 @@
     get_centre_tpi_stats,
     get_zone_ids_in_centre,
 )
-from app.db.models.auto_spatial_advisory import RunTypeEnum
+from app.db.models.auto_spatial_advisory import RunTypeEnum, TPIClassEnum
 from app.schemas.fba import (
     AdvisoryCriticalHours,
     ClassifiedHfiThresholdFuelTypeArea,
     FireCenterListResponse,
+    FireCentreTPIResponse,
     FireShapeAreaListResponse,
     FireShapeArea,
     FireZoneTPIStats,
@@ -166,32 +169,64 @@ async def get_fire_zone_tpi_stats(fire_zone_id: int, run_type: RunType, run_date
     async with get_async_read_session_scope() as session:
         stats = await get_zonal_tpi_stats(session, fire_zone_id, run_type, run_datetime, for_date)
         square_metres = math.pow(stats.pixel_size_metres, 2) if stats is not None else None
+        tpi_fuel_stats = await get_fire_zone_tpi_fuel_areas(session, fire_zone_id)
+        valley_bottom_tpi = None
+        mid_slope_tpi = None
+        upper_slope_tpi = None
+
+        for tpi_fuel_stat in tpi_fuel_stats:
+            if tpi_fuel_stat.tpi_class == TPIClassEnum.valley_bottom:
+                valley_bottom_tpi = tpi_fuel_stat.fuel_area
+            elif tpi_fuel_stat.tpi_class == TPIClassEnum.mid_slope:
+                mid_slope_tpi = tpi_fuel_stat.fuel_area
+            elif tpi_fuel_stat.tpi_class == TPIClassEnum.upper_slope:
+                upper_slope_tpi = tpi_fuel_stat.fuel_area
         return FireZoneTPIStats(
             fire_zone_id=fire_zone_id,
-            valley_bottom=stats.valley_bottom * square_metres if stats is not None else None,
-            mid_slope=stats.mid_slope * square_metres if stats is not None else None,
-            upper_slope=stats.upper_slope * square_metres if stats is not None else None,
+            valley_bottom_hfi=stats.valley_bottom * square_metres if stats is not None else None,
+            valley_bottom_tpi=valley_bottom_tpi,
+            mid_slope_hfi=stats.mid_slope * square_metres if stats is not None else None,
+            mid_slope_tpi=mid_slope_tpi,
+            upper_slope_hfi=stats.upper_slope * square_metres if stats is not None else None,
+            upper_slope_tpi=upper_slope_tpi,
         )
 
 
-@router.get("/fire-centre-tpi-stats/{run_type}/{for_date}/{run_datetime}/{fire_centre_name}", response_model=dict[str, List[FireZoneTPIStats]])
+@router.get("/fire-centre-tpi-stats/{run_type}/{for_date}/{run_datetime}/{fire_centre_name}", response_model=FireCentreTPIResponse)
 async def get_fire_centre_tpi_stats(fire_centre_name: str, run_type: RunType, run_datetime: datetime, for_date: date, _=Depends(authentication_required)):
     """Return the elevation TPI statistics for each advisory threshold for a fire centre"""
     logger.info("/fba/fire-centre-tpi-stats/")
     async with get_async_read_session_scope() as session:
         tpi_stats_for_centre = await get_centre_tpi_stats(session, fire_centre_name, run_type, run_datetime, for_date)
+        tpi_fuel_stats = await get_fire_centre_tpi_fuel_areas(session, fire_centre_name)
 
-        data = []
+        hfi_tpi_areas_by_zone = []
         for row in tpi_stats_for_centre:
+            fire_zone_id = row.source_identifier
             square_metres = math.pow(row.pixel_size_metres, 2)
-
-            data.append(
+            tpi_fuel_stats_for_zone = [stats for stats in tpi_fuel_stats if stats[2] == fire_zone_id]
+            valley_bottom_tpi = None
+            mid_slope_tpi = None
+            upper_slope_tpi = None
+
+            for tpi_fuel_stat in tpi_fuel_stats_for_zone:
+                if tpi_fuel_stat[0] == TPIClassEnum.valley_bottom:
+                    valley_bottom_tpi = tpi_fuel_stat[1]
+                elif tpi_fuel_stat[0] == TPIClassEnum.mid_slope:
+                    mid_slope_tpi = tpi_fuel_stat[1]
+                elif tpi_fuel_stat[0] == TPIClassEnum.upper_slope:
+                    upper_slope_tpi = tpi_fuel_stat[1]
+
+            hfi_tpi_areas_by_zone.append(
                 FireZoneTPIStats(
-                    fire_zone_id=row.source_identifier,
-                    valley_bottom=row.valley_bottom * square_metres,
-                    mid_slope=row.mid_slope * square_metres,
-                    upper_slope=row.upper_slope * square_metres,
+                    fire_zone_id=fire_zone_id,
+                    valley_bottom_hfi=row.valley_bottom * square_metres,
+                    valley_bottom_tpi=valley_bottom_tpi,
+                    mid_slope_hfi=row.mid_slope * square_metres,
+                    mid_slope_tpi=mid_slope_tpi,
+                    upper_slope_hfi=row.upper_slope * square_metres,
+                    upper_slope_tpi=upper_slope_tpi,
                 )
             )
 
-        return {fire_centre_name: data}
+    return FireCentreTPIResponse(fire_centre_name=fire_centre_name, firezone_tpi_stats=hfi_tpi_areas_by_zone)

api/app/schemas/fba.py

@@ -124,9 +124,17 @@ class FireZoneTPIStats(BaseModel):
     """Classified TPI areas of the fire zone contributing to the HFI >4k. Each area is in square metres."""
 
     fire_zone_id: int
-    valley_bottom: Optional[int]
-    mid_slope: Optional[int]
-    upper_slope: Optional[int]
+    valley_bottom_hfi: Optional[int]
+    valley_bottom_tpi: Optional[float]
+    mid_slope_hfi: Optional[int]
+    mid_slope_tpi: Optional[float]
+    upper_slope_hfi: Optional[int]
+    upper_slope_tpi: Optional[float]
+
+
+class FireCentreTPIResponse(BaseModel):
+    fire_centre_name: str
+    firezone_tpi_stats: List[FireZoneTPIStats]
 
 
 class FireZoneElevationStatsByThreshold(BaseModel):