Flood PAHM Disruption - MWI | raster analysis

Flood PAHM Disruption - MWI | raster analysis#

Health Emergencies Preparedness and Response Program (HEPR)#

This notebook is a template workflow to collect data and prepare the main data to perform a baseline physical accessibility analysis to health facilities.
It uses various tools developed by the World Bank’s Geospatial Operations Support Team (GOST).

This notebook focuses on a raster-based implementation of market access, using the motorized Global Friction Surface from the Malaria Atlas Project.
Additionaly, it uses population data from World Pop (Unconstrained UN-Adjusted 2020, 1km resolution).

Data Download Links#

Setup#

Import packages required for the analysis

# System
import sys
import os
from os.path import join, expanduser
from pathlib import Path
# Avoid warnings to pop up
import warnings
warnings.filterwarnings('ignore')

# Visualization tools
# import folium as flm
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.gridspec as gridspec
from rasterio.plot import plotting_extent
from rasterio.plot import show
from mpl_toolkits.axes_grid1 import make_axes_locatable
import contextily as ctx
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import seaborn as sns

os.environ['CARTOPY_USER_BACKGROUNDS'] = '/home/jupyter-wb618081/Python/Backgrounds/'

# Processing 
import numpy as np
import geopandas as gpd
import pandas as pd
from gadm import GADMDownloader

# Raster
import rasterio as rio
from rasterio.features import shapes
from shapely.geometry import box
from rasterio.features import geometry_mask
from rasterstats import zonal_stats
from shapely.geometry import Polygon, box, Point
from shapely.geometry import mapping
import skimage.graph as graph
from scipy.signal import convolve2d

# Graph
import pickle
import networkx as nx
import osmnx as ox

# for facebook data
from pyquadkey2 import quadkey

# Climate/Flood
# import xarray as xr

# Define your path to the Repositories

sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'gostrocks', 'src'))
sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'GOSTNets_Raster', 'src'))
sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'GOSTnets'))
sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'GOST_Urban', 'src', 'GOST_Urban'))
sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'health-equity-diagnostics', 'src', 'modules'))
sys.path.append(join(expanduser("/home/jupyter-wb618081"), 'Repos', 'INFRA_SAP'))

import GOSTnets as gn
from GOSTnets.load_osm import *
import GOSTRocks.rasterMisc as rMisc
from GOSTRocks.misc import get_utm
import GOSTNetsRaster.market_access as ma
import UrbanRaster as urban

from infrasap import aggregator
from infrasap import osm_extractor as osm   

from utils import download_osm_shapefiles

# auto reload
%load_ext autoreload
%autoreload 2

Define below the local folder where you are located

data_dir = join(expanduser("/home/jupyter-wb618081"), 'data')
scratch_dir = join(expanduser("/home/jupyter-wb618081"), 'Health-Access-Metrics')
out_path = join(expanduser("/home/jupyter-wb618081"), 'Health-Access-Metrics', 'Output')

Data Preparation#

Administrative boundaries#

epsg = "EPSG:4326"
epsg_utm = "EPSG:32736"

country = 'Malawi'
iso = 'MWI'
downloader = GADMDownloader(version="4.0")
adm0 = downloader.get_shape_data_by_country_name(country_name=country, ad_level=0)
adm1 = downloader.get_shape_data_by_country_name(country_name=country, ad_level=1)
adm2 = downloader.get_shape_data_by_country_name(country_name=country, ad_level=2)

# iso = 'MWI'
# adm0_path = join(expanduser("R:/"), 'Data', 'GLOBAL/ADMIN', f'Admin0_Polys.shp')
# adm0 = gpd.read_file(adm0_path)
# adm0 = adm0[adm0["ISO3"] == "MWI"].to_crs(4326)

# iso = 'MWI'
# adm2_path = join(expanduser("R:/"), 'Data', 'GLOBAL/ADMIN', f'Admin2_Polys.shp')
# adm2 = gpd.read_file(adm2_path)
# adm2 = adm2[adm2["ISO3"] == "MWI"].to_crs(4326)

Population (origin)#

# wp_path = join(expanduser("R:/"), 'Data', 'GLOBAL/Population/WorldPop_PPP_2020/MOSAIC_ppp_prj_2020', f'ppp_prj_2020_{iso}.tif') # Download from link above
wp_path = join(data_dir, f'ppp_2020_1km_Aggregated.tif') # Download from link above

pop_surf = rio.open(wp_path)

Health Facilities (destinations)#

hf_path = join(data_dir,'MWI','HF_Malawi.xlsx')
df_hf = pd.read_excel(hf_path)

display('The following categories and numbers of Health Facilities are considered to perform the analysis: ')
display(df_hf["Facility Type"].value_counts())

# Consider all health facilities and hospitals

df_hf_hosp = df_hf.loc[df_hf['Facility Type'] == "Hospital"]

'The following categories and numbers of Health Facilities are considered to perform the analysis: '

Facility Type
Outreach          5090
Village Clinic    3542
Health Centre      542
Health Post        152
Dispensary          87
Hospital            85
Name: count, dtype: int64

