476 analysis comparison; COUNTRY=myanmar (#880)

* 476 add analysis menus

* add config in json

* Update index.test.tsx.snap

* Update index.tsx

* Update index.tsx

* 476 finalyze ui

* 476 fix test

* 476 namming and format intersect value

* 476 fix build

* 476 add percent in table

* 476 throw error if calculation from client side for intersect_percentage

* Fix intersect_percentage API (#893)

* 476 change % for legend

* 476 add stats_intersect_area

* change units

* 476 fix sort

* Update data-utils.ts

* 476 fix units

* 476 fix intersect area unit in table

* Update zonal_stats.py

* Update zonal_stats.py

* 476 fix legend label

* 476 add % in threshold

* 476 fix tests

* 476 fix useCallback

* Update zonal_stats.py

* update default rounding value to 2

* Update test_calculate.py

* add a warning

* temp revert filter in the backend

* add 'values' for SPI layer config

* Revert "temp revert filter in the backend"

This reverts commit 4069a2d.

* Add values to SPI legends

---------

Co-authored-by: Eric Boucher <[email protected]>
Co-authored-by: Amit Wadhwa <[email protected]>

api/app/main.py

@@ -18,7 +18,7 @@
 from app.models import AcledRequest, RasterGeotiffModel
 from app.timer import timed
 from app.validation import validate_intersect_parameter
-from app.zonal_stats import GroupBy, calculate_stats, get_wfs_response
+from app.zonal_stats import DEFAULT_STATS, GroupBy, calculate_stats, get_wfs_response
 from fastapi import Depends, FastAPI, HTTPException, Path, Query, Response
 from fastapi.middleware.cors import CORSMiddleware
 from fastapi.responses import JSONResponse
@@ -146,7 +146,7 @@ def stats(stats_model: StatsModel) -> list[dict[str, Any]]:
     features = _calculate_stats(
         zones,
         geotiff,
-        stats=" ".join(["min", "max", "mean", "median", "sum", "std"]),
+        stats=" ".join(DEFAULT_STATS),
         prefix="stats_",
         group_by=group_by,
         geojson_out=geojson_out,
@@ -356,7 +356,7 @@ def stats_demo(
     features = _calculate_stats(
         zones_filepath,
         geotiff,
-        stats=" ".join(["min", "max", "mean", "median", "sum", "std"]),
+        stats=" ".join(DEFAULT_STATS),
         prefix="stats_",
         group_by=group_by,
         geojson_out=geojson_out,

api/app/raster_utils.py

@@ -5,7 +5,7 @@
 
 import rasterio
 from app.timer import timed
-from rasterio.warp import Resampling, calculate_default_transform, reproject
+from rasterio.warp import CRS, Resampling, calculate_default_transform, reproject
 
 from .models import FilePath
 
@@ -99,3 +99,17 @@ def reproj_match(
                     dst_crs=dst_crs,
                     resampling=resampling_mode,
                 )
+
+
+def calculate_pixel_area(geotiff_file):
+    with rasterio.open(geotiff_file) as dataset:
+        crs = dataset.crs
+        # Get pixel width and height in the CRS units
+        _, pixel_width, pixel_height = calculate_default_transform(
+            crs, CRS.from_epsg(4326), dataset.width, dataset.height, *dataset.bounds
+        )
+
+        # Convert pixel width and height to square kilometers
+        area = pixel_width * pixel_height / 1e6
+
+        return area

api/app/tests/test_calculate.py

@@ -64,7 +64,6 @@ def test_calculate_stats_with_group_by():
         geotiff,
         group_by="ADM1_PCODE",
         geojson_out=False,
-        intersect_comparison=(operator.lt, 1),
     )
     assert len(features) == 4
 

api/app/zonal_stats.py

@@ -4,26 +4,25 @@
 from datetime import datetime
 from json import dump, load
 from pathlib import Path
-from typing import Any, Optional
+from typing import Any, NewType, Optional
 from urllib.parse import urlencode
 
 import numpy as np
 import rasterio  # type: ignore
 from app.caching import CACHE_DIRECTORY, cache_file, get_json_file, is_file_valid
 from app.models import (
     FilePath,
-    FilterProperty,
     GeoJSON,
     GeoJSONFeature,
-    Geometry,
     GroupBy,
     WfsParamsModel,
     WfsResponse,
 )
-from app.raster_utils import gdal_calc, reproj_match
+from app.raster_utils import calculate_pixel_area, gdal_calc, reproj_match
 from app.timer import timed
 from app.validation import VALID_OPERATORS
 from fastapi import HTTPException
+from osgeo import gdal
 from rasterio.warp import Resampling
 from rasterstats import zonal_stats  # type: ignore
 from shapely.geometry import mapping, shape  # type: ignore
@@ -32,7 +31,10 @@
 logger = logging.getLogger(__name__)
 
 
-DEFAULT_STATS = ["min", "max", "mean", "median"]
+DEFAULT_STATS = ["min", "max", "mean", "median", "sum", "std", "nodata", "count"]
+
+AreaInSqKm = NewType("AreaInSqKm", float)
+Percentage = NewType("Percentage", float)
 
 
 def get_wfs_response(wfs_params: WfsParamsModel) -> WfsResponse:
@@ -269,6 +271,8 @@ def calculate_stats(
     prefix: Optional[str] = "stats_",
     geojson_out: bool = False,
     wfs_response: Optional[WfsResponse] = None,
+    # WARNING - currently, when intersect_comparison is used,
+    # regions with 0 overlap are excluded from the results.
     intersect_comparison: Optional[tuple] = None,
     mask_geotiff: Optional[str] = None,
     mask_calc_expr: Optional[str] = None,
@@ -335,19 +339,24 @@ def calculate_stats(
     # Add function to calculate overlap percentage.
     add_stats = None
     if intersect_comparison is not None:
+        pixel_area = calculate_pixel_area(geotiff)
 
-        def intersect_percentage(masked) -> float:
-            # Get total number of elements in the boundary.
+        def intersect_pixels(masked) -> float:
+            # Get total number of elements matching our operator in the boundary.
             intersect_operator, intersect_baseline = intersect_comparison  # type: ignore
+            # If total is 0, no need to count. This avoids returning NaN.
             total = masked.count()
-            # Avoid dividing by 0
             if total == 0:
-                return 0
-            percentage = (intersect_operator(masked, intersect_baseline)).sum()
-            return percentage / total
+                return 0.0
+            return float((intersect_operator(masked, intersect_baseline)).sum())
+
+        def intersect_area(masked) -> AreaInSqKm:
+            # Area in sq km per pixel value
+            return AreaInSqKm(intersect_pixels(masked) * pixel_area)
 
         add_stats = {
-            "intersect_percentage": intersect_percentage,
+            "intersect_pixels": intersect_pixels,
+            "intersect_area": intersect_area,
         }
 
     try:
@@ -359,15 +368,58 @@ def intersect_percentage(masked) -> float:
             geojson_out=geojson_out,
             add_stats=add_stats,
         )
+
     except rasterio.errors.RasterioError as error:
         logger.error(error)
         raise HTTPException(
             status_code=500, detail="An error occured calculating statistics."
         ) from error
     # This Exception is raised by rasterstats library when feature collection as 0 elements within the feature array.
-    except ValueError as error:
+    except ValueError:
         stats_results = []
 
+    # cleanup data and remove nan values
+    # add intersect stats if requested
+    clean_results = []
+    for result in stats_results:
+        stats_properties: dict = result if not geojson_out else result["properties"]
+
+        # clean results
+        clean_stats_properties = {
+            k: 0 if str(v).lower() == "nan" else v for k, v in stats_properties.items()
+        }
+
+        # calculate intersect_percentage
+        if intersect_comparison is not None:
+            safe_prefix = prefix or ""
+            total = (
+                float(clean_stats_properties[f"{safe_prefix}count"])
+                + clean_stats_properties[f"{safe_prefix}nodata"]
+            )
+            if total == 0:
+                intersect_percentage = Percentage(0.0)
+            else:
+                intersect_percentage = Percentage(
+                    clean_stats_properties[f"{safe_prefix}intersect_pixels"] / total
+                )
+
+            clean_stats_properties = {
+                **clean_stats_properties,
+                f"{safe_prefix}intersect_percentage": intersect_percentage,
+            }
+
+        # merge properties back at the proper level
+        # and filter out properties that have "no" intersection
+        # by setting a limit at 0.005 (0.5%).
+        if intersect_comparison is not None and intersect_percentage < 0.005:
+            continue
+        if not geojson_out:
+            clean_results.append(clean_stats_properties)
+        else:
+            clean_results.append({**result, "properties": clean_stats_properties})
+
+    stats_results = clean_results
+
     if not geojson_out:
         feature_properties = _extract_features_properties(zones_filepath)
 

frontend/src/components/MapView/Layers/AnalysisLayer/index.tsx

@@ -9,7 +9,7 @@ import {
   isAnalysisLayerActiveSelector,
 } from 'context/analysisResultStateSlice';
 import { legendToStops } from 'components/MapView/Layers/layer-utils';
-import { LegendDefinition } from 'config/types';
+import { AggregationOperations, LegendDefinition, units } from 'config/types';
 import {
   BaselineLayerResult,
   ExposedPopulationResult,
@@ -18,6 +18,7 @@ import {
 import { getRoundedData } from 'utils/data-utils';
 import { useSafeTranslation } from 'i18n';
 import { LayerDefinitions } from 'config/utils';
+import { formatIntersectPercentageAttribute } from 'components/MapView/utils';
 
 function AnalysisLayer({ before }: { before?: string }) {
   // TODO maybe in the future we can try add this to LayerType so we don't need exclusive code in Legends and MapView to make this display correctly
@@ -102,18 +103,37 @@ function AnalysisLayer({ before }: { before?: string }) {
           const statisticKey = analysisData.statistic;
           const precision =
             analysisData instanceof ExposedPopulationResult ? 0 : undefined;
+          const formattedProperties = formatIntersectPercentageAttribute(
+            evt.features[0].properties,
+          );
           dispatch(
             addPopupData({
               [analysisData.getStatTitle(t)]: {
-                data: getRoundedData(
-                  get(evt.features[0], ['properties', statisticKey]),
+                data: `${getRoundedData(
+                  formattedProperties[statisticKey],
                   t,
                   precision,
-                ),
+                )} ${units[statisticKey] || ''}`,
                 coordinates,
               },
             }),
           );
+          if (statisticKey === AggregationOperations['Area exposed']) {
+            dispatch(
+              addPopupData({
+                [`${
+                  analysisData.getHazardLayer().title
+                } (Area exposed in km²)`]: {
+                  data: `${getRoundedData(
+                    formattedProperties.stats_intersect_area || null,
+                    t,
+                    precision,
+                  )} ${units.stats_intersect_area}`,
+                  coordinates,
+                },
+              }),
+            );
+          }
         }
 
         if (analysisData instanceof BaselineLayerResult) {