Commit 8706cf29 authored by Cecilia Nievas's avatar Cecilia Nievas
Browse files

Add enter to print statement in refactoring of SERA_mapping_admin_units_to_cells.py

parent 08395980
......@@ -24,9 +24,9 @@ Section 2.6: Seismic Hazard and Risk Dynamics
SERA_mapping_admin_units_to_cells
=================================
This code determines the mapping between the administrative units used in SERA
and the grid cells we are defining with a 10 arcsec resolution. It does not
carry out the assignment of SERA buildings to each cell but only creates this
mapping, as there can be different ways of doing the assignment of SERA
and the grid cells we are defining with a 10 arcsec resolution. It does not
carry out the assignment of SERA buildings to each cell but only creates this
mapping, as there can be different ways of doing the assignment of SERA
buildings that do not need this mapping to be re-done each time.
"""
......@@ -37,7 +37,7 @@ import geopandas as gpd
import datetime
import logging
import GDE_TOOLS_world_grid as gdet_grid
import GDE_TOOLS_read_SERA as gdet_sera
import GDE_TOOLS_read_SERA as gdet_sera
import GDE_TOOLS_general as gdet_gral
import GDE_TOOLS_psql as gdet_psql
import GDE_TOOLS_read_config_file as gdet_conf
......@@ -49,8 +49,8 @@ def run_this_file(config_dict):
####################################################
print('Processing configuration parameters...')
# Path for output:
out_path = config_dict['File Paths']['out_path']
# SERA models path:
out_path = config_dict['File Paths']['out_path']
# SERA models path:
sera_models_path = config_dict['File Paths']['sera_models_path']
# SERA boundaries shapefiles path:
sera_shp_path = config_dict['File Paths']['sera_boundaries_path']
......@@ -58,14 +58,14 @@ def run_this_file(config_dict):
DB_username_grid = config_dict['Tiles Database']['db_tiles_username']
DB_name_grid = config_dict['Tiles Database']['db_tiles_name']
DB_schema_name_grid = config_dict['Tiles Database']['db_tiles_schema_name']
DB_table_name_grid = config_dict['Tiles Database']['db_tiles_table_name']
DB_table_name_grid = config_dict['Tiles Database']['db_tiles_table_name']
# Countries to process:
countries = config_dict['SERA_mapping_admin_units_to_cells']['countries'].split(', ')
####################################################
# SET UP LOGGER
####################################################
now = datetime.datetime.now()
now_str = '{:s}_{:s}_{:s}_{:s}_{:s}_{:s}'.format(str(now.year),
now_str = '{:s}_{:s}_{:s}_{:s}_{:s}_{:s}'.format(str(now.year),
str(now.month).zfill(2),
str(now.day).zfill(2),
str(now.hour).zfill(2),
......@@ -80,77 +80,77 @@ def run_this_file(config_dict):
####################################################
num_countries = len(countries)
for k,country_name in enumerate(countries):
for result in process_country(country_name, k, num_countries,
for result in process_country(country_name, k, num_countries,
sera_models_path, sera_shp_path, logger):
# Unpack results:
cell_id, adm_level, admin_id, area_intersect, geom_intersect, cases_in_lev_grs, shp_not_found = result
# Write entry to database:
err_msg = gdet_psql.insert_data_in_cell_admin_unit_mapping(DB_name_grid,
DB_username_grid,
DB_schema_name_grid,
DB_table_name_grid,
cell_id,
country_name,
adm_level,
admin_id,
area_intersect,
geom_intersect,
err_msg = gdet_psql.insert_data_in_cell_admin_unit_mapping(DB_name_grid,
DB_username_grid,
DB_schema_name_grid,
DB_table_name_grid,
cell_id,
country_name,
adm_level,
admin_id,
area_intersect,
geom_intersect,
cases_in_lev_grs)
if len(err_msg)>0:
logger.info(country_name+': '+err_msg)
logger.debug('\r %s' % (gdet_gral.add_to_log_shapefiles_found(country_name,
logger.debug('\r %s' % (gdet_gral.add_to_log_shapefiles_found(country_name,
shp_not_found)))
def process_country(country_str, country_pos, number_countries,
def process_country(country_str, country_pos, number_countries,
sera_mods_path, sera_shape_path, logger):
"""
country_str = country name
sera_mods_path = path to the SERA model files
sera_shape_path = path to the shapefiles used by SERA
"""
logger.info('\nWorking on %s, country %i of %i\n'
% (country_str, country_pos+1, number_countries))
logger.info('\nWorking on %s, country %i of %i\n'
% (country_str, country_pos+1, number_countries))
# There is repetition of admin levels across Res/Com/Ind.
# Pre-analyse the case so as to carry out spatial intersections once.
packed_res = gdet_sera.get_admin_level_definition(country_str,
packed_res = gdet_sera.get_admin_level_definition(country_str,
sera_mods_path,
full_files=True,
output_dfs=True)
admin_levels, level_groups, level_groups_levels, dfs_dict = packed_res
# Go by group of any combination of Res/Com/Ind with the same admin level
# definition (given by level_groups_levels):
for j in range(len(level_groups)):
for j in range(len(level_groups)):
try:
adm_boundary = gpd.GeoDataFrame.from_file(os.path.join(sera_shape_path,
'Adm%i_%s.shp'
'Adm%i_%s.shp'
% (level_groups_levels[j], country_str)))
except:
logger.info(' Shapefile not found for Admin Level %i'
% (level_groups_levels[j]))
# Keep track of admin IDs not found in the shapefile:
shp_not_found= 0
# Keep track of admin units that checked, because not all units might
# exist in each case (Res, Com, Ind):
admin_units_investigated= []
admin_units_investigated= []
# Go by case (Res, Com or Ind):
for case in level_groups[j]:
for case in level_groups[j]:
#unique_admins= np.unique(dfs_dict[case]['id_%i' % (level_groups_levels[j])].values)
unique_admins= np.unique(dfs_dict[case]['id_%i' % (level_groups_levels[j])].values[np.where(dfs_dict[case]['id_2']==4)[0]])
# Go one by one the admin units named in the exposure CSV file:
for m, admin_id in enumerate(unique_admins):
for m, admin_id in enumerate(unique_admins):
if admin_id in admin_units_investigated:
continue
admin_units_investigated.append(admin_id)
admin_units_investigated.append(admin_id)
# If the admin id cannot be found in the corresponding shapefile:
if admin_id not in adm_boundary['ID_%i' % (level_groups_levels[j])].values:
if admin_id not in adm_boundary['ID_%i' % (level_groups_levels[j])].values:
shp_not_found= shp_not_found + 1
continue
......@@ -170,22 +170,22 @@ def process_country(country_str, country_pos, number_countries,
cell_id, area_intersect, geom_intersect = results_adm
yield cell_id, level_groups_levels[j], admin_id, area_intersect, geom_intersect, level_groups[j], shp_not_found
def process_admin_unit(adm_boundary_gdf, which_row_of_gdf, logger):
def process_admin_unit(adm_boundary_gdf, which_row_of_gdf, logger):
"""
adm_boundary_gdf = GeoPandas DataFrame with admin boundaries
which_row_of_gdf = row of adm_boundary_gdf to be processed
"""
# Bounding box of the admin unit:
# Bounding box of the admin unit:
bbox = adm_boundary_gdf['geometry'].values[which_row_of_gdf[0]].bounds
# GeoPandas DF with the cells associated with this admin unit:
gdf_cells = gdet_grid.get_cells_in_bbox_as_geopandas(bbox[0], bbox[1],
bbox[2], bbox[3])
gdf_cells = gdet_grid.get_cells_in_bbox_as_geopandas(bbox[0], bbox[1],
bbox[2], bbox[3])
# gdf_intersect has as many rows as cells do intersect this admin unit:
gdf_intersect = gpd.sjoin(gdf_cells,
adm_boundary_gdf.iloc[[which_row_of_gdf[0]],:],
how='right')
gdf_intersect = gpd.sjoin(gdf_cells,
adm_boundary_gdf.iloc[[which_row_of_gdf[0]],:],
how='right')
# Loop through all cells that intersect the admin boundary:
for i in range(gdf_intersect.shape[0]):
print('\r Intersection of cell '+str(gdf_intersect['cell_id'].values[i]), end='')
......@@ -196,19 +196,19 @@ def process_admin_unit(adm_boundary_gdf, which_row_of_gdf, logger):
geom_intersection = geom_adm_unit.intersection(geom_cell)
area_intersection = gdet_grid.get_area_polygon_on_Earth(geom_intersection,
units='km2')
print('\n')
yield gdf_intersect['cell_id'].values[i], area_intersection, geom_intersection
if __name__=='__main__':
# This code needs to be run from the command line as python3 namefile.py configfile.ini
# sys.argv retrieves all the commands entered in the command line;
# sys.argv retrieves all the commands entered in the command line;
# position [0] is this code, position [1] will be the config file name
config_filename = sys.argv[1]
section_names_to_validate = ['File Paths',
'Tiles Database',
config_filename = sys.argv[1]
section_names_to_validate = ['File Paths',
'Tiles Database',
os.path.basename(__file__).split('.')[0]]
config_dict = gdet_conf.read_config_parameters(os.path.join(os.getcwd(),
config_filename),
config_dict = gdet_conf.read_config_parameters(os.path.join(os.getcwd(),
config_filename),
section_names_to_validate)
run_this_file(config_dict)
run_this_file(config_dict)
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