Skip to content

Commit

Permalink
update F3 script
Browse files Browse the repository at this point in the history
  • Loading branch information
@matamadio
matamadio committed Feb 12, 2024
1 parent 9c951d8 commit 1b05882
Show file tree
Hide file tree
Showing 3 changed files with 22 additions and 22 deletions.
2 changes: 1 addition & 1 deletion Top-down/FathomV3/.env
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

# Location of data files
# Use absolute paths, and keep the trailing slash
DATA_DIR = X:/Work/Geodata/CCDR/Python/data
DATA_DIR = X:/YourWorkDir/data

# Location to store downloaded rasters and other data
# for the analysis notebooks
Expand Down
24 changes: 11 additions & 13 deletions Top-down/FathomV3/FathomV3.py
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,23 @@
# Importing the required packages
import sys
from common import *
from damageFunctions import damage_factor_FL_builtup, damage_factor_TC_builtup, damage_factor_agri, mortality_factor
from damageFunctions import mortality_factor, damage_factor_builtup, damage_factor_agri
from runAnalysis import run_analysis

# Defining the main function
def main():
# Defining the initial parameters
country = 'MYS' # ISO3166-a3 code
haz_cat = 'FL' # Function analysis:'FL' for floods; 'CF' for coastal floods; 'TC' for tropical cyclones; 'HS' for heat stress
# Class analysis: 'LS' for landslides. 'HS' for heat stress.
return_periods = [5, 10, 20, 50, 75, 100, 200, 250, 500, 1000] # example for Fathom FL [5, 10, 20, 50, 75, 100, 200, 250, 500, 1000]
# example for Deltares CF [5, 10, 25, 50, 100, 250]
# example for Storm TC [10, 20, 50, 100, 200, 500]
# example for Vito HS [5, 20, 100]
min_haz_slider = 0.1 # FL 0.05 # TC 25.0
exp_cat_list = ['POP', 'BU', 'AGR'] # ['POP', 'BU', 'AGR']
exp_nam_list = ['POP', 'BU', 'AGR'] # if None, the default applies: 'Population':'POP', 'Built-up':'BU', 'Agricultural land':'AGR' - If not None, expect a list of same length of exp_cat_list
country = 'MYS' # ISO3166-a2 code
haz_cat = 'COASTAL_UNDEFENDED' # Hazard type:'FLUVIAL_UNDEFENDED'; 'FLUVIAL_DEFENDED', 'PLUVIAL_DEFENDED'; 'COASTAL_UNDEFENDED'; 'COASTAL_DEFENDED'
period = '2050' # Period of the analysis: '2020', '2030', '2050', '2080'
scenario = 'SSP3_7.0' # Climate scenario: 'SSP1_2.6', 'SSP2_4.5', 'SSP3_7.0', 'SSP5_8.5'
return_periods = [5, 10, 20, 50, 100, 200, 500, 1000] # example for Fathom [5, 10, 20, 50, 100, 200, 500, 1000]
exp_cat_list = ['POP', 'BU'] # ['POP', 'BU', 'AGR']
exp_nam_list = ['POP', 'BU'] # if None, the default applies: 'Population':'POP', 'Built-up':'BU', 'Agricultural land':'AGR' - If not None, expect a list of same length of exp_cat_list
adm = 'ADM2' #['ADM1', 'ADM2', 'ADM3']
analysis_app = 'Function' #['Classes', 'Function']
class_edges = [0.05, 0.25, 0.50, 1.00, 2.00] # FL [0.05, 0.25, 0.50, 1.00, 2.00] # TC [17.0, 32.0, 42.0, 49.0, 58.0, 70.0] #HS [23.0, 28.0, 30.0] #LS [0.001, 0.01]
min_haz_slider = 0.1 # Data to be ignored below threshold
class_edges = [0.05, 0.25, 0.50, 1.00, 2.00] # FL [0.05, 0.25, 0.50, 1.00, 2.00]
save_check_raster = False

