Created SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles.py

Created alternative function to retrieve GHS built-up area from OBM
Tiles database in GDE_TOOLS_psql.py (connection to the database opened
before calling the function). Updated config file template and reader
accordingly.

GDE_TOOLS_psql.py

@@ -623,7 +623,10 @@ def get_ghs_built_up_area_from_obm_tiles(grid_cell_id, db_hostname, db_name,
     """This function retrieves the GHS built-up area from the OBM tiles database for a specified
     zoom level 18 quadkey.
 
-    Note that a password is required to access the database.
+    Note that a password is required to access the database. If not provided as input, this
+    function will prompt the user to indicate it. If the connection to the database and a cursor
+    are opened before calling this function, consider using the function
+    get_ghs_built_up_area_from_obm_tiles_opened_connection() instead.
 
     Args:
         grid_cell_id (str): Quadkey of the zoom level 18 target tile.
@@ -634,9 +637,9 @@ def get_ghs_built_up_area_from_obm_tiles(grid_cell_id, db_hostname, db_name,
         source_id (int): Integer indicating the source ID of the GHS built-up area. Deafault: 1.
     
     Returns:
-        ghs_built_up (float): Built-up area in the target tile according to the Global Human
-                              Settlement GHSL-BUILT multitemporal dataset available in 30-m
-                              resolution at https://ghsl.jrc.ec.europa.eu/download.php.
+        ghs_built_up (float): Built-up area (m2) in the target tile according to the Global
+                              Human Settlement GHSL-BUILT multitemporal dataset available in
+                              30-m resolution at https://ghsl.jrc.ec.europa.eu/download.php.
     """
 
     # Connect to the database:
@@ -659,8 +662,52 @@ def get_ghs_built_up_area_from_obm_tiles(grid_cell_id, db_hostname, db_name,
     elif len(exec_result) == 1:
         ghs_built_up = exec_result[0][0]
     else:  # More than one entries found, this is an error
-        print("ERROR IN get_ghs_built_up_area_from_obm_tiles: MORE THAN ONE ENTRIES FOUND FOR CELL ID %s"
+        print("ERROR IN get_ghs_built_up_area_from_obm_tiles: "
+              "MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
               % (grid_cell_id))
         ghs_built_up = -999.9
 
     return ghs_built_up
+
+
+def get_ghs_built_up_area_from_obm_tiles_opened_connection(grid_cell_id, open_cursor,
+                                                           db_table, source_id=1):
+    """This function retrieves the GHS built-up area from the OBM tiles database for a specified
+    zoom level 18 quadkey.
+
+    Note that this function requires that a connection to the database has been established and
+    a cursor opened. The function get_ghs_built_up_area_from_obm_tiles() can be used if this is
+    not the case.
+
+    Args:
+        grid_cell_id (str): Quadkey of the zoom level 18 target tile.
+        open_cursor: Cursor connecting to the OBM tiles database (already open).
+        db_table (str): Name of the table where this information is stored.
+        source_id (int): Integer indicating the source ID of the GHS built-up area. Deafault: 1.
+    
+    Returns:
+        ghs_built_up (float): Built-up area (m2) in the target tile according to the Global
+                              Human Settlement GHSL-BUILT multitemporal dataset available in
+                              30-m resolution at https://ghsl.jrc.ec.europa.eu/download.php.
+    """
+
+    # Query the database:
+    sql_query = "SELECT built_area_size FROM %s WHERE (quadkey='%s' AND source_id=%s);" \
+                % (db_table, grid_cell_id, source_id)
+    open_cursor.execute(sql_query)
+    exec_result = open_cursor.fetchall()
+
+    if len(exec_result) == 0:  # If the quadkey is not found the built-up area is zero
+        ghs_built_up = 0.0
+    elif len(exec_result) == 1:
+        ghs_built_up = exec_result[0][0]
+    else:  # More than one entries found, this is an error
+        print("ERROR IN get_ghs_built_up_area_from_obm_tiles: "
+              "MORE THAN ONE ENTRY FOUND FOR CELL ID %s"
+              % (grid_cell_id))
+        ghs_built_up = -999.9
+
+    return ghs_built_up
+
+
+

GDE_TOOLS_read_config_file.py

@@ -296,6 +296,27 @@ def check_parameters(config, section_name):
             "SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5", "ghs_input_filename"
         ):
             raise IOError("ERROR!! ghs_input_filename PARAMETER MISSING FROM CONFIG FILE!!")
+    elif section_name == "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles":
+        if not config.has_option(
+            "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles", "obm_tiles_db_hostname"
+        ):
+            raise IOError("ERROR!! obm_tiles_db_hostname PARAMETER MISSING FROM CONFIG FILE!!")
+        if not config.has_option(
+            "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles", "obm_tiles_db_name"
+        ):
+            raise IOError("ERROR!! obm_tiles_db_name PARAMETER MISSING FROM CONFIG FILE!!")
+        if not config.has_option(
+            "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles", "obm_tiles_db_username"
+        ):
+            raise IOError("ERROR!! obm_tiles_db_username PARAMETER MISSING FROM CONFIG FILE!!")
+        if not config.has_option(
+            "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles", "obm_tiles_db_table"
+        ):
+            raise IOError("ERROR!! obm_tiles_db_table PARAMETER MISSING FROM CONFIG FILE!!")
+        if not config.has_option(
+            "SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles", "obm_tiles_source_id"
+        ):
+            raise IOError("ERROR!! obm_tiles_source_id PARAMETER MISSING FROM CONFIG FILE!!")
     elif section_name == "SERA_mapping_admin_units_to_cells_add_Sat":
         if not config.has_option("SERA_mapping_admin_units_to_cells_add_Sat", "apply_model"):
             raise IOError("ERROR!! apply_model PARAMETER MISSING FROM CONFIG FILE!!")

GDE_config_file_TEMPLATE.ini

@@ -202,6 +202,13 @@ occupancy_cases = Res, Com, Ind
 # Countries to process (if many, separate with comma and space):
 countries = Greece
 
+[SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles]
+obm_tiles_db_hostname = nozomi.gfz-potsdam.de
+obm_tiles_db_name = obm_tiles
+obm_tiles_db_username = consumer
+obm_tiles_db_table = obm_built_area_assessments
+obm_tiles_source_id = 1
+
 [SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5]
 ghs_input_filename = FILENAME.hdf5
 

SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles.py

+"""
+Copyright (C) 2021
+  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/>.
+
+
+Global Dynamic Exposure Model
+Helmholtz Centre Potsdam
+GFZ German Research Centre for Geosciences
+Section 2.6: Seismic Hazard and Risk Dynamics
+
+SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles
+========================================================
+
+"""
+
+import sys
+import os
+import shapely
+import datetime
+import h5py
+import getpass
+import psycopg2
+import numpy as np
+import GDE_TOOLS_psql as gdet_psql
+import GDE_TOOLS_world_grid as gdet_grid
+import GDE_TOOLS_general as gdet_gral
+import GDE_TOOLS_read_config_file as gdet_conf
+
+
+def run_this_file(config_dict):
+    ####################################################
+    # READ CONFIGURATION PARAMETERS
+    ####################################################
+    print('Processing configuration parameters...')
+    # Path for output:
+    out_path = config_dict['File Paths']['out_path']
+    # Path for GHS HDF5 file:
+    ghs_path = config_dict['File Paths']['ghs_path'] 
+    # Database from which the GHS built-up values are read:
+    self_dict = config_dict['SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles']
+    obm_tiles_db_hostname = self_dict['obm_tiles_db_hostname']
+    obm_tiles_db_name = self_dict['obm_tiles_db_name']
+    obm_tiles_db_username = self_dict['obm_tiles_db_username']
+    obm_built_area_assessments = self_dict['obm_tiles_db_table']
+    obm_tiles_source_id = int(self_dict['obm_tiles_source_id'])
+    # Tiles database:
+    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'] 
+    # This code will go cell by cell. Define the list of cells from the config file:  
+    grid_cells_list, grid_cells_def_method = gdet_gral.define_grid_cell_list_to_process(
+        config_dict,
+        'Cells to Process')
+    # Cell IDs with numbers below this one will be ignored
+    # (useful while running pieces of regions):
+    if 'min_grid_cell_id' in config_dict['Cells to Process']:
+        min_grid_cell_id = int(config_dict['Cells to Process']['min_grid_cell_id'])
+    else:
+        min_grid_cell_id= np.min(grid_cells_list) - 1
+    ####################################################
+    # START
+    ####################################################
+    log = [] # list that I fill in for logging purposes and then dump to a txt file
+    if grid_cells_def_method == 'country_in_DB':
+        print('Defining grid cells list using method %s... Country: %s'
+              % (grid_cells_def_method, config_dict['Cells to Process']['country_name']))
+    elif grid_cells_def_method == 'admin_ID_in_country': 
+        print('Defining grid cells list using method %s... Country: %s... Admin ID: '
+              % (grid_cells_def_method,
+                 config_dict['Cells to Process']['country_name'],
+                 config_dict['Cells to Process']['admin_id_in_country']))
+    elif grid_cells_def_method == 'random_from_country':
+        print('Defining grid cells list using method %s... Country: %s... '
+              'Number of random cells: %s'
+              % (grid_cells_def_method,
+                 config_dict['Cells to Process']['country_name'],
+                 config_dict['Cells to Process']['num_random_cells']))
+    else:
+        print('Defining grid cells list using method %s...'
+              % (grid_cells_def_method)) 
+    # Connect to the database:
+    connection = psycopg2.connect('host=%s dbname=%s user=%s password=%s'
+                                  % (obm_tiles_db_hostname,
+                                     obm_tiles_db_name,
+                                     obm_tiles_db_username,
+                                     getpass.getpass()))
+    connection.set_session(autocommit=True)
+    # Open a cursor to perform database operations
+    cursor = connection.cursor()
+    for i, grid_cell_id in enumerate(grid_cells_list):
+        if grid_cell_id > min_grid_cell_id:
+            print('\r Working on quadkey %s. Tile %s  of %s.'
+                  % (str(grid_cell_id), str(i+1), str(len(grid_cells_list))), end='')
+            # Coordinates of this particular cell:
+            coords = gdet_grid.get_coordinates_from_cell_id(str(grid_cell_id))
+            lon_w_i, lat_s_i, lon_e_i, lat_n_i = coords
+            # Area of the whole cell:
+            cell_poly = shapely.geometry.Polygon([(lon_w_i,lat_s_i),
+                                                 (lon_e_i,lat_s_i),
+                                                 (lon_e_i,lat_n_i),
+                                                 (lon_w_i,lat_n_i)])
+            area_i = gdet_grid.get_area_polygon_on_Earth(cell_poly, units='km2')
+            # Built-up area in m2:
+            ghs_built_up_i = gdet_psql.get_ghs_built_up_area_from_obm_tiles_opened_connection(
+                str(grid_cell_id),
+                cursor,
+                obm_built_area_assessments,
+                obm_tiles_source_id)
+            print(str(grid_cell_id)+": "+str(ghs_built_up_i))
+            # Write GHS built up area to database (distributed proportionally to the area of
+            # the cell pieces):
+            error_str_i = gdet_psql.write_parameter_proportional_to_area(DB_name_grid,
+                                                                        DB_username_grid,
+                                                                        DB_schema_name_grid,
+                                                                        DB_table_name_grid,
+                                                                        grid_cell_id, area_i, 
+                                                                        ghs_built_up_i/1000.0,
+                                                                        'ghs_km2')
+            if error_str_i != '':
+                log.append(error_str_i)
+    # Close connection to the database:
+    connection.close()
+    # Dump the log list to a file:
+    now = datetime.datetime.now()
+    now_str = "%s_%s_%s_%s_%s_%s" \
+              % (str(now.year),
+                 str(now.month),
+                 str(now.day),
+                 str(now.hour),
+                 str(now.minute),
+                 str(now.second))
+    gdet_gral.write_log_file(log,
+                             os.path.join(out_path,
+                                          'log_file_adding_GHS_%s.txt' % (now_str)))
+    print('\n')
+    print('Done!')
+
+
+
+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; 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',
+                                'Cells to Process',
+                                os.path.basename(__file__).split('.')[0]]
+    config_dict= gdet_conf.read_config_parameters(os.path.join(os.getcwd(),
+                                                               config_filename),
+                                                  section_names_to_validate)
+    run_this_file(config_dict)

docs/guide_to_use_of_config_file.csv

@@ -18,7 +18,8 @@ SERA_exploration_investigate_macro_taxonomies,Y,Y,N,N,N,N,N,N
 SERA_exploration_longest_strings,N,Y,N,N,N,N,N,N
 SERA_exploration_summarise_columns_in_full_CSV_files,Y,Y,N,N,N,N,N,N
 SERA_mapping_admin_units_to_cells,Y,Y,N,Y,N,N,N,N
-SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5,N,Y,N,Y,N,N,N,Y
+SERA_mapping_admin_units_to_cells_add_GHS_from_HDF5,Y,Y,N,Y,N,N,N,Y
+SERA_mapping_admin_units_to_cells_add_GHS_from_OBM_tiles,Y,Y,N,Y,N,N,N,Y
 SERA_mapping_admin_units_to_cells_add_GPW,N,Y,N,Y,N,N,N,Y
 SERA_mapping_admin_units_to_cells_add_Sat,Y,Y,N,Y,N,N,N,Y
 SERA_testing_compare_visual_output_vs_OQ_input_files,Y,Y,N,N,N,N,N,N