In dimitrisk/goal: GeOgraphical Analysis Library

knitr::opts_chunk$set(warning = FALSE, message = FALSE) 

Introduction

Use of goal package and more specificaly the 'osm' group of functions for the analysis of OpenStreetMap POIs (Points of Interest) and road network. The following two links, can help on identifying OSM entities (Municipalities, cities, ...)

• Search OSM

• Bounding BOX

POIs analysis

Aim: Analyse the location of Amenities and Shops in this municipality. First install goal library

library(devtools)
#devtools::install_github("dimitrisk/goal", quiet=T)
library(goal)
#library(osmdata)
library(raster)

Amenities

Download all amenities in this municipality (Points or Polygons or Multipolygons)

am = osm.getPOI(inPerioxi = "Mytilene Municipal Unit", inkey = "amenity")
am

Alternatively, we can use a bounding box of an area instead of the name of a municipality.

myt = osm.getPOI_usingbb( c(26.547303,39.101658,26.564641,39.113247), inkey ="amenity" )
myt

After downloading the POIs (am), we merge all POIs together (Points, Polygons, Multipolygons), and convert them all to points.

amenities = osm.combineAmenities(am)
amenities

Frequency table of amenities types in this municipality

freq1 = osm.getFrequency(amenities, inword = "amenity", removeNA = F)
freq1

Shops

Get all shops in this municipality (Points, Polygons, Multipolygons)

sh = osm.getPOI(inPerioxi = "Mytilene Municipal Unit", inkey = "shop")
sh

#' Merge all shops (Points, Polygons, Multipolygons), and convert them to points.
shops = osm.combineShops(sh)
shops

Frequency table of shop types in this municipality

freq2 = osm.getFrequency(shops, inword = "shop", removeNA = F)
freq2

Density

Create empty raster of this area and get the polygon of this area.

r = osm.CreateEmptyRaster(inPerioxi = "Mytilene Municipal Unit")

#' Get polygon of the area
mypol = osmdata::getbb("Mytilene Municipal Unit", format_out = "sf_polygon")$multipolygon

Sum of shops per cell

density1 = raster::rasterize(shops, r, field=1, fun=sum) # sum
plot(density1)
plot(mypol, add=T)

Sum of amenities per cell

density2 = raster::rasterize(amenities, r, field=1, fun=sum) # sum
plot(density2)
plot(mypol, add=T)

Presense/Absense

Presense/Absence of shops in cells

presense1 = raster::rasterize(shops, r, field=1) # presense
plot(presense1)
plot(mypol, add=T)

Presense/Absence of shops in cells

presense2 = raster::rasterize(amenities, r, field=1) # presense
plot(presense2)
plot(mypol, add=T)

Web Map

#' web map1
library(leaflet)
leaflet(amenities) %>% addTiles() %>%  leaflet::addMarkers(popup = ~as.character(amenity), label = ~as.character(amenity))

#' web map2
library(simplevis)
library(dplyr)
library(palmerpenguins)
library(stars)
simplevis::leaf_sf_col(amenities,  col_var = amenity)

Custom bounding box

We now use a custom bounding box and then download data only for this area.

am = osm.getPOI_usingbb( c(26.547303,39.101658,26.564641,39.113247), inkey ="amenity" ) # Download
mypol_bb = osm.osmdata_result_2_bbox_pol(am) # get polygon of this bounding box.
amenities = osm.combineAmenities(am) # combine all amenities into points.

Now get the general polygon of this municipality (coastline)

mypol2 = osmdata::getbb("Mytilene Municipal Unit", format_out = "sf_polygon")$multipolygon

q=c(26.545029,39.088569,26.570177,39.116810)
poly = goal::osm.bb_2_pol(inVec=q, outcrs=2100) 
#plot(poly)

What we have so far:

plot(mypol2, border="red")
plot(mypol_bb, add=T)
plot(amenities, add=T)

zoomed in:

plot(mypol_bb )
plot(mypol2, border="red", add=T)
plot(amenities, add=T)

We need to calculate the intersection between municipality polygon (coastline) and the custom square bounding box we have defined.

clip1 = sf::st_intersection(mypol_bb,mypol2)

show the area

plot(clip1)
plot(mypol2, border="red", add=T)
plot(amenities, add=T)

Create an empty raster of this intersection:

r = raster(clip1 )

Calculate the density of amenities (points) in this clipped area:

density1 = raster::rasterize(amenities, r, field=1, fun=sum) # sum
plot(density1)
plot(clip1, add=T)
plot(as_Spatial(amenities), add=T, cex=0.5)

Road network analysis

Get all roads of this municipality by using just the municipality name:

myr1 = osm.getRoads(inPerioxi = "Mytilene Municipal Unit", 
                    withBB=FALSE, outcrs = 2100)
myr1
plot(myr1, col="blue", main="sfnetworks Roads (by municipality name)")

Get all roads of this area by using a bounding box:

myr2 = osm.getRoads(incoordinates = c(26.545029,39.088569,26.570177,39.116810), 
                    withBB=TRUE, outcrs = 2100 )
myr2
plot(myr2, col="grey", main="sfnetworks Roads (by bbox)")

example

q=c(26.545029,39.088569,26.570177,39.116810)
net2 = osm.getRoads(q, withBB=TRUE, outcrs=4326)
poly = osm.bb_2_pol(q, outcrs =  4326)

net3 = osm.ClipSFnetwork_with_poly(net2, poly)

plot(net3,col="grey", main="Clipped sfnetwork")
plot(poly,add=T)

Walkability

Walkability is a measure of accessibility of amenities by foot. A possible proxy way of measuring walkability in a city can be considered the total length of footways of the city. The following function, calculates the total length (and some other statistics) of footways in a city:

library(goal)
mylength = osm.getLength_footway( place="Mytilene Municipal Unit" )
mylength

dimitrisk/goal documentation built on April 30, 2024, 7:21 p.m.