Introduction to Python and Geopandas for Flooded Area Analysis - Tutorial

Geopandas is one of the most advanced geospatial libraries in Python because it combines the spatial tools of Shapely, it can create and read different OGC vector spatial data, it can couple the Pandas tools to manage, filter, and make operations over the columns of the metadata, it has the capability to plot geospatial data on Matplotlib and even to Folium among other features. We have developed a tutorial of Geopandas applied to the analysis of flooded areas over the Boise city for a return period of 200 years; the tutorial covers introductory concepts of Geopandas, it will work with point, line and polygon vector data, create plots, simplify vertices and perform geospatial queries over inundated facilities and highways.

Tutorial


Code

import geopandas as gpd
import matplotlib.pyplot as plt
#open flooded areas
floodedAreas = gpd.read_file('../Data/Shps/boiseFlood200yr_v4.shp')
floodedAreas.head()

fid STAGE ELEV USGSID GRIDID QCFS STGRATING Comments Area geometry
0 1.0 14.94 2617.892 13206000 13 23900 EXTRAPOLATED, from USACE MODEL 0.5-percent chance flow 29243133 POLYGON ((-116.30349 43.68274, -116.30350 43.6...
#show columns
floodedAreas.columns
Index(['fid', 'STAGE', 'ELEV', 'USGSID', 'GRIDID', 'QCFS', 'STGRATING',
       'Comments', 'Area', 'geometry'],
      dtype='object')
#show bounds per element
floodedAreas.bounds.head()

minx miny maxx maxy
0 -116.516417 43.662396 -116.298905 43.694833
#plot geometry
floodedAreas.plot()
<AxesSubplot:>
#calculate areas
floodedAreasUtm = floodedAreas.to_crs(32611)
floodedAreasUtm['area'] = floodedAreasUtm.area
floodedAreasUtm['area_km2'] = floodedAreasUtm.area/1000000
floodedAreasUtm.head()

fid STAGE ELEV USGSID GRIDID QCFS STGRATING Comments Area geometry area area_km2
0 1.0 14.94 2617.892 13206000 13 23900 EXTRAPOLATED, from USACE MODEL 0.5-percent chance flow 29243133 POLYGON ((556138.632 4836872.248, 556138.373 4... 2.922113e+07 29.221126

Working with linestrings

primaryHighways = gpd.read_file('../Data/Shps/primaryHighways.gpkg')
primaryHighways.head()

full_id osm_id osm_type highway old_name maxspeed:variable maxspeed:conditional turn:lanes:both_ways turn:lanes turn:lanes:backward ... alt_name old_ref ref oneway maxspeed NHS lanes surface name geometry
0 w13691899 13691899 way primary None None None None None None ... None None None None None None None asphalt North Linder Road LINESTRING (-116.41365 43.67338, -116.41365 43...
1 w13692173 13692173 way primary None None None None None None ... None None None None None None 4 None North Ten Mile Road LINESTRING (-116.43362 43.66294, -116.43362 43...
2 w13692295 13692295 way primary None None None None None None ... None None ID 55 yes 55 mph yes None asphalt North Eagle Road LINESTRING (-116.35424 43.64725, -116.35423 43...
3 w13695920 13695920 way primary None None None None None None ... None US 30 None yes 40 mph None None paved East Fairview Avenue LINESTRING (-116.35642 43.61960, -116.35903 43...
4 w13696901 13696901 way primary None None None None None None ... North Horseshoe Bend Road None ID 55 None 55 mph yes None asphalt North Highway 55 LINESTRING (-116.26604 43.76983, -116.26620 43...

5 rows × 35 columns

#show column names
primaryHighways.columns
Index(['full_id', 'osm_id', 'osm_type', 'highway', 'old_name',
       'maxspeed:variable', 'maxspeed:conditional', 'turn:lanes:both_ways',
       'turn:lanes', 'turn:lanes:backward', 'turn:lanes:forward',
       'official_name', 'maxspeed:type', 'cycleway', 'lanes:forward',
       'lanes:both_ways', 'lanes:backward', 'layer', 'bridge',
       'source:maxspeed', 'source:hgv:national_network',
       'hgv:national_network', 'hgv', 'cycleway:right', 'bicycle', 'alt_name',
       'old_ref', 'ref', 'oneway', 'maxspeed', 'NHS', 'lanes', 'surface',
       'name', 'geometry'],
      dtype='object')
#show roads
primaryHighways['name'].value_counts()
West Chinden Boulevard    38
West Fairview Avenue      34
West State Street         24
North Eagle Road          17
East Chinden Boulevard    16
East Fairview Avenue      12
South Eagle Road          12
North Highway 55          10
South Linder Road         10
East State Street          9
North Ten Mile Road        9
North Linder Road          8
North Glenwood Street      8
North Milwaukee Street     4
North Meridian Road        4
Emmett Highway             3
State Hwy 55               1
North Five Mile Road       1
West Cherry Lane           1
Name: name, dtype: int64
# show figure with plotting options
fig, ax = plt.subplots(figsize=(16,16))
primaryHighways.plot(column='name', ax=ax, legend=True, categorical=True)
<AxesSubplot:>
# show figure with plotting options with flooding area
fig, ax = plt.subplots(figsize=(16,16))
primaryHighways.plot(column='name', ax=ax, legend=True, categorical=True)
floodedAreas.plot(ax=ax, alpha=0.5)
<AxesSubplot:>

Working with polygons

primaryFacilities = gpd.read_file('../Data/Shps/schoolsHospitalsPoliceFire.gpkg')
primaryFacilities.head()

amenity name geometry
0 school Pioneer Elementary School MULTIPOLYGON (((-116.34768 43.64752, -116.3475...
1 school Centennial High School MULTIPOLYGON (((-116.33940 43.65290, -116.3389...
2 school Joplin Elementary School MULTIPOLYGON (((-116.33274 43.65530, -116.3300...
3 school Lowell Scott Middle School MULTIPOLYGON (((-116.35409 43.65206, -116.3539...
4 school Cecil D. Andrus Elementary School MULTIPOLYGON (((-116.34585 43.66075, -116.3447...
#show columns
primaryFacilities['amenity'].unique()
array(['school', 'hospital', 'fire_station'], dtype=object)
#plot schools 
fig, ax = plt.subplots(figsize=(16,16))
ax.axis('equal')
primaryFacilities.plot(column='amenity', ax=ax, legend=True, categorical=True, cmap='plasma')
<AxesSubplot:>
#plot schools with flooded areas
fig, ax = plt.subplots(figsize=(16,16))
ax.axis('equal')
primaryFacilities.plot(column='amenity', ax=ax, legend=True, categorical=True, cmap='plasma')
floodedAreas.plot(ax=ax, alpha=0.5)
<AxesSubplot:>

Spatial operations

#show geoframe crs
print(floodedAreas.crs)
print(primaryHighways.crs)
print(primaryFacilities.crs)
epsg:4326
epsg:4326
epsg:4326
#make composite plot
fig, ax = plt.subplots(figsize=(16,12))

flooded = floodedAreas.plot(ax=ax, alpha=0.5)
highways = primaryHighways.plot(column='name', ax=ax, categorical=True)
schools = primaryFacilities.plot(column='amenity', ax=ax, categorical=True, cmap='plasma')

ax.set_aspect('equal')
plt.show()

Plot data with Folium

All data need to be in geographical coordinates

# plot with folium
import folium 

#convert everythin to WGS84
#floodedAreasGeo = floodedAreas.to_crs(4326)
#primaryHighwaysGeo = primaryHighwaysFix.to_crs(4326)
#schoolPointsGeo = schoolPoints.to_crs(4326)
#we need to have a reference centroid for the map

#dissolve flooded areas
floodedDissolved = floodedAreas.dissolve()
floodedDissolved

geometry fid STAGE ELEV USGSID GRIDID QCFS STGRATING Comments Area
0 POLYGON ((-116.30349 43.68274, -116.30350 43.6... 1.0 14.94 2617.892 13206000 13 23900 EXTRAPOLATED, from USACE MODEL 0.5-percent chance flow 29243133
#get reference coords to initialize the map
refX = floodedDissolved.centroid.x
refY = floodedDissolved.centroid.y

#create map
map = folium.Map(location = [refY,refX], tiles='OpenStreetMap', zoom_start = 13)
map
C:\Users\saulm\AppData\Local\Temp\ipykernel_2280\3149012135.py:2: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  refX = floodedDissolved.centroid.x
C:\Users\saulm\AppData\Local\Temp\ipykernel_2280\3149012135.py:3: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  refY = floodedDissolved.centroid.y
Make this Notebook Trusted to load map: File -> Trust Notebook
#check the number of vertices of lines
import numpy as np

nVertexFlooded = [np.array(row.geometry.exterior.coords).shape[0] for index, row in floodedAreas.iterrows()]
nVertexHighways = [np.array(row.geometry.coords).shape[0] for index, row in primaryHighways.iterrows()]

print(sum(nVertexFlooded))
print(sum(nVertexHighways))
41787
2280
#a light version of the flooded areas
floodedAreasLight = floodedAreas.simplify(tolerance=0.0001)
nVertexFlooded = [np.array(row.exterior.coords).shape[0] for row in floodedAreasLight]
print(sum(nVertexFlooded))
3250
#create map
map = folium.Map(location = [refY,refX], tiles='OpenStreetMap', zoom_start = 12)

floodedJson = folium.GeoJson(data=floodedAreasLight.to_json(),
                             style_function=lambda x: {'fillColor': 'royalblue'})
primaryJson = folium.GeoJson(data=primaryHighways.to_json(),
                             style_function=lambda x: {'color': 'darkgreen'}) 
facilitiesJson = folium.GeoJson(data=primaryFacilities.to_json(),
                                style_function=lambda x: {'color': 'crimson'}) 

floodedJson.add_to(map)
primaryJson.add_to(map)
facilitiesJson.add_to(map)
map
Make this Notebook Trusted to load map: File -> Trust Notebook

Clip and export flooded schools

#analyze schools that are flooded
for index, facility in primaryFacilities.iterrows():
    floodedStatus = floodedDissolved.intersects(facility.geometry)[0]
    print(floodedStatus)
    primaryFacilities.loc[index,'Flooded']=floodedStatus
False
False
False
False
False
False
False
False
False
True
False
False
False
False
False
False
False
False
False
False
True
False
False
False
False
False
floodedFacilities = primaryFacilities[primaryFacilities['Flooded']==True]
floodedFacilities.head()

amenity name geometry Flooded
9 hospital Saint Alphonsus Eagle Health Plaza MULTIPOLYGON (((-116.34898 43.68815, -116.3488... True
20 hospital None MULTIPOLYGON (((-116.31779 43.68126, -116.3163... True
floodedFacilities.dtypes
amenity       object
name          object
geometry    geometry
Flooded       object
dtype: object
#export as geopackage
floodedFacilities['geometry'].to_file('../Out/floodedSchools.gpkg')

Input data

You can download the input data from this link.

Comment

Saul Montoya

Saul Montoya es Ingeniero Civil graduado de la Pontificia Universidad Católica del Perú en Lima con estudios de postgrado en Manejo e Ingeniería de Recursos Hídricos (Programa WAREM) de la Universidad de Stuttgart con mención en Ingeniería de Aguas Subterráneas y Hidroinformática.

 

Suscribe to our online newsletter

Subscribe for free newsletter, receive news, interesting facts and dates of our courses in water resources.