No more GIS digitizing…

How boring, manual and prone to errors is GIS layers digitizing. Points of Interest (POI), roads or buildings are three very useful vector maps examples used in many applications such as map visualization, logistics, disaster management, marketing, sales and many others. However, creating these layers is a time consuming task for experts.

Fortunately, there are million of people who have made the job just for you. That is OSM, the biggest open source reservoir for GIS vector data in the world.

This article summarizes the way to extract poi, roads and buildings from OSM, and shows a python code example to easy perform the task.

Downloading Python Libraries

There are different libraries that help in the task, the main one is the Overpass API.

- Overpass

This api will connect directly to OSM, perform the query to download needed layers in a specific area of interest, defined by geographic coordinates. See details and other examples in https://python-overpy.readthedocs.io/

- Pandas and geopandas

Pandas is a very well known python library for data analysis, geopandas is the pandas library specific for geospatial data analysis. https://geopandas.org/en/stable/getting_started.html

- Shapely

It is used for the management of geometry of GIS layers, for instance creating point, line and polygons. https://shapely.readthedocs.io/en/stable/manual.html

- Json

This library allows to manage information in the widely used GeoJSON format. https://docs.python.org/3/library/json.html

- Folium

It makes it easy to visualize data that’s been manipulated in Python on an interactive leaflet map. It enables both the binding of data to a map for choropleth visualizations as well as passing rich vector/raster/HTML visualizations as markers on the map. https://python-visualization.github.io/folium/

Jupyter Notebooks or Google Collab

These are two options to run the python code, the very handy Jupyter Notebooks run locally, and the Google Collab Cyber Infrastructure.

For this example, Collab will be used. The first step is to connect personal Google Drive

from google.colab import drive
drive.mount('/content/drive', force_remount=True)
root_dir = "/content/drive/My Drive/"
base_dir = root_dir + 'myname_dir/'
print("base_dir is: ", base_dir)

Installing python libraries on collab via Pip

!pip install overpass
!pip install geopandas
!pip install shapely
!pip install folium

Now querying roads from a specific area of interest (a rectangular area in geographic coordinates WGS84). Pois and buildings can be obtained in a similar fashion.

# This code downloads roads for a AOI defined by: Xmin, Ymin, Xmax, Ymax

import os
import overpass
import geopandas as gpd
import time
import json

myname = "roads_Retiro.geojson"
mynameout = "roads_Retiro.shp"
#myname = "roads_UIUC.geojson"
#mynameout = "roads_UIUC.shp"

map_extent = 'way["highway"](6.053347, -75.509213, 6.067174, -75.495738);'
#map_extent = 'way["highway"](40.102735, -88.234170, 40.109743, -88.224600);'

mytime = time.time()
api = overpass.API()
data = api.get(map_extent, verbosity='geom', responseformat="geojson")
[f for f in data.features  if f.geometry['type'] == "LineString"]

with open(os.path.join(base_dir, myname), 'w', encoding='utf-8') as myfile:
    json.dump(data, myfile, ensure_ascii=False, indent=4)

print("Roads have been obtained from OSM: ", myname, 'in', base_dir)

gdf = gpd.read_file(os.path.join(base_dir, myname))
gdf.to_file(os.path.join(base_dir, mynameout), driver='Shapefile')

mytimef = time.time()
print("-> GeoJSON Roads are converted to shapefile sucessfully:", 
round((mytimef-mytime),2), 'seconds')

The AOI polygon can be created in a similar way.

from shapely.geometry import Polygon

myAOI = base_dir + 'AOI.geojson'

# lat_points = [ymax, ymin, ymin, ymax]    
lat_points = [6.0584, 6.0547, 6.0547, 6.0584]

# lon_points = [xmin, xmin, xmax, xmax]
lon_points = [-75.5028, -75.5028, -75.4964, -75.4964]

polygon_geom = Polygon(zip(lon_points, lat_points))
crs = {'init': 'epsg:4326'}
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])       
print(polygon.geometry)

#polygon.to_file(myAOI, driver="ESRI Shapefile")
polygon.to_file(myAOI, driver='GeoJSON')

print('xmin,ymin', lon_points[0], lat_points[1])
print('xmax,ymax', lon_points[2], lat_points[0])

Now, road vector layers obtained from OSM within a specific AOI can be visualized using an interactive leaflet map.

myroads = base_dir + myname
#mypath = base_dir + 'roads_UIUC.geojson'

mymap = folium.Map([6.0534443537329, -75.5057858485094],zoom_start=15)
#mymap = folium.Map([40.10273, -88.234170],zoom_start=15)

folium.GeoJson(myroads).add_to(mymap)
folium.GeoJson(myAOI).add_to(mymap)

# Adding pinpoint coordinates to map
folium.LatLngPopup().add_to(mymap)

# save the map in a html file
## myhtml = base_dir + 'roads.html'
## mymap.save(myhtml)

# open the previous html file in a browser
### import webbrowser
### webbrowser.open(myhtml)

mymap

Visualizing Roads obtained from OSM and exported AOI polygon from coordinates

Revisiting Code

• Change the word “way” by “poi” or “building” to obtain corresponding vector layers from OSM in the specific AOI.
• Uncomment code to obtain roads around the University of Illinois at Urbana-Champaing (UIUC).
• Uncomment the double “##” in the folium code and see the resultant separate html file.
• Also, if html file needs to be opened right after saved, uncomment triple “###”.

References

• A GIS Pipeline to Produce GeoAI Datasets from Drone Overhead Imagery https://www.mdpi.com/2220-9964/11/10/508

Support me

Enjoying my work? Show your support with Buy me a coffee, a simple way for you to encourage me and others to write. If you feel like it, just click the next link and I will enjoy a cup of coffee!