feat(building-comparison): support multiple ref data (#768)

add Microsoft Building Footprints dataset --------- Co-authored-by: Gigaszi <[email protected]> Co-authored-by: hn437 <[email protected]>
GIScience · Feb 3, 2024 · 8335cfb · 8335cfb
1 parent d072bd4
commit 8335cfb
Show file tree

Hide file tree

Showing 19 changed files with 592 additions and 385 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,5 +1,13 @@
 # Changelog
 
+## Current Main
+
+### New Features
+
+- building-comparison: support comparison with multiple datasets ([#768])
+
+[#768]: https://github.com/GIScience/ohsome-quality-api/pull/768
+
 ## Release 1.2.0
 
 ### New Features

diff --git a/ohsome_quality_api/definitions.py b/ohsome_quality_api/definitions.py
@@ -22,6 +22,7 @@
         "VNL": "Earth Observation Group Nighttime Light Data",
         "EUBUCCO": "European building stock characteristics in a common and open "
         + "database",
+        "Microsoft Buildings": "Microsoft Building Footprints (ODbL)",
     }
 )
 
@@ -106,6 +107,6 @@ def get_project_keys() -> Iterable[str]:
 
 def get_attribution(data_keys: list) -> str:
     """Return attribution text. Individual attributions are separated by semicolons."""
-    assert set(data_keys) <= {"OSM", "GHSL", "VNL", "EUBUCCO"}
+    assert set(data_keys) <= {"OSM", "GHSL", "VNL", "EUBUCCO", "Microsoft Buildings"}
     filtered = dict(filter(lambda d: d[0] in data_keys, ATTRIBUTION_TEXTS.items()))
     return "; ".join([str(v) for v in filtered.values()])
diff --git a/ohsome_quality_api/geodatabase/client.py b/ohsome_quality_api/geodatabase/client.py
@@ -19,7 +19,7 @@
 import asyncpg
 import geojson
 from asyncpg import Record
-from geojson import Feature, FeatureCollection
+from geojson import Feature, FeatureCollection, MultiPolygon
 
 from ohsome_quality_api.config import get_config_value
 
@@ -69,50 +69,44 @@ async def get_shdi(bpoly: Feature | FeatureCollection) -> list[Record]:
         return await conn.fetch(query, geom)
 
 
-async def get_building_area(bpoly: Feature) -> list[Record]:
-    """Get area of building footprints for a bounding polygon."""
-    file_path = os.path.join(WORKING_DIR, "select_building_area.sql")
+# TODO: Check calls of the function
+async def get_reference_coverage(table_name: str) -> Feature:
+    """Get reference coverage for a bounding polygon."""
+    file_path = os.path.join(WORKING_DIR, "select_coverage.sql")
     with open(file_path, "r") as file:
         query = file.read()
-    geom = str(bpoly.geometry)
     async with get_connection() as conn:
-        return await conn.fetch(query, geom)
+        result = await conn.fetch(query.format(table_name=table_name))
+    return Feature(geometry=geojson.loads(result[0]["geom"]))
 
 
-async def get_eubucco_coverage(inverse: bool) -> list[Record]:
-    file_path = os.path.join(WORKING_DIR, "select_eubucco_coverage.sql")
-    with open(file_path, "r") as file:
-        query = file.read()
-    if inverse:
-        table_name = "eubucco_v0_1_coverage_inversed"
-    else:
-        table_name = "eubucco_v0_1_coverage_simple"
-    query = query.format(table_name=table_name)
-    async with get_connection() as conn:
-        return await conn.fetch(query)
-
-
-async def get_eubucco_coverage_intersection_area(bpoly: Feature) -> list[Record]:
+async def get_intersection_area(bpoly: Feature, table_name: str) -> float:
     """Get ratio of AOI area to intersection area of AOI and coverage geometry.
 
     The result is the ratio of area within coverage (between 0-1) or an empty list if
     AOI lies outside of coverage geometry.
     """
-    file_path = os.path.join(WORKING_DIR, "select_check_eubucco_coverage.sql")
+    file_path = os.path.join(WORKING_DIR, "select_intersection.sql")
     with open(file_path, "r") as file:
         query = file.read()
     geom = str(bpoly.geometry)
     async with get_connection() as conn:
-        return await conn.fetch(query, geom)
+        result = await conn.fetch(query.format(table_name=table_name), geom)
+        if result:
+            return result[0]["area_ratio"]
+        else:
+            return 0.0
 
 
-async def get_eubucco_coverage_intersection(bpoly: Feature) -> Feature:
+async def get_intersection_geom(bpoly: Feature, table_name: str) -> Feature:
     """Get intersection geometry of AoI and coverage geometry."""
-    file_path = os.path.join(WORKING_DIR, "get_coverage_intersection.sql")
+    file_path = os.path.join(WORKING_DIR, "select_intersection.sql")
     with open(file_path, "r") as file:
         query = file.read()
     geom = str(bpoly.geometry)
     async with get_connection() as conn:
-        result = await conn.fetch(query, geom)
-        bpoly["geometry"] = geojson.loads(result[0]["geom"])
-        return bpoly
+        result = await conn.fetch(query.format(table_name=table_name), geom)
+        if result:
+            return Feature(geometry=geojson.loads(result[0]["geom"]))
+        else:
+            return Feature(geometry=MultiPolygon(coordinates=[]))
diff --git a/ohsome_quality_api/geodatabase/get_coverage_intersection.sql b/ohsome_quality_api/geodatabase/get_coverage_intersection.sql
diff --git a/ohsome_quality_api/geodatabase/select_building_area.sql b/ohsome_quality_api/geodatabase/select_building_area.sql
@@ -3,9 +3,8 @@ WITH bpoly AS (
         ST_Setsrid (ST_GeomFromGeoJSON (%s), 4326) AS geom
 )
 SELECT
-    SUM(eubucco.area) as area
-FROM
-    eubucco,
+    SUM({table_name}.area) as area
+FROM {table_name},
     bpoly
 WHERE
-    ST_Intersects (eubucco.centroid, bpoly.geom);
+    ST_Intersects ({table_name}.centroid, bpoly.geom);
diff --git a/...i/geodatabase/select_eubucco_coverage.sql → ...ality_api/geodatabase/select_coverage.sql b/...i/geodatabase/select_eubucco_coverage.sql → ...ality_api/geodatabase/select_coverage.sql
diff --git a/...atabase/select_check_eubucco_coverage.sql → ...y_api/geodatabase/select_intersection.sql b/...atabase/select_check_eubucco_coverage.sql → ...y_api/geodatabase/select_intersection.sql
@@ -4,9 +4,10 @@ WITH bpoly AS (
 )
 SELECT
     -- ratio of area within coverage (empty if outside, between 0-1 if intersection)
-    ST_Area (ST_Intersection (bpoly.geom, coverage.geom)) / ST_Area (bpoly.geom) as area_ratio
+    ST_Area (ST_Intersection (bpoly.geom, coverage.geom)) / ST_Area (bpoly.geom) as area_ratio,
+    ST_AsGeoJSON (ST_Intersection (bpoly.geom, coverage.geom)) AS geom
 FROM
     bpoly,
-    eubucco_v0_1_coverage_simple coverage
+    {table_name} coverage
 WHERE
     ST_Intersects (bpoly.geom, coverage.geom)
diff --git a/ohsome_quality_api/indicators/base.py b/ohsome_quality_api/indicators/base.py
@@ -2,7 +2,7 @@
 from abc import ABCMeta, abstractmethod
 
 import plotly.graph_objects as go
-from geojson import Feature, MultiPolygon, Polygon
+from geojson import Feature, Polygon
 
 from ohsome_quality_api.definitions import get_attribution, get_metadata
 from ohsome_quality_api.indicators.models import IndicatorMetadata, Result
@@ -102,22 +102,26 @@ def attribution(cls) -> str:
         return get_attribution(["OSM"])
 
     @classmethod
-    async def coverage(cls, inverse=False) -> Polygon | MultiPolygon:
-        """Return coverage geometry. Default is global coverage."""
+    async def coverage(cls, inverse=False) -> list[Feature]:
+        """Return coverage geometries. Default is global coverage."""
         if inverse is False:
-            return Polygon(
-                coordinates=[
-                    [
-                        (-180, 90),
-                        (-180, -90),
-                        (180, -90),
-                        (180, 90),
-                        (-180, 90),
-                    ]
-                ]
-            )
+            return [
+                Feature(
+                    geometry=Polygon(
+                        coordinates=[
+                            [
+                                (-180, 90),
+                                (-180, -90),
+                                (180, -90),
+                                (180, 90),
+                                (-180, 90),
+                            ]
+                        ]
+                    )
+                )
+            ]
         else:
-            return Polygon(coordinates=[])
+            return [Feature(Polygon(coordinates=[]))]
 
     @abstractmethod
     async def preprocess(self) -> None:

diff --git a/ohsome_quality_api/indicators/building_comparison/datasets.yaml b/ohsome_quality_api/indicators/building_comparison/datasets.yaml
@@ -0,0 +1,23 @@
+EUBUCCO:
+  name: EUBUCCO
+  link: https://docs.eubucco.com/
+  date: Nov 3, 2022
+  description: >-
+        EUBUCCO is a dataset of building footprints for Europe.
+        It is derived from administrative datasets.
+  color: PURPLE
+  coverage:
+    simple: eubucco_v0_1_coverage_simple
+    inversed: eubucco_v0_1_coverage_inversed
+
+Microsoft Buildings:
+    name: Microsoft Building Footprints
+    link: https://planetarycomputer.microsoft.com/dataset/ms-buildings
+    date: July 5, 2022
+    description: >-
+        Microsoft Building Footprints is a dataset of building footprints for the world.
+        It is derived from satellite imagery.
+    color: ORANGE
+    coverage:
+      simple: microsoft_buildings_coverage_simple
+      inversed: microsoft_buildings_coverage_inversed