# Convert from pandas.Dataframe to Geopandas.dataframe

geodf_hf = gpd.GeoDataFrame(
    df_hf, geometry=gpd.points_from_xy(df_hf.Eastings, df_hf.Northings), crs=epsg
)

geodf_hf_hosp = gpd.GeoDataFrame(
    df_hf_hosp, geometry=gpd.points_from_xy(df_hf_hosp.Eastings, df_hf_hosp.Northings), crs=epsg
)

# Clean the geodf
geodf_hf = geodf_hf[['Facility Name', 'Facility Type','District', 'TA', 'geometry']]; geodf_hf.loc[:, 'ID'] = df_hf.index
geodf_hf_hosp = geodf_hf_hosp[['Facility Name', 'Facility Type', 'District', 'TA', 'geometry']]; geodf_hf_hosp.loc[:, 'ID'] = df_hf_hosp.index

Assure correspondence of ADM1 names in Health facilities and official Administrative Units

geodf_hf.rename(columns={'District': 'ADM1', 'TA': 'ADM2'}, inplace=True)
geodf_hf_hosp.rename(columns={'District': 'ADM1', 'TA': 'ADM2'}, inplace=True)
adm1.rename(columns={"NAME_1":"ADM1"}, inplace=True) 

# The following ADM1 are not corresponding
miss_adm1 = np.setdiff1d(np.sort(geodf_hf.ADM1.unique()), np.sort(adm1.ADM1.unique()))
display(miss_adm1)
display(np.setdiff1d(np.sort(adm1.ADM1.unique()), np.sort(geodf_hf.ADM1.unique())))

array(['Mzimba North', 'Mzimba South'], dtype=object)

array(['Mzimba'], dtype=object)

for adm in miss_adm1:
    for idx in (geodf_hf[geodf_hf.ADM1 == adm].index):
            geodf_hf.loc[idx, 'ADM1'] = 'Mzimba'

    for idx_hosp in (geodf_hf_hosp[geodf_hf_hosp.ADM1 == adm].index):
            geodf_hf_hosp.loc[idx_hosp, 'ADM1'] = 'Mzimba'

Flood#

Here, we import Fathom flood data (.tif) of fluvial floods with different return periods.
This data represent and mimic the climate impact on infrastructure and the disruption of the accessibility to health facilities.

# Import multiple rasterio .tif file as a dictionary
# Keys are return periods
# Values are rasterio arrays

# inland waters and oceans: 999
# not-flooded areas: -999 (Fluvial)
# not-flooded areas: 0 (Pluvial)
# Other values represent the flood depth (in m)

flood_fluvial_path = join(data_dir, iso,'FLOOD_SSBN','fluvial_undefended')
flood_pluvial_path = join(data_dir, iso,'FLOOD_SSBN','pluvial')

files=os.listdir(flood_fluvial_path)
flood_dict_fluvial = {}
for file in files:
    key = file.split('_')[1].split('.')[0]
    value = rio.open(join(flood_fluvial_path,file)) #.read(1)
    flood_dict_fluvial[key] = value

files=os.listdir(flood_pluvial_path)
flood_dict_pluvial = {}
for file in files:
    key = file.split('_')[1].split('.')[0]
    value = rio.open(join(flood_pluvial_path,file)) #.read(1)
    flood_dict_pluvial[key] = value

# Preserve the maximum flood depth
flood_dict = {}
for f,key in enumerate(flood_dict_pluvial.keys()):
    out_flood_path = join(data_dir, iso,'FLOOD_SSBN', 'Fmax_' + key +'.tif')
    if os.path.isfile(out_flood_path):
        value = rio.open(out_flood_path) 
        flood_dict[key] = value
    else:
        out_meta = flood_dict_pluvial[key].meta
        flood_max = np.fmax(flood_dict_fluvial[key].read(1),flood_dict_pluvial[key].read(1))
        flood_dict[key] = flood_max
        # flood_dict[key][flood_dict[key] == 0] = -999
        # Write the output raster
        out_flood_path = join(data_dir, iso,'FLOOD_SSBN', 'Fmax_' + key +'.tif')
        with rio.open(out_flood_path, 'w', **out_meta) as dst:
            dst.write(flood_max, 1)
        # Read the output raster
        value = rio.open(out_flood_path) 
        flood_dict[key] = value

Friction Surface#

Process the travel cost surface from the Malaria Atlas Project, clip the raster to our region of interest.

# Only the first time, clip the travel friction surface to the country of interest
out_travel_surface = join(data_dir, iso, f"travel_surface_motorized_{iso}.tif")

if not os.path.isfile(out_travel_surface):
    gfs_path = join(data_dir, '2020_motorized_friction_surface.geotiff')
    gfs_rio = rio.open(gfs_path)
    rMisc.clipRaster(gfs_rio, adm0, out_travel_surface, crop=False)

# Import the clipped friction surface
travel_surf = rio.open(out_travel_surface) #.read(1)

print(travel_surf.res)
print(pop_surf.res)

