Paris: Cut and Load OSM Polygon - plot with TMAP
Jesper Hybel
18/10/2020
Here I show how to make a spatial selection from OSM and load the data into R.
The method I use is described by James Chevalier here:
https://github.com/jameschevalier/cities
The first step is to
- Get The Poly File For A Specific City
There are two ways to do this:
(1.1) Go here https://extract.bbbike.org/ and make a polygon and press extract to have it sent to your mail.
(1.2) Or go here http://nominatim.openstreetmap.org/ and enter the name of the area. Then click on “details” and get the OSM-relation number. With the OSM relation number go here http://polygons.openstreetmap.fr/. Enter the relation number and select one of the poly files. The browser opens and shows numbers defining the polygon. Right click in the browser and choose save. Alter the extension of the file to .poly. You then need to do (1.2.1)-(1.2.3) which is installation of osmosis.
(1.2.1) Install osmosis and add it to your path variable.
(1.2.2) Download the region as .pbf file and containing the spatial area you wanted. So if the spatial area is Paris download for example France as .pbf file.
(1.2.3) Finally execute osmosis command in the directory where the .pdf file is
osmosis --read-pbf-fast file="YOUR-REGION-latest.osm.pbf" --bounding-polygon file="CITY-NAME_STATE.poly" --write-xml file="CITY-NAME_STATE.osm"
For example:
osmosis --read-pbf-fast file="france-latest.osm.pbf" --bounding-polygon file="paris.poly" --write-xml file="holyoke_ma.osm"this should then generate an extract file .pbf in the directory and this can then be used.
The .pbf file can be converted using the package OSMEXTRACT. The package is not CRAN so install by
library(osmextract)
mypath <- "C://users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//test"
mydata <- oe_get(place="Denmark",download_directory=mypath)and then load
library(osmextract)## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Any product made from OpenStreetMap must cite OSM as the data source.
## Geofabrik data are taken from https://download.geofabrik.de/
## For usage details of bbbike data see https://download.bbbike.org/osm/
## OpenStreetMap_fr data are taken from http://download.openstreetmap.fr/to read the .pbf file simply use
mypath <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//paris"
setwd(mypath)
mydata_line <- oe_read("paris.osm.pbf")
mydata_line## Simple feature collection with 414779 features and 9 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 2.008637 ymin: 48.63208 xmax: 2.663 ymax: 49.08814
## geographic CRS: WGS 84
## First 10 features:
## osm_id name highway waterway aerialway barrier
## 1 2569 Rue Nationale tertiary <NA> <NA> <NA>
## 2 2570 Avenue des Gobelins secondary <NA> <NA> <NA>
## 3 2573 Place des Alpes residential <NA> <NA> <NA>
## 4 2574 Rue Taine tertiary <NA> <NA> <NA>
## 5 286314 Rue Saint-Louis en l'Île residential <NA> <NA> <NA>
## 6 286537 Quai de Bourbon living_street <NA> <NA> <NA>
## 7 286538 Rue Jean du Bellay living_street <NA> <NA> <NA>
## 8 286539 Quai de Béthune living_street <NA> <NA> <NA>
## 9 286540 Rue Saint-Louis en l'Île living_street <NA> <NA> <NA>
## 10 286541 Rue Saint-Louis en l'Île living_street <NA> <NA> <NA>
## man_made z_order
## 1 <NA> 4
## 2 <NA> 6
## 3 <NA> 3
## 4 <NA> 4
## 5 <NA> 3
## 6 <NA> 0
## 7 <NA> 0
## 8 <NA> 0
## 9 <NA> 0
## 10 <NA> 0
## other_tags
## 1 "oneway"=>"yes","surface"=>"asphalt","maxspeed"=>"50"
## 2 "lit"=>"yes","oneway"=>"yes","surface"=>"paving_stones","maxspeed"=>"50","wikidata"=>"Q2874102","source:maxspeed"=>"FR:urban"
## 3 "oneway"=>"yes","surface"=>"asphalt","maxspeed"=>"30","wikidata"=>"Q3390483"
## 4 "lit"=>"yes","surface"=>"asphalt","maxspeed"=>"50","wikidata"=>"Q3450115","busway:right"=>"lane","cycleway:right"=>"share_busway"
## 5 "lit"=>"yes","lanes"=>"2","surface"=>"asphalt","maxspeed"=>"50","source:maxspeed"=>"sign"
## 6 "lanes"=>"1","oneway"=>"yes","surface"=>"asphalt","cycleway"=>"opposite","maxspeed"=>"20","wikidata"=>"Q778365","zone:maxspeed"=>"FR:20"
## 7 "lanes"=>"1","oneway"=>"yes","bicycle"=>"yes","surface"=>"sett","cycleway"=>"opposite","maxspeed"=>"20","zone:maxspeed"=>"FR:20"
## 8 "lanes"=>"1","oneway"=>"yes","surface"=>"asphalt","cycleway"=>"opposite","maxspeed"=>"20","wikidata"=>"Q3412680","zone:maxspeed"=>"FR:20"
## 9 "lit"=>"yes","lanes"=>"1","oneway"=>"yes","surface"=>"asphalt","cycleway"=>"opposite","maxspeed"=>"20","wikidata"=>"Q3449922","zone:maxspeed"=>"FR:20"
## 10 "lanes"=>"1","oneway"=>"yes","surface"=>"asphalt","cycleway"=>"opposite","maxspeed"=>"20","wikidata"=>"Q3449922","zone:maxspeed"=>"FR:20"
## geometry
## 1 LINESTRING (2.369443 48.821...
## 2 LINESTRING (2.355171 48.832...
## 3 LINESTRING (2.358575 48.832...
## 4 LINESTRING (2.39454 48.8394...
## 5 LINESTRING (2.360093 48.850...
## 6 LINESTRING (2.357066 48.852...
## 7 LINESTRING (2.353301 48.852...
## 8 LINESTRING (2.355918 48.851...
## 9 LINESTRING (2.356547 48.851...
## 10 LINESTRING (2.356547 48.851...which should create a new file in the mypath directory called paris.gpkg.
By default oe_read() gets the lines from “paris.osm.pbf” and some of these lines are the streets.
Streets are stored under the key called “highway” and to get a view of the particular values make a table of the highway variable.
table(mydata_line$highway)##
## * abandoned access
## 1 3 1
## access_ramp bridleway bus
## 5 384 1
## bus_guideway bus_stop construction
## 117 1 422
## corridor cycleway disused
## 203 5880 10
## disused:footway elevator emergency_access_point
## 1 95 1
## footpath footway lane
## 1 71307 5
## living_street motorway motorway_link
## 2378 1612 1458
## parh path pedestrian
## 1 16870 3784
## platform primary primary_link
## 13 10537 1283
## proposed raceway residential
## 31 6 67649
## road secondary secondary_link
## 124 13070 537
## service services steps
## 68573 10 12227
## tertiary tertiary_link track
## 12373 288 6920
## traffic_island trunk trunk_link
## 1 1142 1097
## unclassified virtual
## 10587 12The road-network can then be plotted using the package TMAP. This is done here:
library(tmap)
ht = c("primary","secondary","primary_link","secondary_link")
index <- which(mydata_line$highway%in%ht)
tm_shape(mydata_line[index,"highway"])+tm_lines()## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
setwd(mypath)
mydata_poly <- oe_read("paris.osm.pbf",layer="multipolygons")index1 <- which(mydata_poly$admin_level==7)
tmap_mode("plot")## tmap mode set to plottingtm_shape(mydata_poly[index1,"admin_level"]) + tm_fill("seashell3")+tm_borders("white", lwd = 3) +
tm_shape(mydata_line[index,"highway"])+ tm_lines(col="grey29") 