Commit 0ac53958 authored by Simantini Shinde's avatar Simantini Shinde
Browse files

Imported districts for industrial aem

parent f5ba9af8
Pipeline #46007 passed with stage
in 1 minute and 51 seconds
#!/usr/bin/env python3
# Copyright (C) 2022:
# Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero
# General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.
import logging
from exposurelib.database import SpatialiteDatabase
import shapely.wkt
import shapely.wkb
import pyproj
from shapely.ops import transform
from exposurejapan import constants
# Initialize log
logger = logging.getLogger(__name__)
class BaseAggregatedExposure(SpatialiteDatabase):
"""
The BaseAggregatedExposure class represents a Spatialite database for the aggregated
exposure model for Japan. It is derived from the generic Database class.
Args:
database_filepath (str):
Filepath for the Spatialite database file.
spatialite_filepath (str):
Filepath of the Spatialite extension.
Attributes:
database_filepath (str):
Spatialite database filepath.
spatialite_filepath (str):
Filepath to the Spatialite extension.
"""
def __init__(self, database_filepath, spatialite_filepath="mod_spatialite"):
SpatialiteDatabase.__init__(self, database_filepath, spatialite_filepath)
def create_district_table(self):
"""
Creates the `District` table in the database.
District : Stores the districts with their IDs, names and geometries.
"""
sql_statement = "CREATE TABLE District ("
sql_statement += "id INTEGER PRIMARY KEY, "
sql_statement += "name TEXT, "
sql_statement += "area_size REAL)"
self.connection.execute(sql_statement)
sql_statement = "SELECT AddGeometryColumn('District', 'geom', 4326, "
sql_statement += "'MULTIPOLYGON', 'XY')"
self.connection.execute(sql_statement)
logger.debug("Table District created")
@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.
More information on the Albers Equal-Area projection can be found at:
https://pro.arcgis.com/en/pro-app/latest/help/mapping/properties/albers.htm
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)
# Get the bounding box of the polygon and sort latitudes and longitudes
bbox = input_polygon.bounds
minx, maxx = sorted([bbox[0], bbox[2]])
miny, maxy = sorted([bbox[1], bbox[3]])
# Use the resulting coordinates to perform the transformation
project_albers = pyproj.Proj(
"+proj=aea +lat_1={} +lat_2={} +lat_0={} +lon_0={}".format(
miny, maxy, (miny + maxy) / 2.0, (minx + maxx) / 2.0
)
)
geometry = transform(project_albers, 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
def import_districts_and_boundaries(self, district_boundary_filepath):
"""
Imports all districts and boundaries from a prepared GeoPackage file.
('estat_bound_municipal'). The function opens this file as Spatialite database
and copies all necessary district and boundary information into the `District` table.
Args:
district_boundary_filepath (str):
Filepath to the boundary file
"""
# Attach GeoPackage boundary database to this database
sql_statement = "SELECT EnableGpkgAmphibiousMode();"
self.connection.execute(sql_statement)
sql_statement = "ATTACH '%s' AS BoundaryDatabase" % district_boundary_filepath
logger.debug(sql_statement)
self.connection.execute(sql_statement)
# From the boundary database select `id`, `name` and `geometry`
# and insert into the `District` table
sql_statement = (
"INSERT INTO District (id, name, geom) "
"SELECT key_code_ward, CITY_NAME, "
"ST_Transform(CastToMultiPolygon(GeomFromGPB(geom)), "
"%d) " % constants.WGS84
)
sql_statement += "FROM BoundaryDatabase.Boundary"
logger.debug(sql_statement)
self.cursor.execute(sql_statement)
self.connection.commit()
# Calculate the area for each district using Albers Equal-Area projection
# and insert into the `District` table
sql_statement = "SELECT id, AsBinary(geom) FROM District"
self.cursor.execute(sql_statement)
districts = self.cursor.fetchall()
for district in districts:
district_geom = shapely.wkb.loads(district[1])
district_size = (
# Reproject district geometry to Albers Equal-Area projection
# and extract its area
self.reproject_polygon(district_geom, "epsg:%d" % constants.WGS84, "aea").area
* 0.000001 # Convert from square meters to square kilometers
)
sql_statement = "UPDATE District SET area_size = %f WHERE id = %s" % (
district_size,
district[0],
)
logger.debug(sql_statement)
self.cursor.execute(sql_statement)
self.connection.commit()
# Detach GeoPackage boundary database
sql_statement = "DETACH DATABASE 'BoundaryDatabase'"
logger.debug(sql_statement)
self.connection.execute(sql_statement)
logger.info("Districts and boundaries added")
......@@ -20,7 +20,8 @@
import logging
import sys
import sqlite3
from residential import ResidentialAggregatedExposure
from exposurejapan.residential import ResidentialAggregatedExposure
from exposurejapan.industrial import IndustrialAggregatedExposure
# Add a logger printing error, warning, info and debug messages to the screen
......@@ -54,6 +55,18 @@ def main():
res_db.import_gem_building_class_mapping("mapping_files/GEM_Tax_Building_class.xlsx")
res_db.assign_building_assets()
res_db.export_to_csv("output.csv")
# Run the industrial aggregated exposure model
ind_db = IndustrialAggregatedExposure("aggregated_industrial.sqlite")
try:
ind_db.connect()
except sqlite3.OperationalError:
logger.warning("Spatialite extension cannot be loaded. Exiting ...")
exit()
ind_db.create_district_table()
ind_db.import_districts_and_boundaries("data/Boundary.gpkg")
# Leave the program
sys.exit()
......
#!/usr/bin/env python3
# Copyright (C) 2022:
# Helmholtz-Zentrum Potsdam Deutsches GeoForschungsZentrum GFZ
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero
# General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.
import logging
from exposurejapan.baseaggregatedexposure import BaseAggregatedExposure
# Initialize log
logger = logging.getLogger(__name__)
class IndustrialAggregatedExposure(BaseAggregatedExposure):
"""
The IndustrialDatabase class represents a Spatialite database for the industrial
aggregated exposure model for Japan. It is derived from the generic `BaseAggregatedExposure`
class.
Args:
database_filepath (str):
Filepath for the Spatialite database file.
spatialite_filepath (str):
Filepath of the Spatialite extension (optional).
Attributes:
database_filepath (str):
Spatialite database filepath.
spatialite_filepath (str, optional):
Filepath to the Spatialite extension (default: 'mod_spatialite').
"""
def __init__(self, database_filepath, spatialite_filepath="mod_spatialite"):
BaseAggregatedExposure.__init__(self, database_filepath, spatialite_filepath)
......@@ -19,7 +19,7 @@
import logging
from exposurelib.database import SpatialiteDatabase
import pandas
import constants
from exposurejapan import constants
import shapely.wkt
import shapely.wkb
import pyproj
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment