1. Introduction

This is an R Markdown Notebook which illustrates how to use geospatial web services from R. In particular, it shows how to scrap geographic features from any web feature service (WFS) hosted by ArcGIS Server. This notebook aims to help Geomatica Basica students at Universidad Nacional de Colombia to get started with Rgis.

We will the esri2sf package which enables users to get vector data delivered by an ArcGIS Server through the server’s REST API. The package allows users to download geographic features from the server as simple features (sf). In turn, the sf and geospatial libraries offer several functionalities for spatial analysis and thematic mapping.

2. Preliminaries

Install packages (if this is your first time with this notebook):

  knitr::opts_chunk$set(echo = TRUE) 
  ##knitr::read_chunk("MyExternalRCode.R") 
#library(devtools)
#install_github("yonghah/esri2sf")

Clean the workspace:

rm(list=ls())

Load packages:

library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(leaflet)
library(leaflet.extras)
library(htmltools)
library(esri2sf)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

3. Getting simple features from the web

How it works

The esri2sf program sends a request to an ArcGIS Server and gets json responses containing coordinates of geometries of which format is not the same as geojson. So it converts the json into simple feature geometries from the response. Then it combines attribute data to the geometries to create sf dataframe. Often ArcGIS servers limits the maximum number of rows in the result set. So this program creates 500 features per request and automatically re-send requests until it gets all features in the dataset.

Point data

We can get the Punto Geodesico layer from the Mapa de Referencia provided by IDECA. By the way, what is IDECA? Which web services provide IDECA? You can find infomation here.

Let’s know what is a punto geodesico:

Let’s get the puntos geodesicos from IDECA WFS:

library("esri2sf")
url5 <- "http://oaiee.scj.gov.co/agc/rest/services/Ideca/Mapa_Referencia_Ideca/FeatureServer/5"
pg <- esri2sf(url5)
## [1] "Feature Layer"
## [1] "esriGeometryPoint"

What is pg? A simple feature collection with POINT geometry type. How many puntos geodesicos are located in Bogota? Let’s see:

pg

Let’s visualize the puntos geodesicos:

# get matrix of coordinates (lat, long, alt)
coords <- st_coordinates(pg)
# create vector of coordinates
# what are vectors in R
# https://datascienceplus.com/vectors-and-functions-in-r/#:~:text=A%20vector%20is%20the%20simplest,a%20vector%20of%20logical%20values.
lat = coords[, 2]
long = coords[,1]
mapa1 <- leaflet() 
mapa1 <-   addTiles(mapa1)  # Add default OpenStreetMap map tiles
mapa1 <-   addMarkers(mapa1, lng=long, lat=lat, popup=pg$PGENOMBRE)
# show map
mapa1

We can save the data for later usage:

st_write(obj = pg, dsn = "./datos/ptosgeodesicos.gpkg", append=FALSE)
## Deleting layer `ptosgeodesicos' using driver `GPKG'
## Writing layer `ptosgeodesicos' to data source `./datos/ptosgeodesicos.gpkg' using driver `GPKG'
## Writing 48 features with 10 fields and geometry type Point.

Polyline data

We can get polyline data such as the Bogota’s malla vial. Let’s get the malla vial from IDECA WFS. Be patient that it takes several minutes to download about 137000 features.

url14 <- "http://oaiee.scj.gov.co/agc/rest/services/Ideca/Mapa_Referencia_Ideca/FeatureServer/14"
mv <- esri2sf(url14)
## [1] "Feature Layer"
## [1] "esriGeometryPolyline"

What is mv? A simple feature collection with MULTILINESTRING geometry type.

mv

Let’s select only the major streets, that is, those classified as MVITCLA = 1.

avenidas = filter(mv, MVITCLA == 1)
avenidas

Let’s visualize the major streets:

mapa2 <- leaflet()
mapa2 <- addTiles(mapa2)
mapa2 <- addPolylines(mapa2, data = avenidas)
mapa2

We can save the data for later usage:

st_write(obj = avenidas, dsn = "./datos/avenidas.gpkg", append=FALSE)
## Deleting layer `avenidas' using driver `GPKG'
## Writing layer `avenidas' to data source `./datos/avenidas.gpkg' using driver `GPKG'
## Writing 17187 features with 13 fields and geometry type Multi Line String.

Polygon data

We can get polygon data such as the Bogota’s localidades.

Let’s get the localidades from IDECA WFS:

url47 <- "http://oaiee.scj.gov.co/agc/rest/services/Ideca/Mapa_Referencia_Ideca/FeatureServer/47"
loc <- esri2sf(url47)
## [1] "Feature Layer"
## [1] "esriGeometryPolygon"

What is loc? What is the attribute used to store area of each localidad?

loc
loc$area <- loc$LOCAREA/1000000
loc

Let’s create a palette to visualize the localidades We will use the area attribute to define the polygons’ fill color:

bins <- c(0, 20, 50, 100, 200, 300, 400, 500, 1000, Inf)
pal <- colorBin("YlOrRd", domain = loc$area, bins = bins)

Let’s create labels:

labels <- sprintf("<strong>%s</strong><br/>%g km<sup>2</sup>", loc$LOCNOMBRE, loc$area) %>% lapply(htmltools::HTML)

Now, it’s mapping time:

mapa3 <- leaflet(loc) %>% addProviderTiles(providers$OpenStreetMap) %>%
                    addPolygons(
                    fillColor = ~pal(area),
                    weight = 2,
                    opacity = 1,
                    color = "white",
                    dashArray = "3",
                    fillOpacity = 0.7,
                    highlight = highlightOptions(
                   weight = 5,
                   color = "#666",
                   dashArray = "",
                  fillOpacity = 0.7,
                   bringToFront = TRUE), 
                    ##label= ~LOCNOMBRE,
                    ##label = labels,
                    ##labelOptions = labelOptions(
                    ##style = list('font-weight' = 'normal', padding = '3px 8px'),
                    ##textsize = '15px',
                    ##direction = 'auto'
                    popup = paste0("Localidad: ", loc$LOCNOMBRE, "<br>",
                                   "Area: ", loc$area)) %>%
                   addLegend(pal = pal, values = ~area, opacity = 0.7, title = "Area [km2]", position = "bottomright")

mapa3

Note that, albeit the mapa3 displays perfectly as a local html, it refuses to get published at RPubs. Thus, I decided not to evaluate the previous code chunk and to include a static map.

We can save the data for later usage:

st_write(obj = loc, dsn = "./datos/localidades.gpkg", append=FALSE)
## Deleting layer `localidades' using driver `GPKG'
## Writing layer `localidades' to data source `./datos/localidades.gpkg' using driver `GPKG'
## Writing 20 features with 6 fields and geometry type Multi Polygon.

Now it’s your turn to try these functionalities. Good luck!!!

4. Reproducibility

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_1.0.2          esri2sf_0.1.1        htmltools_0.5.0     
## [4] leaflet.extras_1.0.0 leaflet_2.0.3.9000   sf_0.9-6            
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5         RColorBrewer_1.1-2 compiler_3.6.3     pillar_1.4.6      
##  [5] class_7.3-17       tools_3.6.3        digest_0.6.25      jsonlite_1.7.1    
##  [9] evaluate_0.14      lifecycle_0.2.0    tibble_3.0.4       pkgconfig_2.0.3   
## [13] rlang_0.4.8        DBI_1.1.0          crosstalk_1.1.0.1  curl_4.3          
## [17] yaml_2.2.1         xfun_0.18          e1071_1.7-3        stringr_1.4.0     
## [21] httr_1.4.2         knitr_1.30         generics_0.0.2     vctrs_0.3.4       
## [25] htmlwidgets_1.5.1  classInt_0.4-3     grid_3.6.3         tidyselect_1.1.0  
## [29] glue_1.4.2         R6_2.4.1           rmarkdown_2.3.5    farver_2.0.3      
## [33] purrr_0.3.4        magrittr_1.5       scales_1.1.1       ellipsis_0.3.1    
## [37] units_0.6-7        colorspace_1.4-1   KernSmooth_2.23-17 stringi_1.5.3     
## [41] munsell_0.5.0      crayon_1.3.4