Included external datasets management

.gitlab-ci.yml

@@ -3,7 +3,7 @@ image: debian:bullseye-slim
 # Make pip cache the installed dependencies
 variables:
   PIP_CACHE_DIR: "$CI_PROJECT_DIR/.cache/pip"
-  OBMGAPANALYSIS_BASE_PATH: "$CI_PROJECT_DIR/tests/data"
+  TEST_BASE_PATH: "$CI_PROJECT_DIR/tests/data"
 cache:
   paths:
     - .cache/pip

obmgapanalysis/database.py

@@ -127,9 +127,7 @@ class Database:
                     with all features and attributes that intersect the given quadtile.
         """
         geometry = tile.geometry
-        print(str(geometry))
-        print(tile.crs)
-        print("epsg:{}".format(crs_number))
+
         if tile.crs != "epsg:{}".format(crs_number):
             project = pyproj.Transformer.from_proj(
                 pyproj.Proj(tile.crs),
@@ -137,7 +135,6 @@ class Database:
                 always_xy=True,
             )
             geometry = transform(project.transform, geometry)
-            print(str(geometry))
 
         sql_query = (
             "SELECT * "

obmgapanalysis/obmgapanalysis.py

@@ -45,8 +45,6 @@ def main():
 
     logger.info("obmgapanalysis has started")
 
-    # If valid instantiate the GHSL dataset as a DataSource class
-
     datasource = DataSource(**datasource_config)
 
     logger.info(
@@ -117,7 +115,22 @@ def main():
                 len(roads_in_tile.index), tile.quadkey
             )
         )
+        # Process the roads query
+        roads_processed = TileProcessor.process_dataframe_with_tile(
+            roads_in_tile, tile=tile, buffer_magnitude=3.0
+        )
+
+        # Substract result from built area
+        refined_built_area = TileProcessor.polygon_difference(
+            clip_built_geometry, roads_processed
+        )
 
+        area = TileProcessor.albers_area_calculation(refined_built_area, tile.crs)
+        logger.info(
+            "Within tile {}, there are {} squared meters of built area.".format(
+                tile.quadkey, area
+            )
+        )
     roads_database.connection.close()
 
     # Leave the program

obmgapanalysis/processor.py

@@ -22,6 +22,9 @@ import rasterio
 from rasterio.mask import mask
 from shapely.geometry import mapping
 from rasterio import features
+import pyproj
+from shapely.ops import transform
+from functools import partial
 import geopandas
 import os
 
@@ -116,7 +119,7 @@ class TileProcessor:
 
         Returns:
             built_area_geometry (shapely.geometry.multipolygon.MultiPolygon): Shapely
-                            polygon or multipolygon
+                            polygon or multipolygon.
         """
         built_pixel_values = datasource.built_pixel_values
         results = []
@@ -142,19 +145,130 @@ class TileProcessor:
     @staticmethod
     def clip_to_tile_extent(polygon, tile):
         """Returns a Shapely geometry object of a (multi)polygon based on the
-        intersection of the given geometry (built areas) with the tile.geometry
+        intersection of the given geometry (built areas) with the tile.geometry.
 
         Args:
             polygon (shapely.geometry.multipolygon.MultiPolygon): Shapely polygon
-                    or multipolygon
+                    or multipolygon.
 
-            tile (tile.Tile): Tile object with quadkey, crs and geometry attributes
+            tile (tile.Tile): Tile object with quadkey, crs and geometry attributes.
 
         Returns:
             polygon (shapely.geometry.multipolygon.MultiPolygon): Shapely polygon
-                    or multipolygon clipped to the shape of tile.geometry
+                    or multipolygon clipped to the shape of tile.geometry.
         """
 
         clipped_polygon = polygon.intersection(tile.geometry)
 
         return clipped_polygon