(0.008333333333333333, 0.008333333333333333)
(0.0083333333, 0.0083333333)

fig, ax = plt.subplots(figsize=(6,8))
ax.set_title("Travel Time", fontsize=12, horizontalalignment='center')
im = ax.imshow(travel_surf.read(1)*1000/60, norm=colors.PowerNorm(gamma=0.5), cmap='viridis')
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size="4%", pad=0.1)
cb = fig.colorbar(im, cax=cax, orientation='vertical')
cb.set_label("(hrs. for 1 km cell)")

../_images/1539e4aee2b877c4a2963ea884820b4946513ff58368ca953dab4aa423d8cc57.png

Preprocessing#

Align the POPULATION & FLOOD raster to the friction surface, ensuring that they have the same extent and resolution.

# If the Standardized data are already present, skip, else generate them

def checkDir(out_path):
    if not os.path.exists(out_path):
        os.makedirs(out_path)

out_pop_surface_std = join(out_path, iso, "WP_2020_1km_STD.tif")
if not os.path.isfile(out_pop_surface_std):
    rMisc.standardizeInputRasters(pop_surf, travel_surf, out_pop_surface_std, resampling_type="nearest")

checkDir(join(scratch_dir, 'data', iso, 'flood'))

for f,key in enumerate(flood_dict.keys()):
    # out_flood_path = join(data_dir, iso,'FLOOD_SSBN', 'Fmax_' + key)
    out_flood_std = join(out_path, iso, 'flood', "STD_" + 'Fmax_'+ key +'.tif')
    if os.path.isfile(out_flood_std):
        None
    else:
        rMisc.standardizeInputRasters(flood_dict[key], travel_surf, out_flood_std, resampling_type="nearest")

Correct the reprojected FLOOD raster extension, ensuring that lands are not covered by inland and open waters

# Import multiple rasterio .tif file as a dictionary
# Keys are return periods
# Values are rasterio arrays

flood_path = join(out_path, iso, 'flood')
files=os.listdir(flood_path)

flood_dict_std = {}
for file in files:
    if file.startswith('STD') == True:
        key = file.split('_')[2].split('.')[0]
        value = rio.open(join(flood_path,file)) #.read(1)
        flood_dict_std[key] = value

display(flood_dict_std["1in10"].read_crs())

CRS.from_epsg(4326)

# # Correct Flood Raster extension, by replacing pixels with value == 999 (inland/open waters) is they overlap a land geometry

# def flood_correction(arr, land_shp, out_path_flood, exclude_value):
#     # arr: rasterio dataset (not readed yet)
    
#     flood_data = arr.read(1).copy()
    
#     # Mask land geometries (True inside)
#     geom_mask = geometry_mask([mapping(geometry) for geometry in land.geometry], 
#                               invert=True, transform=arr.transform, 
#                               all_touched=True, out_shape=(arr.height, arr.width))
    
#     flood_data[geom_mask] = -999
#     correct_mask = np.logical_and(flood_data == -999, arr.read(1) == exclude_value)
    
#     # Create a kernel to look at all neighboring cells
#     kernel = np.ones((3, 3))
#     kernel[1, 1] = 0  # Do not include the cell itself
    
#     # Temporarily replace 999 with np.nan for mean calculation
#     flood_data_with_nan = arr.read(1).copy()
#     flood_data_with_nan[flood_data_with_nan == exclude_value] = np.nan

#     # Temporarily replace np.nan with 0 for convolution and count valid (non-np.nan) neighbors
#     temp_data = np.nan_to_num(flood_data_with_nan, nan=0)
#     valid_neighbors = np.where(flood_data_with_nan != np.nan, 1, 0)

    
#     # Calculate sum and count of neighbors excluding np.nan (previously 999)
#     neighbor_sum = convolve2d(temp_data, kernel, mode='same', boundary='fill', fillvalue=0)
#     neighbor_count = convolve2d(valid_neighbors, kernel, mode='same', boundary='fill', fillvalue=0)
    
#     # Avoid division by zero for locations with all neighbors as np.nan (or originally 999)
#     with np.errstate(divide='ignore', invalid='ignore'):
#         neighbor_mean = neighbor_sum / neighbor_count
#         neighbor_mean[neighbor_count == 0] = np.nan  # Re-assign NaN to locations with no valid neighbors
    
#     # Replace the original values with the computed means where correct_mask is True
#     arr_corrected = arr.read(1).copy()
#     arr_corrected[correct_mask] = neighbor_mean[correct_mask]
    
#     # If necessary, revert np.nan back to 999 in flood_data
#     # flood_data[np.isnan(flood_data)] = exclude_value
    
#     # Save the corrected flood raster
#     out_meta = arr.meta
#     with rio.open(out_path_flood, 'w', **out_meta) as dst:
#             dst.write(flood_data, 1)

#     # Open the corrected flood raster
#     # value = rio.open(out_path_flood) 
#     # flood_data = value
#     return arr_corrected

# land = gpd.read_file(join(data_dir, iso, "geoBoundaries", "geoBoundaries-"+iso+"-ADM1_simplified.shp"))

