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))

#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.

Blog

Introduction to Python and Geopandas for Flooded Area Analysis - Tutorial