Added the function to insert buildings and associate them with a location

README.md

@@ -6,6 +6,8 @@ Create exposure model for Coquimbo, Chile.
 
 * `python3`
 * `shapely`
+* `geopandas`
+* `pygeotile`
 
 ## Copyright and copyleft
 

exposurecoquimbo/database.py

@@ -20,8 +20,10 @@
 import logging
 import sqlite3
 import csv
+import geopandas
 from collections import namedtuple
 from shapely import wkt
+from pygeotile.tile import Tile
 
 # Initialize log
 logger = logging.getLogger(__name__)
@@ -144,9 +146,9 @@ class CoquimboDatabase(Database):
 
         # Create table Buildings
         sql_statement = "CREATE TABLE Buildings ("
-        sql_statement += "idx       INTEGER PRIMARY KEY AUTOINCREMENT, "
-        sql_statement += "location  INTEGER, "
-        sql_statement += "quadkey   TEXT)"
+        sql_statement += "idx           INTEGER PRIMARY KEY AUTOINCREMENT, "
+        sql_statement += "location_idx  INTEGER, "
+        sql_statement += "quadkey       TEXT)"
         self.connection.execute(sql_statement)
         sql_statement = "SELECT AddGeometryColumn('Buildings', "
         sql_statement += "'geom', 4326, 'MULTIPOLYGON', 'XY')"
@@ -198,6 +200,8 @@ class CoquimboDatabase(Database):
                  number of buildings,
                  structural value of buildings (total value),
                  number of people inside buildings during night (total value))
+        Return:
+            Index of inserted dataset (Buildings.idx)
         """
 
         sql_statement = "INSERT INTO Assets "
@@ -205,6 +209,24 @@ class CoquimboDatabase(Database):
         sql_statement += "VALUES (%s, '%s', %s, %s, %s)" % dataset
         self.cursor.execute(sql_statement)
 
+    def insert_building(self, dataset):
+        """
+        Inserts a dataset in the Buildings table.
+
+        Args:
+            dataset (tuple):
+                (location_index corresponsing to Locations.idx,
+                 quadkey of the level 18 tile of the building centroid,
+                 building_geometry,
+                 centroid_geometry)
+        """
+        sql_statement = "INSERT INTO Buildings "
+        sql_statement += "(location_idx, quadkey, geom, centroid) "
+        sql_statement += "VALUES (%d, %s, %s, %s)" % dataset
+        self.cursor.execute(sql_statement)
+
+        return self.cursor.lastrowid
+
     def read_sara_exposure(self, exposure_filepath):
         """
         Reads a SARA exposure file into the tables Locations and Assets.
@@ -299,3 +321,54 @@ class CoquimboDatabase(Database):
             logger.debug("Voronoi cell assigned to location %d." % index)
 
         self.connection.commit()
+
+    def add_buildings(self, buildings_filepath):
+        """
+        Reads in OpenStreetMap buildings from a geopackage file. For each building,
+        its centroid is computed and the Voronoi cell is identified in which the
+        building is located. Then the building is associated with the location of
+        SARA data to which the Voronoi cell belongs.
+
+        Args:
+            buildings_filepath (str):
+                File path to a geopackage with OpenStreetMap buildings.
+        """
+
+        buildings = geopandas.read_file(buildings_filepath)
+        for index, building in buildings.iterrows():
+            sql_statement = "SELECT location_idx, geom FROM Voronoi WHERE "
+            sql_statement += "within(GeomFromText('%s', 4326), " % building.geometry.centroid
+            sql_statement += "voronoi.geom)"
+            self.cursor.execute(sql_statement)
+            voronoi = self.cursor.fetchall()
+            try:
+                location_index = voronoi[0][0]
+                centroid = "POINT(%f %f)" % (
+                    building.geometry.centroid.x,
+                    building.geometry.centroid.y,
+                )
+            except IndexError:
+                # If centroid lies on a boundary of the Voronoi cells, move the centroid
+                centroid = "POINT(%f %f)" % (
+                    building.geometry.centroid.x + 0.000001,
+                    building.geometry.centroid.y + 0.000001,
+                )
+                sql_statement = "SELECT location_idx, geom FROM Voronoi WHERE "
+                sql_statement += "within(GeomFromText('%s', 4326), " % centroid
+                sql_statement += "voronoi.geom)"
+                self.cursor.execute(sql_statement)
+                voronoi = self.cursor.fetchall()
+                location_index = voronoi[0][0]
+
+            quadkey = Tile.for_latitude_longitude(
+                longitude=float(building.geometry.centroid.x),
+                latitude=float(building.geometry.centroid.y),
+                zoom=18,
+            ).quad_tree
+            building_geometry = "GeomFromText('%s', 4326)" % building.geometry
+            centroid_geometry = "GeomFromText('%s', 4326)" % centroid
+            dataset = (location_index, quadkey, building_geometry, centroid_geometry)
+            building_idx = self.insert_building(dataset)
+            self.cursor.execute(sql_statement)
+            logger.debug("Building %d added to location %d " % (building_idx, location_index))
+        self.connection.commit()

exposurecoquimbo/exposurecoquimbo.py

@@ -38,6 +38,13 @@ def main():
     parser = argparse.ArgumentParser(description=description)
 
     parser.add_argument("command", type=str, help="'create' to create the exposure database")
+    parser.add_argument(
+        "-b",
+        "--buildings-filepath",
+        required=True,
+        type=str,
+        help="File path of the buildings dataset",
+    )
     parser.add_argument(
         "-d",
         "--database",
@@ -47,14 +54,14 @@ def main():
     )
     parser.add_argument(
         "-e",
-        "--exposure_filepath",
+        "--exposure-filepath",
         required=False,
         type=str,
         help="File path of a SARA exposure-model file (any occupancy type)",
     )
     parser.add_argument(
         "-s",
-        "--spatialite_extension",
+        "--spatialite-extension",
         required=False,
         type=str,
         default="mod_spatialite",
@@ -68,6 +75,7 @@ def main():
     database_filepath = args.database
     exposure_filepath = args.exposure_filepath
     spatialite_filepath = args.spatialite_extension
+    buildings_filepath = args.buildings_filepath
 
     db = CoquimboDatabase(database_filepath, spatialite_filepath)
 
@@ -83,6 +91,8 @@ def main():
         db.create_tables()
         db.read_sara_exposure(exposure_filepath)
         db.add_voronoi_cells()
+        db.add_buildings(buildings_filepath)
+        logger.info("Exposure model created.")
     else:
         logger.warning("Command not recognized. Exiting ...")
 

setup.py

@@ -10,7 +10,7 @@ setup(
     version="0.1",
     description="Create exposure model for Coquimbo",
     license="AGPLv3+",
-    install_requires=["shapely"],
+    install_requires=["shapely", "geopandas", "pygeotile"],
     extras_require={
         "tests": tests_require,
         "linters": linters_require,