+
+    @staticmethod
+    def reproject_polygon(input_polygon, input_crs, target_crs):
+        """
+        Returns a (multi)polygon object transformed from input_crs to target_crs. if
+        "aea" is specified as target crs, returns the Albers equal area transformed
+        polygon.
+
+        Args:
+            input_polygon (shapely.geometry.polygon.Polygon): Polygon object to reproject.
+
+            input_crs (str): Initial crs of the input_polygon (e.g. "epsg:4326").
+
+            target_crs (str): Target crs for the final polygon object. If "aea" is
+                    specified, process the complete Albers equal area transformation.
+
+        Returns:
+            geometry (shapely.geometry.polygon.Polygon): Polygon object reprojected to
+                    target_crs.
+        """
+        if target_crs == "aea":
+            project = pyproj.Transformer.from_proj(
+                pyproj.Proj(input_crs),
+                pyproj.Proj("epsg:4326"),
+                always_xy=True,
+            )
+            input_polygon = transform(project.transform, input_polygon)
+
+            project = partial(
+                pyproj.transform,
+                pyproj.Proj("epsg:4326"),
+                pyproj.Proj(
+                    proj="aea", lat_1=input_polygon.bounds[1], lat_2=input_polygon.bounds[3]
+                ),
+            )
+            geometry = transform(project, input_polygon)
+
+        else:
+            project = pyproj.Transformer.from_proj(
+                pyproj.Proj(input_crs),
+                pyproj.Proj(target_crs),
+                always_xy=True,
+            )
+            geometry = transform(project.transform, input_polygon)
+
+        return geometry
+
+    @staticmethod
+    def process_dataframe_with_tile(input_dataframe, tile, buffer_magnitude=False):
+        """
+        Returns a (multi)polygon processed with a tile object and, if desired, buffered
+        by a certain magnitude.
+
+        Args:
+            input_dataframe (geopandas.geodataframe.GeoDataFrame): GeoPandas dataframe to be
+                    processed. May contain buildings, roads or other objects.
+
+            tile (tile.Tile): Tile object with quadkey, crs and geometry attributes.
+
+            buffer_magnitude (float): Numeric magnitude for the polygon buffer (units are
+                    equal to the coordinate system units). Default = False
+        Returns:
+            geometry (shapely.geometry.polygon.Polygon): Polygon object resulting from the
+                    geometric operations.
+        """
+        geometry = input_dataframe.unary_union
+        geometry = TileProcessor.reproject_polygon(geometry, input_dataframe.crs, tile.crs)
+        if buffer_magnitude:
+            geometry = geometry.buffer(buffer_magnitude)
+        geometry = TileProcessor.clip_to_tile_extent(geometry, tile)
+
+        return geometry
+
+    @staticmethod
+    def polygon_difference(reference_polygon, secondary_polygon):
+        """
+        Returns a (multi)polygon object based on the geometric difference between two
+        polygons.
+
+        Args:
+            reference_polygon (shapely.geometry.polygon.Polygon): Geometry of
+                    reference polygon.
+            secondary_polygon (shapely.geometry.polygon.Polygon): Geometry to
+                    compare to the reference geometry
+        Returns:
+            output_polygon (shapely.geometry.polygon.Polygon): Polygon object
+                    representing the difference between the two input polygons.
+        """
+
+        output_polygon = reference_polygon.difference(secondary_polygon)
+
+        return output_polygon
+
+    @staticmethod
+    def albers_area_calculation(input_polygon, polygon_crs):
+        """Return the area of a polygon in the Albers equal area projection.
+
+        Args:
+            input_polygon (shapely.geometry.polygon.Polygon): Input polygon object for
+                    area calculation.
+
+            polygon_crs (str): Initial crs associated to input_polygon.
+
+        Returns:
+            float: Area measured in squared meters for the input_polygon projected to the
+                    Albers equal area system.
+        """
+
+        polygon = TileProcessor.reproject_polygon(input_polygon, polygon_crs, "aea")
+
+        return polygon.area

tests/test_datasource.py

@@ -22,11 +22,15 @@ import geopandas
 import os
 
 
+if "TEST_BASE_PATH" not in os.environ:
+    os.environ["TEST_BASE_PATH"] = "./tests/data"
+
+
 def test_init():
 
     datasource = DataSource(
         crs="epsg:3857",
-        pathname=os.environ["OBMGAPANALYSIS_BASE_PATH"],
+        pathname=os.environ["TEST_BASE_PATH"],
         raster_files_index="GHS_TEST_INDEX.shp",
         built_pixel_values=[6, 5, 4, 3],  # Built pixels in GHSL
     )

tests/test_processor.py

@@ -23,16 +23,24 @@ from obmgapanalysis.datasource import DataSource
 import os
 import numpy
 import affine
+import geopandas
 
 
+if "TEST_BASE_PATH" not in os.environ:
+    os.environ["TEST_BASE_PATH"] = "./tests/data"
+
 # Create an instance for a datasource and a tile
 datasource = DataSource(
     crs="epsg:3857",
-    pathname=os.environ["OBMGAPANALYSIS_BASE_PATH"],
+    pathname=os.environ["TEST_BASE_PATH"],
     raster_files_index="GHS_TEST_INDEX.shp",
     built_pixel_values=[6, 5, 4, 3],
 )
-tile = Tile("122100200320321022", "epsg:3857")
+
+tile = Tile("122100200320321022", datasource.crs)
+roads_query_file = geopandas.read_file(
+    os.path.join(os.environ["TEST_BASE_PATH"], "roads_query.gpkg")
+)
 
 
 def test_intersect_datasource_with_tile():
@@ -81,6 +89,7 @@ def test_poligonize_array():
         "2550000 4629690, 2550000 4629720, 2549970 4629720, "
         "2549970 4629780)))"
     )
+
     built_geometry = TileProcessor.polygonize_array(
         pixel_values, pixel_georeferences, datasource
     )
@@ -108,7 +117,92 @@ def test_clip_to_tile_extent():
     built_geometry = TileProcessor.polygonize_array(
         pixel_values, pixel_georeferences, datasource
     )
-
     clipped_built_geometry = TileProcessor.clip_to_tile_extent(built_geometry, tile)
 
     assert str(clipped_built_geometry) == expected_clipped_geometry