# flood_dict_correct = {}
# for f,key in enumerate(flood_dict.keys()):
#     out_flood_std = join(out_path, iso, 'flood', "STD_" + 'Fmax_'+ key +'_correct1.tif')
    
#     flood_dict_correct[key] = flood_correction(flood_dict_std[key], land, out_flood_path, 999)

# flood_dict_correct["1in100"] = flood_correction(flood_dict_std["1in100"], land, join(out_path, iso, 'flood', "STD_" + 'Fmax_'+ "1in100" +'_codio.tif'), 999)

Origins#

Prepare a standard grid (pandas.Dataframe) using each cell from the 1km World Pop raster.

pop_surf = rio.open(out_pop_surface_std)
pop = pop_surf.read(1, masked=False)

pop_copy = pop.copy()
pop_copy[pop_copy==0] = np.nan

fig, ax = plt.subplots(figsize=(6,8))
ax.set_title("Population Density", fontsize=12, horizontalalignment='center')
im = ax.imshow(pop_copy, norm=colors.PowerNorm(gamma=0.5), cmap='viridis')
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size="4%", pad=0.1)
cb = fig.colorbar(im, cax=cax, orientation='vertical')
cb.set_label("(n° people per 1 km cell) ")

pop_surf.xy(0,0)

(32.67083333333333, -9.3625)

../_images/4c13d939c6d6881d93b89a804b36f23de9ec823d1dfef80dd5a337aad3b16e6c.png

# Create a population df from population surface 
indices = list(np.ndindex(pop.shape))
xys = [Point(pop_surf.xy(ind[0], ind[1])) for ind in indices]
res_df = pd.DataFrame({
    'spatial_index': indices,
    'xy': xys,
    'pop': pop.flatten()
})
res_df['pointid'] = res_df.index
res_df

	spatial_index	xy	pop	pointid
0	(0, 0)	POINT (32.67083333333333 -9.3625)	67.614548	0
1	(0, 1)	POINT (32.67916666666666 -9.3625)	70.830971	1
2	(0, 2)	POINT (32.68749999999999 -9.3625)	96.481911	2
3	(0, 3)	POINT (32.695833333333326 -9.3625)	111.579651	3
4	(0, 4)	POINT (32.70416666666666 -9.3625)	161.928329	4
...	...	...	...	...
363865	(932, 385)	POINT (35.879166666666656 -17.129166666666666)	32.535667	363865
363866	(932, 386)	POINT (35.88749999999999 -17.129166666666666)	24.958359	363866
363867	(932, 387)	POINT (35.89583333333332 -17.129166666666666)	14.960836	363867
363868	(932, 388)	POINT (35.904166666666654 -17.129166666666666)	14.153514	363868
363869	(932, 389)	POINT (35.91249999999999 -17.129166666666666)	13.631968	363869

363870 rows × 4 columns

Flood impact on Health facilities#

Consider the raw Fmax flood rasters, thus exploiting the original FATHOM dataset resolution
For every Standardized Flood layer (Return Period), extract Flood Depth on Health Facilities location

# Extract flood depth at every facility location

for key in flood_dict.keys():
    coords = [(x,y) for x, y in zip(geodf_hf_hosp.geometry.x, geodf_hf_hosp.geometry.y)]
    geodf_hf_hosp[key] = [x[0] for x in flood_dict[key].sample(coords)]

    coords = [(x,y) for x, y in zip(geodf_hf.geometry.x, geodf_hf.geometry.y)]
    geodf_hf[key] = [x[0] for x in flood_dict[key].sample(coords)]

# Identify Flooded and not Flooded Facilities 

flood_geodf_hf_hosp = {}
dry_geodf_hf_hosp = {}
flood_geodf_hf = {}
dry_geodf_hf = {}
for key in flood_dict.keys():
    flood_geodf_hf_hosp[key] = geodf_hf_hosp[geodf_hf_hosp[key] > 0.2]
    dry_geodf_hf_hosp[key] = geodf_hf_hosp[geodf_hf_hosp[key] == 0]
    
    flood_geodf_hf[key] = geodf_hf[geodf_hf[key] > 0.2]
    dry_geodf_hf[key] = geodf_hf[geodf_hf[key] == 0]

Summary statistics#

For every Flood Scenario (RP):

Number and % of Health facilities disrupted by Facility Type
Number and % of Health facilities disrupted by ADM1 Districts

# Reorder Flood Return Periods in geodf columns

def sort_flood_col(gdf, str_start, str_sep):
    # Define str_start and str_sep as the strings that identify the start of column names and the separator with the numerical value
    fd_columns = [col for col in gdf.columns if col.startswith(str_start)]
    non_fd_columns = [col for col in gdf.columns if not col.startswith(str_start)]
    fd_num = [(int(col.split(str_sep)[1]), col) for col in fd_columns]
    fd_col_sort = [col for _, col in sorted(fd_num)]
    new_column_order = non_fd_columns + fd_col_sort
    # Reorder the DataFrame columns
    gdf = gdf[new_column_order]
    return(gdf)

geodf_hf = sort_flood_col(geodf_hf, "1in", "in")
geodf_hf_hosp = sort_flood_col(geodf_hf_hosp, "1in", "in")

# % of Facilities disrupted, by Facility Type

stats1 = geodf_hf[geodf_hf != 0].groupby("Facility Type").count().drop(columns = ["Facility Name", "ADM1", "ADM2", "geometry"]).rename(columns={'ID': 'Total n°'})
fd_columns = [col for col in stats1.columns if col.startswith("1in")]
for col in fd_columns:
    stats1[col] = (stats1[col]/stats1["Total n°"])*100

stats1[fd_columns] = stats1[fd_columns].round(2)
display(stats1)

	Total n°	1in5	1in10	1in20	1in50	1in75	1in100	1in200	1in250	1in500	1in1000
Facility Type
Dispensary	87	6.90	10.34	13.79	16.09	16.09	17.24	20.69	21.84	24.14	25.29
Health Centre	542	6.27	8.49	9.59	12.73	13.65	14.58	16.42	16.61	18.63	19.56
Health Post	152	14.47	18.42	21.71	22.37	23.68	23.68	25.66	26.32	26.97	29.61
Hospital	84	3.57	3.57	3.57	3.57	3.57	3.57	4.76	4.76	7.14	9.52
Outreach	5090	7.70	10.53	12.73	15.30	16.44	17.43	19.17	19.57	21.96	23.63
Village Clinic	3542	7.40	10.19	12.54	15.19	16.29	17.11	19.11	19.34	21.46	23.57

stats1_long = pd.melt(stats1.reset_index().drop(columns = "Total n°"), id_vars="Facility Type", var_name="scen")
scen = stats1_long.scen.unique()

# % of Facilities disrupted, by ADM1 (Districts)

stats2 = geodf_hf[geodf_hf != 0].groupby("ADM1").count().drop(columns = ["Facility Name", "Facility Type", "ADM2", "geometry"]).rename(columns={'ID': 'Total n°'})
fd_columns = [col for col in stats2.columns if col.startswith("1in")]
for col in fd_columns:
    stats2[col] = (stats2[col]/stats2["Total n°"])*100

stats2[fd_columns] = stats2[fd_columns].round(2)

# Merge with ADM1 geometries
stats2 = stats2.merge(adm1[["ADM1", "geometry"]], on='ADM1', how='left').set_index("ADM1")
stats2 =gpd.GeoDataFrame(stats2, geometry=stats2["geometry"], crs="EPSG:4326")
display(stats2)

	Total n°	1in5	1in10	1in20	1in50	1in75	1in100	1in200	1in250	1in500	1in1000	geometry
