# Import required libraries
from geopandas import read_file, gpd
import folium
from shapely.geometry import Point
from termcolor import colored
# gpd.__version__
# !pip install geopandas==0.10.2
# !pip install folium==0.16.0
# !pip install shapely==1.8.1
# !pip install numpy==1.23.5
# Surge Areas Data
input_file = './Data/StormSurgeRiskAreas/Storm_Surge_Risk_Areas.shp'
surge = gpd.read_file(input_file)
surge['HES_ZONE'].replace(99, 5, inplace=True)
input_file = './Data/Pharmacy_Locations/Pharmacy_Locations.shp'
pharma = gpd.read_file(input_file)
pharma.tail()
# Reproject
print(pharma.crs)
pharma = pharma.to_crs(4326)
Click on the map to select point location for further analysis.
center_lat = pharma.geometry.y.mean()
center_lon = pharma.geometry.x.mean()
m = folium.Map(location=[center_lat, center_lon], tiles="CartoDB positron", zoom_start=12)
pharma.explore(m=m, name="pharmacy locations")
m.add_child(folium.ClickForLatLng(format_str='"[" + lat + "," + lng + "]"', alert=True))
m
def get_category_surge(category):
surge_category = surge[surge['HES_ZONE']<=category]
return surge_category
def get_category_pharma(category):
pharma_category = gpd.clip(pharma, get_category_surge(category))
return pharma_category
def create_buffer(clicked_point):
# Convert the clicked point into a GeoSeries
click_geom = gpd.GeoSeries([Point(clicked_point[1], clicked_point[0])], crs="EPSG:4326")
# Find category of storm that would affect the given location
for idx, row in surge.iterrows():
category = row['HES_ZONE']
if click_geom.geometry[0].within(row['geometry']):
print(colored(f"Your location is affected by storm category {category} or above!", "red", attrs=["bold"]))
break
else:
print(colored("Your location is not affected by storm surge", "green", attrs=["bold"]))
# Plot user-input location
click_geom.explore(m=m, name="selected point", marker_kwds=dict(radius=4), style_kwds={"color":"black", "fillColor":"black"})
# Reproject to UTM for accurate distance-based buffering
click_geom_utm = click_geom.to_crs(pharma.estimate_utm_crs())
# Create a buffer of one mile (m)
buffer = click_geom_utm.buffer(1609) # one mile buffer
# Reproject back to EPSG:4326
buffer_wgs84 = buffer.to_crs(epsg=4326)
# Plot buffer
buffer_wgs84.explore(m=m, name="buffer", style_kwds={"fill":True, "fillOpacity":0.3})
# Find points from pharma within the buffer
points_in_buffer = pharma[pharma.geometry.within(buffer_wgs84.iloc[0])]
# Find points that are inundated
pharma_category = get_category_pharma(category)
buffer_points_in_category = points_in_buffer[points_in_buffer['OBJECTID'].isin(pharma_category['OBJECTID'])]
# Plot points in buffer in green and inundated pharmacies in red
if len(points_in_buffer)>0:
points_in_buffer.explore(m=m, name="pharmacies in buffer zone", style_kwds={"color":"green", "weight":3, "fillColor":"green"})
if len(buffer_points_in_category)>0:
buffer_points_in_category.explore(m=m, name="inundated pharmacies", style_kwds={"color":"red", "weight":3, "fillColor":"red"})
#Plot surge inundated area for the category of storm for given location
surge_category = get_category_surge(category)
surge_category.explore(m=m, name="storm surge areas",
style_kwds=dict(color="black", weight=0.5, opacity=0.3, fillColor='yellow', fillOpacity=0.3))
# Print message if no points found
else:
print('No pharmacies in buffer zone, Expand search zone')
return points_in_buffer, category, m
point = [38.878740,-77.020340]
pts, cat, m = create_buffer(point)
folium.LayerControl().add_to(m)
m
m = folium.Map(location=[center_lat, center_lon], tiles="CartoDB positron", zoom_start=12)
point = [38.939918,-77.041626]
pts, cat, m = create_buffer(point)
folium.LayerControl().add_to(m)
m