Skip to content
Snippets Groups Projects
Commit ca4585ce authored by Danijel Schorlemmer's avatar Danijel Schorlemmer
Browse files

Make the preprocessor able to handle quadtree-tile geometries

parent b3279ca4
No related branches found
No related tags found
1 merge request!65Make the preprocessor being able to handle quadtree-tile geometries
Pipeline #35353 passed
SOURCES=losscalculator tests utils/preprocessor.py setup.py
SOURCES=losscalculator tests utils/preprocessor.py setup.py utils/database.py
LENGTH=96
check: $(SOURCES)
......
......@@ -16,6 +16,7 @@ setup(
"scipy",
"shapely",
"pyproj",
"pygeotile",
],
extras_require={
"tests": tests_require,
......
......@@ -22,6 +22,8 @@ import sqlite3
import glob
import numpy as np
import shapely.geometry as shgeometry
from shapely.geometry import box
from pygeotile.tile import Tile
# Initialize log
logger = logging.getLogger(__name__)
......@@ -305,7 +307,7 @@ class ExposureDatabase(Database):
sql_statement += "VALUES (%d, %s, %d)" % (building_id, geom, tile_id)
self.cursor.execute(sql_statement)
def insert_tile(self, tile_id):
def insert_tile(self, tile_id, use_tiles=True):
"""
Inserts a tile ID and its respective geometry to the Tile table.
The geometry is computed from the tile ID before inserting.
......@@ -313,13 +315,26 @@ class ExposureDatabase(Database):
Args:
tile_id (int):
ID of the tile
use_tiles (boolean):
Use Quadtree tile IDs instead of IDs of 10-arcsecond cells
"""
if use_tiles:
tile = Tile.from_quad_tree(str(tile_id))
tile_polygon = box(
tile.bounds[0].longitude,
tile.bounds[0].latitude,
tile.bounds[1].longitude,
tile.bounds[1].latitude,
)
tile_geometry = tile_polygon.wkt
else:
tile_geometry = get_geometry_wkt_of_cell(tile_id)
sql_statement = "INSERT INTO Tile "
sql_statement += "(id, geom) "
sql_statement += "VALUES (%d, %s)" % (
tile_id,
"GeomFromText('%s', 4326)" % get_geometry_wkt_of_cell(tile_id),
"GeomFromText('%s', 4326)" % tile_geometry,
)
self.cursor.execute(sql_statement)
......@@ -357,7 +372,9 @@ class ExposureDatabase(Database):
sql_statement += "VALUES (%d, '%s', '%s')" % (district_id, admin_id, admin_name)
self.cursor.execute(sql_statement)
def import_exposure(self, import_search_pattern, building_geometries_filepath):
def import_exposure(
self, import_search_pattern, building_geometries_filepath, use_tiles=True
):
"""
Import an exposure model given as OpenQuake-type exposure CSV files to the database.
In case the model contains assets describing buildings, the building polygons
......@@ -366,12 +383,13 @@ class ExposureDatabase(Database):
Args:
import_search_pattern (str):
The pattern of the exposure model files to import (e.g. 'data/*.csv').
building_geometries_filepath (str):
File path of the building geometries (in EPSG:3857) as CSV file. The file
includes the building IDs (labeled `origin-id`), their corresponding geometries
and IDs of the tiles in which the building centroids are located.
See the example test for a sample version of this file.
use_tiles (boolean):
Use Quadtree tile IDs instead of IDs of 10-arcsecond cells
"""
self.create_tables()
......@@ -382,21 +400,25 @@ class ExposureDatabase(Database):
next(reader)
for building in reader:
building_id = int(building[0][4:])
tile_id = int(
building[2][5:] # Removes the first 5 characters (i.e. the part 'cell_').
)
# Removes the first 5 characters (i.e. the part 'cell_').
tile_id = int(building[2][5:])
if tile_id not in tile_ids:
tile_ids.append(tile_id)
self.insert_building(
building_id,
"ST_Transform(CastToMultiPolygon(GeomFromText('%s', 3857)), 4326)"
% building[1],
# If the geometry is given in 3857, it needs to be transformed into 4326
# "ST_Transform(CastToMultiPolygon(GeomFromText('%s', 3857)), 4326)"
# % building[1],
"GeomFromText('%s', 4326)" % building[1],
tile_id,
)
logger.info("Buildings imported")
self.connection.commit()
else:
print("Please provide the input `-g`,`--building-geometries` if there are any building assets in the input CSV format exposure file. ")
print(
"Please provide the input `-g`,`--building-geometries` if there are any "
+ "building assets in the input CSV format exposure file."
)
# Create temporary storage for taxonomy strings and districts
taxonomy_list = []
district_list = []
......@@ -465,7 +487,7 @@ class ExposureDatabase(Database):
logger.debug("Districts stored")
for tile_id in tile_ids:
self.insert_tile(tile_id)
self.insert_tile(tile_id, use_tiles=use_tiles)
self.connection.commit()
logger.debug("Tiles generated and stored")
......@@ -511,10 +533,6 @@ class ExposureDatabase(Database):
)
writer.writerow(fieldnames)
#if exists BuildingAsset:
# Building asset exposure output
sql_statement = "SELECT BuildingAsset.id, "
sql_statement += "X(ST_Centroid(Building.geom)), "
......@@ -596,8 +614,8 @@ class ExposureDatabase(Database):
asset[9],
asset[10],
asset[11],
'-1',
'POINT EMPTY',
"-1",
"POINT EMPTY",
]
writer.writerow(write_asset)
line_number += 1
......
......@@ -148,6 +148,13 @@ def main():
default="mod_spatialite",
help="File path of the spatialite extension (default: mod_spatialite)",
)
parser.add_argument(
"-c",
"--cells",
required=False,
action="store_true",
help="Use 10-arcsecond cells instead of Quadtree tiles",
)
args = parser.parse_args()
......@@ -158,6 +165,7 @@ def main():
building_geometries_filepath = args.building_geometries
overwrite_exposure_database = args.overwrite
spatialite_filepath = args.spatialite_extension
use_tiles = not args.cells
if command == "import":
if os.path.exists(exposure_database):
......@@ -172,7 +180,9 @@ def main():
except sqlite3.OperationalError:
logger.warning("Spatialite extension cannot be loaded. Exiting ...")
exit()
db.import_exposure(import_search_pattern, building_geometries_filepath)
db.import_exposure(
import_search_pattern, building_geometries_filepath, use_tiles=use_tiles
)
elif command == "export":
db = ExposureDatabase(exposure_database, spatialite_filepath)
try:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment