2.4. Case Study Solution#
This notebook presents the solution to the fictional case study on the earthquake impact in San Juan, Argentina.
The goal is to provide the estimated population affected within a 15 km buffer around the epicenter, the number of healthcare facilities located within that area, and a visual representation of the affected zone, including population distribution and health facilities by H3 level 6 hexagons.
import geopandas as gpd
from shapely.geometry import Point, Polygon
import h3
import rasterio
from rasterio.mask import mask
import matplotlib.pyplot as plt
import folium
import branca.colormap as cm
import numpy as np
import tqdm
2.4.1. Estimate population within 15 km of the earthquake epicenter#
We start by defining the epicenter as a point and projecting it to a metric CRS (for San Juan, EPSG:5344) so we can create a 15 km buffer in meters.
# Define the epicenter point
epicenter = gpd.GeoSeries([Point(-68.529, -31.553)], crs="EPSG:4326")
# Project to the correct CRS to build the buffer
epicenter_utm = epicenter.to_crs(epsg=5344)
# Create 15 km buffer (in meters)
buffer = epicenter_utm.buffer(15000)
Next, we load the WorldPop raster, reproject the buffer to the raster’s CRS, mask the raster with the buffer, and sum valid pixel values to estimate the affected population.
path_to_files = '/home/sol/gitrepo/alternative-data-for-crisis/data/review-notebooks/'
with rasterio.open(path_to_files + "arg_ppp_2020_constrained.tif") as src: # Load the raster
# Get raster CRS and reproject buffer to match
buffer_in_raster_crs = buffer.to_crs(src.crs)
# Mask raster with buffer
out_image, _ = mask(src, buffer_in_raster_crs.geometry, crop=True)
data = out_image[0] # Use first band
# Sum all positive values
affected_population = data[data > 0].sum()
print(f"Estimated affected population within 15 km from the epicenter: {int(affected_population):,}")
Estimated affected population within 15 km from the epicenter: 563,361
The same result can be achived using rasterstats package.
from rasterstats import zonal_stats
raster_stats = zonal_stats(buffer_in_raster_crs,
path_to_files + "arg_ppp_2020_constrained.tif",
stats = ['sum'])
print(f"Estimated affected population within 15 km from the epicenter: {raster_stats[0]['sum']}")
Estimated affected population within 15 km from the epicenter: 563361.625
2.4.2. Calculate the number of healthcare facilities within 15 km of the epicenter#
To estimate the number of affected healthcare facilities, we load the shapefile with healthcare locations, reproject the points to match the buffer CRS, and check which facilities fall within the buffer.
# Load the facilities shapefile
facilities_gdf = gpd.read_file(path_to_files + "health_care_facilities/shapefiles/argentina.shp")
# Reproject the facilities to the same CRS as the buffer
facilities_proj = facilities_gdf.to_crs(epsg=5344)
# Select facilities within the 15 km buffer (buffer is a GeoSeries with one geometry; we use .iloc[0] to extract the polygon)
facilities_within_buffer = facilities_proj[facilities_proj.within(buffer.iloc[0])]
print(f"Number of healthcare facilities within 15 km buffer: {len(facilities_within_buffer)}")
Number of healthcare facilities within 15 km buffer: 294
Another approach is to use a spatial join between the projected facilities and the buffer. This may be more efficient for large datasets or multiple buffers.
# Convert the buffer to a GeoDataFrame
buffer_gdf = gpd.GeoDataFrame(geometry=buffer, crs="EPSG:5344")
# Perform spatial join
facilities_within_buffer = gpd.sjoin(facilities_proj, buffer_gdf, predicate="within", how="inner")
print(f"Number of healthcare facilities within 15 km buffer: {len(facilities_within_buffer)}")
Number of healthcare facilities within 15 km buffer: 294
2.4.3. Create an interactive map showing the 15 km buffer, H3-level population, and healthcare facilities#
Before developing the map, we first need to calculate the population for each hexagon that intersects the buffer. We start by loading the GeoJSON file with the hexagons and reprojecting it to match the raster’s CRS.
# Load hexagon layer
hex_gdf = gpd.read_file(path_to_files + "san_juan_hexagons.geojson")
# Open population raster
raster_path = path_to_files + "arg_ppp_2020_constrained.tif"
src = rasterio.open(raster_path)
# Reproject hexagons to match raster CRS
hex_gdf = hex_gdf.to_crs(src.crs)
Next, we filter the hexagons that intersect the buffer. For each of these hexagons, we mask the raster within its geometry and sum the valid pixel values to estimate the population.
If a hexagon has no valid raster data, its population is set to 0.
# Reproject buffer to raster CRS
buffer_proj = buffer_gdf.to_crs(src.crs)
# Select only hexagons intersecting the buffer
hex_buffer = hex_gdf[hex_gdf.intersects(buffer_proj.geometry.iloc[0])].copy()
pop_list = []
with rasterio.open(raster_path) as src:
for geom in hex_buffer.geometry:
out_image, _ = mask(src, [geom], crop=True)
data = out_image[0]
masked_data = np.ma.masked_equal(data, src.nodata)
pop_list.append(masked_data.sum() if masked_data.count() > 0 else 0)
hex_buffer['population'] = pop_list
Another approach is by using rasterstats packages. The value of the above method is that it allows you to see the complete process while rasterstats is more like a blackbox.
raster_stats = zonal_stats(hex_buffer,
path_to_files + "arg_ppp_2020_constrained.tif",
stats = ['sum'])
raster_stats
[{'sum': 5924.9931640625},
{'sum': 4297.50634765625},
{'sum': 8240.6806640625},
{'sum': 638.0714721679688},
{'sum': 3211.880859375},
{'sum': 19055.7578125},
{'sum': 1944.790283203125},
{'sum': 2416.74609375},
{'sum': 14723.875},
{'sum': 10910.8349609375},
{'sum': 98.37982177734375},
{'sum': 183.72927856445312},
{'sum': 123721.984375},
{'sum': 8337.365234375},
{'sum': 15684.705078125},
{'sum': 1122.9620361328125},
{'sum': 15220.75390625},
{'sum': 2888.13232421875},
{'sum': 3384.66259765625},
{'sum': 100496.921875},
{'sum': 3912.71240234375},
{'sum': 23364.0390625},
{'sum': 404.011962890625},
{'sum': 18881.255859375},
{'sum': 114175.8125},
{'sum': 60529.67578125},
{'sum': 3454.65478515625},
{'sum': 10282.654296875},
{'sum': 7171.0439453125},
{'sum': 2908.416015625},
{'sum': 4802.18896484375},
{'sum': 119.06794738769531}]
Finally, we use Folium to build an interactive map that shows:
A 15 km buffer around the epicenter
Hexagons intersecting the buffer colored by estimated population
Healthcare facilities inside the buffer
Population values are shown using a color scale.
# Reproject to EPSG=4326 for mapping
hex_gdf_latlon = hex_buffer.to_crs(epsg=4326)
buffer_latlon = gpd.GeoDataFrame(geometry=buffer, crs="EPSG:5344").to_crs(epsg=4326)
facilities_within_latlon = facilities_within_buffer.to_crs(epsg=4326)
# Initialize the map
m = folium.Map(location=[-31.553, -68.529], zoom_start=11, tiles='CartoDB positron')
# Color scale for population
max_pop = hex_gdf_latlon['population'].max()
colormap = cm.linear.YlOrRd_09.scale(0, max_pop)
colormap.caption = 'Population per H3 Hexagon'
# Hexagons with population
folium.GeoJson(
hex_gdf_latlon,
name='Population by H3 Hexagons',
style_function=lambda feature: {
'fillColor': colormap(feature['properties']['population']),
'color': 'grey',
'weight': 0.5,
'fillOpacity': 0.7 if feature['properties']['population'] > 0 else 0.1
},
tooltip=folium.GeoJsonTooltip(fields=['h3_index', 'population'], aliases=['H3 Index', 'Population'])
).add_to(m)
# Buffer
folium.GeoJson(
buffer_latlon,
name='15 km Buffer',
style_function=lambda x: {
'fill': False,
'color': "#d11d1d",
'weight': 4,
'opacity': 1
}
).add_to(m)
# Healthcare facilities
facilities_layer = folium.FeatureGroup(name="Healthcare Facilities")
for _, row in facilities_within_latlon.iterrows():
folium.CircleMarker(
location=[row.geometry.y, row.geometry.x],
radius=4,
color="#ED8F32",
fill=True,
fill_opacity=0.8,
tooltip=row.get('name', 'Unnamed Facility')
).add_to(facilities_layer)
facilities_layer.add_to(m)
# Epicenter
folium.Marker(
location=[-31.553, -68.529],
popup="Epicenter",
icon=folium.Icon(color="red", icon="info-sign")
).add_to(m)
# Legend and layer control
colormap.add_to(m)
folium.LayerControl().add_to(m)
m