1. Introduction

This is an R Markdown Notebook which illustrates how to read a text file with geographic coordinates, and convert it into a spatial feature object with point geometries. Then, a shapefile of polygons is read. Finally, a leaflet map with polygons and points is created.

This notebook is written to help Geomatica Basica students at Universidad Nacional de Colombia to get started with geospatial data in R.

As previous attempts to publish the notebook from the ‘Knit to HTML’ RStudio option, failed several times due to the very big size of the html file (about 35 MB in this case), I decided to try simplifying the polygons.

This notebook illustrates the procedure to reduce the polygons size in order to get thin html files (about 1.5 MB in this case) and be able to publish leaflet maps in RPubs.

# The line below removes all variables from the current environment.
# run it from the terminal (NOT FROM HERE) everytime you start a new session
# rm(list=ls())
# install libraries from console -NOT FROM HERE -
# if (!require('devtools')) install.packages('devtools')
# devtools::install_github('rstudio/leaflet')
# load libraries
library(rgeos)
## Loading required package: sp
## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.7.2-CAPI-1.11.2 
##  Linking to sp version: 1.4-1 
##  Polygon checking: TRUE
library(sf)
## Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
## when loading 'dplyr'
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(leaflet)

2. Read table as dataframe

cities <- read.table("../datos/cities.txt", sep=";", quote="", header=TRUE)
# what is cities?
class(cities)
## [1] "data.frame"
# see the first records
head(cities)
##   X.id. X.cod. X.country.        X.name.    X.lat.    X.lon. X.alt.
## 1   "1"   2338  "Colombi"       "Bogota"  4.600000 -74.08334   2620
## 2   "2"   2339  "Colombi"         "Cali"  3.437222 -76.52250    758
## 3   "3"   2340  "Colombi"     "Medellin"  6.291389 -75.53611   2076
## 4   "4"   2341  "Colombi" "Barranquilla" 10.963889 -74.79639     33
## 5   "5"   2342  "Colombi"    "Cartagena" 10.399722 -75.51444     36
## 6   "6"   2343  "Colombi"       "Cucuta"  7.883333 -72.50528    314
colnames(cities) <- c("id", "cod", "country", "name", "lat", "lon", "alt")
# check column names
colnames(cities)
## [1] "id"      "cod"     "country" "name"    "lat"     "lon"     "alt"

3. Convert dataframe into a spatial feature

Let’s use the sf (simple features) library:

# convert dataframe object to simple feature object 
# https://stackoverflow.com/questions/44214487/create-sf-object-from-two-column-matrix
m <- st_as_sf(cities,coords = c(6,5,7))
# This another way of executing the above line
# cities %>%  sf::st_as_sf(coords = c(6,5,7)) -> m
# what is m?
class(m)
## [1] "sf"         "data.frame"
# get matrix of coordinates (lat, long, alt)
coords <- st_coordinates(m)
# 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]
alt = coords[,3]

4. Create leaflet map with departments and points

Use the sf library to read a shapefile of Colombian departments in 2018 –as provided by DANE–:

deptos <- sf::read_sf("../datos/shapefiles/MGN2018_DPTO_POLITICO/MGN_DPTO_POLITICO.shp")

What is the CRS carried by deptos?

st_crs(deptos)
## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

My major aim is to create a leaflet map representing departments as filled color polygons using the attribute Shape_Area. Then, I want to add the Colombian cities as markers.

After several failed attempts to publish the html file using the original geometry, I decided to simplify the polygons’geometry using the sf library functionalities.

In order to understand what is going on under the hood, you will benefit of understanding what are simple features.

First, let’s start by reprojecting the departments (otherwise, it will not be possible):

deptos2 <- deptos %>% st_transform(3116)

What is the CRS carried by deptos2?

st_crs(deptos2)
## Coordinate Reference System:
##   User input: EPSG:3116 
##   wkt:
## PROJCS["MAGNA-SIRGAS / Colombia Bogota zone",
##     GEOGCS["MAGNA-SIRGAS",
##         DATUM["Marco_Geocentrico_Nacional_de_Referencia",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6686"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4686"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",4.596200416666666],
##     PARAMETER["central_meridian",-74.07750791666666],
##     PARAMETER["scale_factor",1],
##     PARAMETER["false_easting",1000000],
##     PARAMETER["false_northing",1000000],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AUTHORITY["EPSG","3116"]]

Then, simplify the geometry:

deptos3  <- sf::st_simplify(deptos2, preserveTopology = TRUE, dTolerance = 1000)

Let’s check sizes of deptos and deptos3:

object.size(deptos)
## 16289960 bytes
object.size(deptos3)
## 178960 bytes

Let’s go back to the WGS84 CRS (this is the CRS used by leaflet):

deptos4 <- deptos3 %>% st_transform(4326)

Let’s apply a step-by-step procedure to create the leaflet map:

mapa <- leaflet(deptos4) 
mapa <- addTiles(mapa)
mapa <-   addPolygons(mapa, color = "#444444", weight = 1, smoothFactor = 0.5,
          opacity = 1.0, fillOpacity = 0.5,
          fillColor = ~colorQuantile("YlOrRd", Shape_Area)(Shape_Area),
          highlightOptions = highlightOptions(color = "red", weight = 2,
          bringToFront = TRUE))  
mapa <-  addMarkers(mapa, lng=long, lat=lat, popup=m$name)
mapa