ADM1
Balaka	320	15.31	20.31	23.75	29.69	30.00	30.31	33.12	33.12	34.06	35.94	MULTIPOLYGON (((35.07923 -15.30382, 35.07925 -...
Blantyre	276	7.61	10.51	11.59	13.04	14.49	15.22	15.94	15.94	17.75	17.75	MULTIPOLYGON (((34.94884 -15.98430, 34.94793 -...
Chikwawa	290	12.41	19.31	23.79	26.90	29.66	34.48	37.24	37.59	43.45	45.86	MULTIPOLYGON (((34.92864 -16.64210, 34.92845 -...
Chiradzulu	156	1.92	7.05	9.62	12.18	12.18	12.18	12.18	12.18	12.18	12.82	MULTIPOLYGON (((35.32407 -15.91358, 35.32412 -...
Chitipa	213	3.76	3.76	4.23	5.63	5.63	6.57	8.92	8.92	8.92	10.80	MULTIPOLYGON (((33.75890 -10.23768, 33.79720 -...
Dedza	389	5.66	7.20	8.23	10.03	10.54	10.80	13.11	13.11	15.68	17.48	MULTIPOLYGON (((34.03575 -14.48153, 34.03575 -...
Dowa	559	2.50	3.04	5.01	7.16	8.23	9.30	9.84	10.38	12.16	13.77	MULTIPOLYGON (((34.18658 -13.69968, 34.18738 -...
Karonga	281	27.76	33.10	38.79	44.48	46.62	48.40	50.89	51.60	53.74	55.87	MULTIPOLYGON (((34.17258 -10.58059, 34.17202 -...
Kasungu	660	0.61	0.91	1.36	2.42	3.64	4.09	4.85	5.00	6.21	6.52	MULTIPOLYGON (((33.36169 -13.58856, 33.36145 -...
Likoma	13	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	MULTIPOLYGON (((34.74764 -12.06736, 34.74708 -...
Lilongwe	572	0.52	0.87	1.75	2.80	2.97	3.67	4.90	4.90	6.82	8.74	MULTIPOLYGON (((33.69741 -14.43242, 33.69709 -...
Machinga	458	10.04	13.32	16.38	20.52	23.58	24.24	26.42	26.86	30.13	32.97	MULTIPOLYGON (((35.34245 -15.19578, 35.34162 -...
Mangochi	472	18.43	23.94	26.48	29.87	31.99	33.47	35.81	37.29	39.83	41.95	MULTIPOLYGON (((35.50280 -14.75274, 35.50346 -...
Mchinji	313	1.28	2.88	4.79	7.67	9.27	9.90	11.50	12.46	16.61	22.68	MULTIPOLYGON (((33.34943 -14.09478, 33.34931 -...
Mulanje	304	10.86	15.79	22.04	25.00	25.99	26.97	32.24	32.24	34.21	35.20	MULTIPOLYGON (((35.34843 -16.17946, 35.34831 -...
Mwanza	72	0.00	0.00	0.00	0.00	0.00	2.78	12.50	12.50	13.89	13.89	MULTIPOLYGON (((34.73567 -15.81710, 34.73438 -...
Mzimba	943	1.38	2.65	3.50	4.56	5.51	5.83	7.00	7.00	9.23	10.60	MULTIPOLYGON (((33.62449 -12.68889, 33.62415 -...
Neno	148	10.14	12.84	13.51	14.86	15.54	16.22	16.22	16.22	16.89	17.57	MULTIPOLYGON (((34.73540 -15.78676, 34.73513 -...
Nkhata Bay	261	6.51	6.51	6.90	9.58	9.58	9.96	10.73	11.11	11.49	14.56	MULTIPOLYGON (((34.01905 -12.23464, 34.01835 -...
Nkhotakota	333	5.71	8.11	9.01	10.21	10.21	10.21	12.61	13.21	15.32	16.22	MULTIPOLYGON (((34.04913 -12.24265, 34.11755 -...
Nsanje	226	18.58	27.43	37.17	42.04	42.92	44.25	45.58	46.90	49.12	50.44	MULTIPOLYGON (((35.28840 -16.75789, 35.28944 -...
Ntcheu	419	7.40	9.07	11.69	13.13	13.84	14.80	16.71	16.95	19.09	19.33	MULTIPOLYGON (((34.84460 -15.22389, 34.84412 -...
Ntchisi	272	0.74	1.10	1.10	2.94	3.31	3.31	3.68	3.68	3.68	5.51	MULTIPOLYGON (((34.04879 -13.53227, 34.04847 -...
Phalombe	191	23.04	25.13	27.75	33.51	35.08	37.70	41.36	41.88	51.31	56.02	MULTIPOLYGON (((35.76186 -15.92668, 35.76157 -...
Rumphi	269	14.50	22.68	27.14	31.60	31.60	31.97	35.32	35.69	39.41	42.75	MULTIPOLYGON (((33.78087 -11.02132, 33.78003 -...
Salima	238	13.03	18.91	24.37	31.51	34.03	36.13	39.92	39.92	42.02	43.28	MULTIPOLYGON (((34.51638 -14.10933, 34.51516 -...
Thyolo	408	3.92	4.41	5.39	6.13	6.37	6.62	6.86	6.86	7.60	8.09	MULTIPOLYGON (((35.11417 -16.31244, 35.11364 -...
Zomba	441	9.52	16.10	17.69	21.54	23.81	25.17	28.34	29.02	32.65	36.51	MULTIPOLYGON (((35.31834 -15.67589, 35.31762 -...

Map HF Impact Results#

fig = plt.figure(figsize=(12, 6))

ax0 = fig.add_subplot(111)
ax0.set_title("Floods impact by Facility Type")
ax0.set_xticklabels(scen)
ax0.set_xlabel("Flood Scenarios (Return Period)")
ax0.set_ylabel("Health Facilities impacted (%)")
ax0.set_ylim(0,40)
ax0.legend(loc='upper left', fontsize = 10, title = "Facility type")
ax0 = sns.barplot(
    data=stats1_long, hue="Facility Type",# errorbar=("pi", 50),
    x="scen", y="value",
    palette='colorblind', alpha=.6, 
)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

../_images/86bfc794d508818809d05ed48741be47f2a648cbce0d93544fe9e12ebc1aa81a.png

figsize = (14, 10)
fig = plt.figure(figsize=figsize)
projection = ccrs.PlateCarree()
gs = gridspec.GridSpec(2, 5)

title = "Flood impact by Districts (ADM1)"
fig.suptitle(title, size=16, y=0.95)

# Define the colormap
cmap = plt.get_cmap('YlOrRd')
# Define the normalization from 0 to 45
norm = colors.Normalize(vmin=0, vmax=45)

for i, flood in enumerate(scen):
    if i < 5:
        ax = fig.add_subplot(gs[0, i], projection=projection)
    else:
        ax = fig.add_subplot(gs[1, i-5], projection=projection)

    ax.set_title(flood)
    ax.get_xaxis().set_visible(True)
    ax.get_yaxis().set_visible(True)

    # Plot the data
    stats2.plot(
        ax=ax, column=flood, cmap=cmap, legend=False,
        alpha=1, linewidth=0.2, edgecolor='black', 
        norm=norm
    )

    ax.background_img(name='NaturalEarthRelief', resolution='high', extent = [32.5, 36, -18, -9])
    
# Add a colorbar to the figure, adjusting its position via `fig.add_axes`
cax = fig.add_axes([0.92, 0.25, 0.015, 0.5])  # Adjust the position and size of the colorbar here
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
fig.colorbar(sm, cax=cax, label="Facilities impacted (%)")

plt.show()
plt.savefig(join(out_path, iso, title + ".png"), dpi=300, bbox_inches='tight', facecolor='white')

../_images/b77ffdb150af483a8759966f8b7d9944ef312c33d3ec83948c304241b8c36659.png

<Figure size 640x480 with 0 Axes>

Flood impact on Friction Surface#

Consider the Standardized Fmax flood rasters, thus exploiting the original FATHOM dataset resolution
We consider different degrees of accessibility/mobility disruption according to different flood depth (FD) levels

If the FD is less than 0.3 meters, than accessibility is preserved but signifcantly impacted (slowered)
If the FD is higher than 0.3 meters, than mobility is interrupted on that specific pathway

Pregnolato et al.PNG

# Import multiple rasterio .tif file as a dictionary
# Keys are return periods
# Values are rasterio arrays

flood_path = join(out_path, iso, 'flood')
files=os.listdir(flood_path)

flood_dict_std = {}
for file in files:
    if file.startswith('STD') == True:
        key = file.split('_')[2].split('.')[0]
        value = rio.open(join(flood_path,file)) #.read(1)
        flood_dict_std[key] = value

display(flood_dict_std["1in10"].read_crs())

CRS.from_epsg(4326)

# Create weights for the Friction Surface
# If water level is less than 30 cm, accessibility is affected (from pregnolato et al., (2017))
# If water level is more than 30 cm, accessibility is disrupted (travel time *50)

def friction_flood_weight(flood_rio):
    ''' Generate flood derived weights for the friction surface 
    
    INPUTS
        flood_rio [rasterio] - raster with flood depths from which to generate weights
  
    RETURNS
        flood_rio [rasterio] - raster with weights 
    ''' 
    
    from scipy.interpolate import interp1d

    # Function for scale Flood Depth array between 1 (min) and 0 (max)
    def scale_array(arr):
        min_val = np.min(arr)
        max_val = np.max(arr)
        scaled_arr = (arr - min_val) / (max_val - min_val)
        scaled_arr = 1 - scaled_arr  # Invert the scaled values
        return scaled_arr
    
    flood_rio = flood_rio.read(1)
    flood_rio_orig = flood_rio.copy()

    # Apply the quadratic (parabolic) relationship from Pregnolato et al., (2017)
    # Describes the vehicles max speed according to Flood Depth values (in mm) 
    flood_arr = flood_rio[(flood_rio_orig > 0) & (flood_rio_orig < 0.3)]*1000   # flood depth from m to mm
    depth_val = np.arange(1,300,0.1)     # Define discrete Flood Depth values between 1 mm and 300 mm
    func = scale_array(9*10**(-4)*(depth_val**2) - 0.5529*depth_val + 86.9448)   # eq.4 from Pregnolato et al., (2017)
    # Scale the function with the maximum vehicles speed
    func = func*func.max()**-1
    interp_func = interp1d(depth_val, func, kind='linear', fill_value='extrapolate')
    flood_arr_weight = interp_func(flood_arr)

    # Flood Depth < 0 or 999 (open waters) -> 1
    # Flood Depth > 30 cm                  -> 50
    # Flood Depth > 0 & < 30 cm            -> inverse of quadratic function (between 1 and 20 circa)
    flood_rio[(flood_rio_orig > 0) & (flood_rio_orig < 0.3)] = 1/flood_arr_weight
    flood_rio[(flood_rio_orig <= 0) | (flood_rio_orig == 999)] = 1
    flood_rio[(flood_rio_orig >= 0.3) & (flood_rio_orig < 999)] = 50

    # To ensure that weights are always >1 (avoid decresing travel time)
    flood_rio[(flood_rio <1)] = 1
    
    del flood_rio_orig
    return flood_rio

flood_dict_fd = flood_dict_std.copy()

for key in flood_dict_fd.keys():

    flood_dict_fd[key] = friction_flood_weight(flood_dict_fd[key])

Create new Friction Surface that accounts for Floods disruption

# Weight the Friction Surface with the discrete flood depth array

travel_surf_flood = dict()
for key in flood_dict_fd.keys():
    travel_surf_flood[key] = travel_surf.read(1) * flood_dict_fd[key]

Create an MCP graph object from the flooded friction surfaces.

# convert friction surface to traversal time (lazily). Original data are minutes to travel 1 m:
# We convert it to minutes to cross the cell (1000m). This could be revised

inG_data = dict()
inG_data["baseline"] = travel_surf.read(1) * 1000
for key in travel_surf_flood.keys():
    inG_data[key] = travel_surf_flood[key] * 1000

# Correct no data values. Not needed but good to check
# inG_data[inG_data < 0] = 99999999
# inG_data[inG_data < 0] = np.nan

mcp = dict()
mcp["baseline"] = graph.MCP_Geometric(inG_data["baseline"])
for key in inG_data.keys():
    mcp[key] = graph.MCP_Geometric(inG_data[key])

# Check descriptive statistics
print(f"Mean: {np.mean(inG_data['1in5'])}")
print(f"Max: {np.max(inG_data['1in5'])}")
print(f"Min: {np.min(inG_data['1in5'])}")
print(f"Std: {np.std(inG_data['1in5'])}")

Mean: 16.180652618408203
Max: 8287.5048828125
Min: 0.5
Std: 91.2672348022461

Data analysis#

Indicators of interest

Percentage of population within 2h of driving to the nearest primary care facility (population level, and by SES quintile).
Percentage of population within 2h of driving to the nearest district hospital (population, and by SES quintile).
Percentage of health facilities with direct access to an all season road.
Percentage of health facilities within 2km of an all season road.

# Calculate the travel time from each grid-cell to the nearest not-flooded destination
# According to the Flood scenario, we'll a different number of flooded facilities  

res_hf = dict()
res_hf_hosp = dict()
res_hf["baseline"] = ma.calculate_travel_time(travel_surf, mcp["baseline"], geodf_hf)[0]
res_hf_hosp["baseline"] = ma.calculate_travel_time(travel_surf, mcp["baseline"], geodf_hf_hosp)[0]
for key in travel_surf_flood.keys():
    res_hf[key] = ma.calculate_travel_time(travel_surf, mcp[key], dry_geodf_hf[key])[0]
    res_hf_hosp[key] = ma.calculate_travel_time(travel_surf, mcp[key], dry_geodf_hf_hosp[key])[0]

# Check dimension
display(len(res_hf[key].flatten()), len(dry_geodf_hf_hosp[key]))

# Decide to consider all facilities or Hospitals
res_df.loc[:, 'tt_hf_' + "baseline"] = res_hf["baseline"].flatten()
res_df.loc[:, 'tt_hosp_' + "baseline"] = res_hf_hosp["baseline"].flatten()
for key in travel_surf_flood.keys():
    res_df.loc[:, 'tt_hf_' + key] = res_hf[key].flatten()
    res_df.loc[:, 'tt_hosp_' + key] = res_hf_hosp[key].flatten()

# remove values where pop is 0 or nan
res_df = res_df.loc[res_df['pop']!=0].copy()
res_df = res_df.loc[~(res_df['pop'].isna())].copy()

res_df.loc[:,'xy'] = res_df.loc[:,'xy'].apply(lambda x: Point(x))

res_df.head(2)

	spatial_index	xy	pop	pointid	tt_hf_baseline	tt_hosp_baseline	tt_hf_1in100	tt_hosp_1in100	tt_hf_1in1000	tt_hosp_1in1000	...	tt_hf_1in250	tt_hosp_1in250	tt_hf_1in20	tt_hosp_1in20	tt_hf_1in5	tt_hosp_1in5	tt_hf_1in50	tt_hosp_1in50	tt_hf_1in200	tt_hosp_1in200
0	(0, 0)	POINT (32.67083333333333 -9.3625)	67.614548	0	205.364383	205.364383	212.504426	212.504426	217.231925	217.231925	...	211.680068	211.680068	209.776903	209.776903	209.552599	209.552599	213.104327	213.104327	211.806125	211.806125
1	(0, 1)	POINT (32.67916666666666 -9.3625)	70.830971	1	203.864383	203.864383	211.004426	211.004426	215.731925	215.731925	...	210.180068	210.180068	208.276903	208.276903	208.052599	208.052599	211.604327	211.604327	210.306125	210.306125

2 rows × 26 columns

Create Geodataframe for the population grid

res_gdf = gpd.GeoDataFrame(res_df, geometry='xy', crs=epsg)
res_gdf.rename(columns={'xy':'geometry'}, inplace=True)
res_gdf.set_geometry('geometry', inplace=True)

# convert travel time from minutes to hours
for col in res_gdf.columns:
    if 'tt_' in col:
        res_gdf.loc[:, col] = res_gdf.loc[:, col] / 60
        
res_gdf.head(2)

	spatial_index	geometry	pop	pointid	tt_hf_baseline	tt_hosp_baseline	tt_hf_1in100	tt_hosp_1in100	tt_hf_1in1000	tt_hosp_1in1000	...	tt_hf_1in250	tt_hosp_1in250	tt_hf_1in20	tt_hosp_1in20	tt_hf_1in5	tt_hosp_1in5	tt_hf_1in50	tt_hosp_1in50	tt_hf_1in200	tt_hosp_1in200
0	(0, 0)	POINT (32.67083 -9.36250)	67.614548	0	3.42274	3.42274	3.54174	3.54174	3.620532	3.620532	...	3.528001	3.528001	3.496282	3.496282	3.492543	3.492543	3.551739	3.551739	3.530102	3.530102
1	(0, 1)	POINT (32.67917 -9.36250)	70.830971	1	3.39774	3.39774	3.51674	3.51674	3.595532	3.595532	...	3.503001	3.503001	3.471282	3.471282	3.467543	3.467543	3.526739	3.526739	3.505102	3.505102