import time
Expand All @@ -32,7 +30,7 @@ def main():
# Defining the list variable to pass to run_analysis
exp_cat = exp_cat_list[i]
exp_nam = exp_nam_list[i]
run_analysis(country, haz_cat, return_periods, min_haz_slider,
run_analysis(country, haz_cat, period, scenario, return_periods, min_haz_slider,
exp_cat, exp_nam, adm, analysis_app, class_edges, save_check_raster)

print("--- %s seconds ---" % (time.time() - start_time))
Expand Down
18 changes: 10 additions & 8 deletions Top-down/FathomV3/runAnalysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ def zonal_stats_parallel(args):
return zonal_stats_partial(*args)

# Defining the main function to run the analysis
def run_analysis(country: str, haz_cat: str, valid_RPs: list[int],
def run_analysis(country: str, haz_cat: str, period: str, scenario: str, valid_RPs: list[int],
min_haz_threshold: float, exp_cat: str, exp_nam: str, adm_name: str,
analysis_type: str, class_edges: list[float],
save_check_raster: bool, n_cores=None):
Expand All @@ -34,6 +34,8 @@ def run_analysis(country: str, haz_cat: str, valid_RPs: list[int],
----------
country : country ISOa2 code
haz_cat : hazard category
period: time period
scenario: SSP scenario
valid_RPs : return period values to be considered
min_haz_threshold : minimum value for hazard values
exp_cat : exposure category
Expand All @@ -51,7 +53,7 @@ def run_analysis(country: str, haz_cat: str, valid_RPs: list[int],

# Defining the location of administrative, hazard and exposure folders
adm_folder = f"{DATA_DIR}/ADM" # Administrative (gpkg) folder
haz_folder = f"{DATA_DIR}/HZD" # Hazard folder
haz_folder = f"{DATA_DIR}/HZD/{country}/{haz_cat}/{period}/{scenario}" # Hazard folder
exp_folder = f"{DATA_DIR}/EXP" # Exposure folder

# If the analysis type is "Classes", then make sure that the user-defined classes are valid
Expand Down Expand Up @@ -244,7 +246,7 @@ def calc_imp_RPs(RPs, haz_folder, analysis_type, country, haz_cat, exp_cat, exp_
# We reproject using WarpedVRT as this applies the operation from disk
# https://github.com/corteva/rioxarray/discussions/207
# https://rasterio.readthedocs.io/en/latest/api/rasterio.vrt.html
with rasterio.open(os.path.join(haz_folder, f"{country}_{haz_cat}_RP{rp}.tif")) as src:
with rasterio.open(os.path.join(haz_folder, f"1in{rp}.tif")) as src:
vrt_options = {
'src_crs': src.crs,
'crs': exp_data.rio.crs,
Expand All @@ -257,7 +259,7 @@ def calc_imp_RPs(RPs, haz_folder, analysis_type, country, haz_cat, exp_cat, exp_
haz_data.rio.write_nodata(-1.0, inplace=True)

except rasterio._err.CPLE_OpenFailedError:
raise IOError(f"Error occurred trying to open raster file: {country}_{haz_cat}_RP{rp}.tif")
raise IOError(f"Error occurred trying to open raster file: 1in{rp}.tif")

# Set values below min threshold to nan
haz_data = haz_data.where(haz_data.data > min_haz_threshold, np.nan)
Expand All @@ -267,7 +269,7 @@ def calc_imp_RPs(RPs, haz_folder, analysis_type, country, haz_cat, exp_cat, exp_
# Assign impact factor (this is F_i in the equations)
haz_data = damage_factor(haz_data)
if save_check_raster:
haz_data.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{rp}_{exp_cat}_haz_imp_factor.tif"))
haz_data.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{period}_{scenario}_{rp}_{exp_cat}_haz_imp_factor.tif"))
elif analysis_type == "Classes":
# Assign bin values to raster data - Follows: x_{i-1} <= x_{i} < x_{i+1}
bin_idx = np.digitize(haz_data, bin_seq)
Expand All @@ -277,7 +279,7 @@ def calc_imp_RPs(RPs, haz_folder, analysis_type, country, haz_cat, exp_cat, exp_
affected_exp = exp_data.where(haz_data.data > 0, np.nan)

if save_check_raster:
affected_exp.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{exp_cat}_{rp}_affected.tif"))
affected_exp.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{period}_{scenario}_{rp}_{exp_cat}_affected.tif"))

# Conduct analyses for classes
if analysis_type == "Classes":
Expand All @@ -303,7 +305,7 @@ def calc_imp_RPs(RPs, haz_folder, analysis_type, country, haz_cat, exp_cat, exp_
impact_exp = affected_exp.data * haz_data
# If save intermediate to disk is TRUE, then
if save_check_raster:
impact_exp.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{exp_cat}_{rp}_impact.tif"))
impact_exp.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{period}_{scenario}_{rp}_{exp_cat}_impact.tif"))
# Compute the impact per ADM level
impact_exp_per_ADM = gen_zonal_stats(vectors=adm_data["geometry"], raster=impact_exp.data, stats=["sum"],
affine=impact_exp.rio.transform(), nodata=np.nan)
Expand Down Expand Up @@ -357,7 +359,7 @@ def calc_EAEI(result_df, RPs, prob_RPs_df, method, analysis_type, country, haz_c
result_df[f"RP{rp}_EAI_tmp"] = result_df[f"RP{rp}_{exp_cat}_imp"] * freq
if save_check_raster:
EAI_i = impact_exp * freq
EAI_i.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{exp_cat}_{rp}_EAI.tif"))
EAI_i.rio.to_raster(os.path.join(OUTPUT_DIR, f"{country}_{haz_cat}_{period}_{scenario}_{rp}_{exp_cat}_EAI.tif"))

# Computing the EAE or EAI, if probabilistic (len(valid_RPs)>1)
if n_valid_RPs_gt_1:
Expand Down

0 comments on commit 1b05882

Please sign in to comment.