Abstract¶
This set of slides is an introduction to geospatial data visualization with roadkill data. It will show how to construct geospatial objects from tabular vector data or shapefiles. To work with these geospatial data in Python, we can use GeoPandas, shapely or related libraries for manipulating and analyzing the data. This slideset will also cover the following parts: loading data from CSV and shapefile, visualizing geospatial data with matplotlib and folium, using spatial operations to make thematic maps, fetching the street newtorks from OpenStreetMap, as well as visualizing them with the osmnx package. All of our datasets come from the research data repository data.depositar.io except street networks data. The techniques presented here can be used for advanced spatial analysis and for exploring spatial relation in various application scenarios.
1. Loading a sample roadkill dataset¶
The Pandas library already has a built-in function read_csv( ) for retrieving data from a public link. We can use it to download the roadkill dataset and view the first 5 items.
rk_species = pd.read_csv('https://data.depositar.io/dataset/b78e1f21-1f5d-4fb6-95b5-619f3cb44797/resource/0158ed0f-1989-4a64-b67c-e7962c353d57/download/export-example.csv')
rk_species.head()
Converting an existing DataFrame into an GeoDataFrame¶
rk_species['紀錄時間'] = pd.to_datetime(rk_species['紀錄時間'])
rk_geo = [gpd.geoseries.Point(loc) for loc in zip(rk_species['經度'], rk_species['緯度'])]
rk_species = gpd.GeoDataFrame(rk_species, geometry=rk_geo, crs={'init': 'epsg:4326'})
rk_species.geometry.head()
rk_species.plot(marker='.', color='b')
Using the Folium library to visualize our roadkill dataset¶
Each observation can be represented as a point on the map. It need to be converted to an appropriate format for the HeatMap plugins in Folium.
osmMap = folium.Map(location=[23.5, 121], zoom_start=7)
osmMap.add_child(plugins.HeatMap([(geo.y, geo.x) for geo in rk_species.geometry],
radius=15, blur=12))
Spliting the dataset into specific group depending on the columns we choose¶
Using the groupby function to split data depending on the scientific names of species and aggregating these data to find out which one is the most observed in our dataset.
rk_species.groupby(['物種名', '物種名（中文）']).size().sort_values(ascending=False).head()
target = rk_species[rk_species['物種名'] == 'Duttaphrynus melanostictus']
targetMap = folium.Map(location=[24, 121], zoom_start=7)
targetMap.add_child(plugins.HeatMap([(geo.y, geo.x) for geo in target.geometry],
radius=15, blur=12))
Checking out the number of observations by months¶
In this part, the column of record time is set to be the index of DataFrame. It help us extract timeseries data from our dataset for different date ranges.
rk_species.index = rk_species['紀錄時間']
rk_species.head()
timeseries_data = rk_species.resample('M').size()
timeseries_data
fig, ax = plt.subplots()
ax.plot_date(timeseries_data.index, timeseries_data, linestyle='-')
ax.set_ylabel('Frequency')
2. Loading the Tainan City (台南市) shapefile¶
link = 'https://data.depositar.io/dataset/796293f0-bee5-4810-918b-5a7df49d958b/resource/be3942ce-b098-49e2-8d54-07258c48b78b/download/2_tainan_3.zip'
def loadshp(link):
zipfile = ZipFile(BytesIO(urlopen(link).read()))
ext = [x for x in zipfile.namelist() if x.endswith(('shp', 'dbf', 'shx', 'prj')) and not x.startswith('__')]
filename = {e.split('.')[1]:e for e in ext}
return [BytesIO(zipfile.read(v)) for k, v in filename.items()]
dbf, prj, shp, shx = loadshp(link)
county = shapefile.Reader(shp=shp, shx=shx, dbf=dbf)
county.numRecords
attributes, geometry = [], []
field_names = [cln[0] for cln in county.fields[1:]]
for r in county.shapeRecords():
geometry.append(shape(r.shape.__geo_interface__))
r.record = map(lambda x: x.decode('big5') if isinstance(x, bytes) else x, r.record)
attributes.append(dict(zip(field_names, r.record)))
attributes[1]
x, y = geometry[1].exterior.xy
plt.figure(figsize=(3,3))
plt.plot(x, y)
gdf = gpd.GeoDataFrame(data=attributes, geometry=geometry, crs={'init': 'epsg:4326'})
gdf.head()
gdf.plot(color='w', edgecolor='b')
Creating a GeoPandas DataFrame from the list of geometry, attributes and setting the CRS associated with the geometry column in GeoPandas object. This make us can be enabled the spatial operation on these objects.
Coordinate reference system (CRS)¶
CRS refers to Coordinate Reference System. You can find the epsg codes for various projections from epsg.io . If your geospatial objects come from different resources with differernt CRS, you need to convert to the same CRS before you apply spaital operators to them.
gdf.head(2)
gdf.to_crs(epsg=3826).head(2)
3. Spatial Join by Geographic Region or Attributes¶
count = gpd.sjoin(rk_species, gdf, op='within').groupby(['townname']).size().to_frame()
count.columns = ['count']
count = count.reset_index()
count.head()
gdf = gdf.merge(count, on='townname', how='left')
gdf = gdf.fillna(0)
gdf.head()
gdf.plot(column='count', figsize=(15, 10), cmap=plt.cm.Blues, legend=True)
4. Fetch OSM Networks from Points in Tainan¶
OpenStreetMap (OSM) is a successful Volunteered Geographic Information (VGI) system. The data provided by the OSM is crowdsourced, completely free to use, can be redistributed to anyone for any purpose. So it's very helpful to create custom maps in Jupyter Notebooks.
city = ox.graph_from_point((gdf.unary_union.centroid.y, gdf.unary_union.centroid.x),
distance=10000,
network_type='drive')
ox.plot_graph(city, fig_height=10, fig_width=10, edge_color='black')
nodes_prj, edges_prj = ox.graph_to_gdfs(city, nodes=True, edges=True)
boundary = box(*edges_prj.unary_union.bounds)
print(boundary)
fig, ax = plt.subplots()
boundarySelection = gpd.GeoDataFrame(geometry=[boundary], crs=rk_species.crs)
rk_species.plot(ax=ax, markersize=5, color='b')
boundarySelection.plot(ax=ax, color='white', edgecolor='r')
rk_species_filter = gpd.sjoin(rk_species, boundarySelection, op='within')
fig, ax = plt.subplots()
edges_prj.plot(ax=ax)
rk_species_filter.plot(ax=ax, markersize=10, color='red')
print(edges_prj.geometry[0])
partcity = ox.graph_from_point((23.225, 120.375), distance=5000, network_type='drive')
partmap = ox.plot_graph_folium(partcity)
partmap
fig, ax = plt.subplots()
edges_prj.plot(ax=ax)
rk_species_filter.buffer(0.0005).plot(ax=ax, color='r') # 50m
rk_species_filter_buf = rk_species_filter.buffer(0.0005)
edges_prj['rkcount'] = 0
addList = []
for p in rk_species_filter_buf:
for r in edges_prj[edges_prj.intersects(p)].iterrows():
addList.append(r[0])
counts = pd.DataFrame(addList)[0].value_counts()
for i in range(len(counts.index)):
edges_prj.at[counts.index[i], 'rkcount'] = counts[counts.index[i]]
edges_prj.plot(column='rkcount', cmap='YlOrRd')
edges_prj[edges_prj['rkcount']>10]