+
+
+def test_process_dataframe_with_tile():
+
+    expected_geometry_roads = (
+        "POLYGON ((2550002.438602296 4629670.059528879, 2550002.454491734 "
+        "4629670.34505833, 2550007.486132718 4629718.894118963, 2550007.534455809 "
+        "4629719.202675823, 2550007.614596096 4629719.504536704, 2550020.082379065 "
+        "4629758.216188851, 2550020.178233149 4629758.472967537, 2550020.297122533 "
+        "4629758.719926265, 2550020.438054846 4629758.955003661, 2550020.599853719 "
+        "4629759.176237528, 2550020.781168612 4629759.381781219, 2550030.443700413 "
+        "4629769.375572084, 2550030.760692946 4629769.659947647, 2550049.072185669 "
+        "4629783.876796038, 2550053.931105311 4629790.803233128, 2550061.2602155 "
+        "4629790.803233128, 2550053.724702562 4629780.061286326, 2550053.541161971 "
+        "4629779.825551415, 2550053.335084511 4629779.609241042, 2550053.108513971 "
+        "4629779.41450048, 2550034.608849079 4629765.051557215, 2550025.57961757 "
+        "4629755.712777978, 2550013.42183824 4629717.963666889, 2550008.437281965 "
+        "4629669.86891498, 2550008.184577072 4629637.929176555, 2550002.184389279 "
+        "4629637.929176555, 2550002.438602296 4629670.059528879))"
+    )
+
+    roads_process = TileProcessor.process_dataframe_with_tile(
+        input_dataframe=roads_query_file, tile=tile, buffer_magnitude=3.0
+    )
+
+    assert str(roads_process) == expected_geometry_roads
+
+
+def test_polygon_difference():
+
+    query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
+    pixel_values, pixel_georeferences = TileProcessor.get_raster_pixels_in_tile(query, tile)
+    built_geometry = TileProcessor.polygonize_array(
+        pixel_values, pixel_georeferences, datasource
+    )
+    clipped_built_geometry = TileProcessor.clip_to_tile_extent(built_geometry, tile)
+
+    roads_process = TileProcessor.process_dataframe_with_tile(
+        input_dataframe=roads_query_file, tile=tile, buffer_magnitude=3.0
+    )
+
+    expected_difference = (
+        "MULTIPOLYGON (((2549970 4629660, 2549940 4629660, 2549939.26359348 4629660, "
+        "2549939.26359348 4629720, 2549940 4629720, 2549970 4629720, 2549970 4629660)), "
+        "((2549970 4629780, 2550000 4629780, 2550044.078820017 4629780, 2550030.760692946 "
+        "4629769.659947647, 2550030.443700413 4629769.375572084, 2550020.781168612 "
+        "4629759.381781219, 2550020.599853719 4629759.176237528, 2550020.438054846 "
+        "4629758.955003661, 2550020.297122533 4629758.719926265, 2550020.178233149 "
+        "4629758.472967537, 2550020.082379065 4629758.216188851, 2550007.614596096 "
+        "4629719.504536704, 2550007.534455809 4629719.202675823, 2550007.486132718 "
+        "4629718.894118963, 2550004.491536504 4629690, 2550000 4629690, 2550000 4629720, "
+        "2549970 4629720, 2549970 4629780)), ((2550013.42183824 4629717.963666889, "
+        "2550025.57961757 4629755.712777978, 2550034.608849079 4629765.051557215, "
+        "2550053.108513971 4629779.41450048, 2550053.335084511 4629779.609241042, "
+        "2550053.541161971 4629779.825551415, 2550053.676985708 4629780, 2550090 "
+        "4629780, 2550090 4629790.803233128, 2550092.13765005 4629790.803233128, "
+        "2550092.13765005 4629750, 2550090 4629750, 2550090 4629720, 2550092.13765005 "
+        "4629720, 2550092.13765005 4629660, 2550090 4629660, 2550090 4629637.929176555, "
+        "2550060 4629637.929176555, 2550060 4629660, 2550030 4629660, 2550030 4629690, "
+        "2550010.523674392 4629690, 2550013.42183824 4629717.963666889)))"
+    )
+
+    output_polygon = TileProcessor.polygon_difference(clipped_built_geometry, roads_process)
+
+    assert str(output_polygon) == expected_difference
+
+
+def test_albers_area_calculation():
+
+    query = TileProcessor.intersect_datasource_with_tile(datasource=datasource, tile=tile)
+    pixel_values, pixel_georeferences = TileProcessor.get_raster_pixels_in_tile(query, tile)
+    built_geometry = TileProcessor.polygonize_array(
+        pixel_values, pixel_georeferences, datasource
+    )
+    clipped_built_geometry = TileProcessor.clip_to_tile_extent(built_geometry, tile)
+    roads_process = TileProcessor.process_dataframe_with_tile(
+        input_dataframe=roads_query_file, tile=tile, buffer_magnitude=3.0
+    )
+
+    polygon_difference = TileProcessor.polygon_difference(clipped_built_geometry, roads_process)
+
+    expected_built_area = 9919.946984796066
+
+    final_built_area = TileProcessor.albers_area_calculation(polygon_difference, tile.crs)
+
+    assert final_built_area == expected_built_area