Reading, filtering and visualizing geospatial vector data in R
Write your name here
Write the publishing date here
1. Introduction
This is an R Markdown Notebook which illustrates basic functions provided by the simple features (sf) and the tidyverse libraries. It aims to diffuse R geospatial functionalities for students of Geomatica Basica at Universidad Nacional de Colombia.
The first step is to install the required libraries:
#install.packages(c("tidyverse", "sf"))Then, we need to load the two libraries:
library(tidyverse)## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.0 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.22. Reading vector data
Now, we can read a shapefile representing Colombian departments which has been previously downloaded from DIVA-GIS.
deptos <- read_sf("./colombia/COL_adm1.shp")We can know what is the contents of the deptos object:
head(deptos)It’s very useful to review basic functions provided by the sf library. For example, to know what is the coordinate reference system of the vector data stored in the deptos object we can use the st_crs function:
st_crs(deptos)## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"3. Using ggplot to visualize geospatial data:
A good starting point to understand data is to plot the data. We can use the ggplot functionalities:
ggplot() + geom_sf(data = deptos) 
It is possible to use any coordinate reference system to plot the data. However, it is not right to use a coordinate reference system which has been explicitly defined for another country or region. More information on EPSG codes can be found here
# The CRS 3978 is used in Canada
ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3978))
Look for the properties of the CRS with EPSG code 32618. It corresponds to UTM 18 N. In case we need to use such CRS,it is necessary to convert the spatial object from EPSG4326 into EPSG:32618.
deptos_utm <- st_transform(deptos, crs = st_crs(32618))deptos_utmggplot() + geom_sf(data = deptos_utm)
4. Filtering geospatial data based on attributes
As we are interested only in one departament, we can filter the data. Review the following chunk of code and change it to match your department:
antioquia <- deptos %>% filter(NAME_1 == "Antioquia")Let’s plot the new object:
ggplot() + geom_sf(data = antioquia) 
We can repeat the previous steps to load the Colombian municipalities and filter the Antioquian ones:
munic <- read_sf("./colombia/COL_adm2.shp")
mun_antioquia <- munic %>% filter(NAME_1 == "Antioquia")
ggplot() + geom_sf(data = mun_antioquia) 
mun_antioquiaantioquia_points<- st_centroid(mun_antioquia)## Warning in st_centroid.sf(mun_antioquia): st_centroid assumes attributes are
## constant over geometries of x## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude dataantioquia_points <- cbind(mun_antioquia, st_coordinates(st_centroid(mun_antioquia$geometry)))## Warning in st_centroid.sfc(mun_antioquia$geometry): st_centroid does not give
## correct centroids for longitude/latitude dataA better plot can be produced using the following chunk:
ggplot(antioquia) +
geom_sf() +
geom_sf(data = antioquia_points, fill = "antiquewhite") +
geom_text(data = antioquia_points, aes(x=X, y=Y,label = ID_2), size = 2) +
coord_sf(xlim = c(-78, -73.5), ylim = c(5.3, 9), expand = FALSE)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorggplot(antioquia) +
geom_sf(data=antioquia_points, aes(x=X, y=Y, fill =
ID_2), color = "black", size = 0.25) +
geom_text(data = antioquia_points, aes(x=X, y=Y,label = ID_2), size = 2) +
theme(aspect.ratio=1)+
scale_fill_distiller(name="ID_2", palette = "YlGn", breaks = pretty_breaks(n = 5))+
labs(title="Another Map of Antioquia")## Warning: Ignoring unknown aesthetics: x, y
However, this visualization is not a real map. Anyway, the output can be saved either as a pdf or as a png.
ggsave("antioquia_municipios.pdf")## Saving 7 x 5 in imageggsave("map_antioquia.png", width = 6, height = 6, dpi = "screen")5. Using leaflet to visualize data
First, let’s install the required library:
#install.packages("leaflet")Now, load the library:
library(leaflet)In order to use the library, we need to convert from simple features to spatial points.
ant_points <- as(antioquia_points, 'Spatial')Uncomment the following command to see what is inside the object:
#head(ant_points)Get area of municipalities:
install.packages("lwgeom")## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)Load a library:
library(lwgeom)## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.5.1, PROJ 4.9.2Let’s calculate the area of each municipality (in square meters):
mun_antioquia$area <- st_area(mun_antioquia) #Take care of unitsNow, let’s create a new field to store area in square kilometers:
mun_antioquia$km2 <- mun_antioquia$area/(1000000)Check the output:
mun_antioquia$km2## Units: [m^2]
## [1] 507.60408 281.74659 131.99541 102.16135 1181.13943 540.10263
## [7] 87.93221 364.57961 1412.77553 222.17398 573.79307 482.78641
## [13] 233.75043 130.72935 221.36748 150.38938 309.02143 225.72397
## [19] 255.26451 409.74933 358.61191 308.97097 2034.38249 465.21183
## [25] 222.97713 152.67026 212.66511 309.89046 78.87443 246.30879
## [31] 170.91354 1320.49864 673.71393 39.00079 262.90943 155.60562
## [37] 201.44867 66.26547 1831.29644 183.02783 220.88909 1927.63127
## [43] 423.67064 64.55167 192.70603 46.77686 225.71517 1315.08305
## [49] 387.27178 101.54419 82.68621 220.66230 82.19749 159.73225
## [55] 81.43762 113.26362 45.62970 23.47969 2129.96167 237.63511
## [61] 177.93366 134.62615 37.04619 200.99640 241.06230 423.52238
## [67] 113.75912 361.37222 80.11718 1574.10821 1200.38203 303.97084
## [73] 491.46259 1160.82451 118.43937 148.43633 392.58891 116.91509
## [79] 1140.18649 543.64472 354.29920 1923.54905 267.50703 215.62236
## [85] 252.58643 21.55404 448.13567 169.30392 764.53825 356.87230
## [91] 151.17818 177.61175 310.80119 443.69131 215.79103 715.29312
## [97] 369.05147 432.29714 204.96476 245.72890 444.90234 871.82295
## [103] 262.35652 1216.27500 1302.79430 188.18380 253.15156 1748.80801
## [109] 139.97416 145.78410 140.09579 3566.77101 242.28989 2642.13520
## [115] 534.07387 121.02665 531.75594 144.63597 1495.85688 487.29555
## [121] 734.33580 896.46719 1825.95416 1021.94188Again, we need a conversion from simple features to spatial polygons:
ant_mun <- as(mun_antioquia, 'Spatial')Uncomment the following code to see what is inside the spatial object:
#head(ant_mun)Now, prepare the plot:
bins <- c(0, 50, 100, 200, 300, 500, 1000, 2000, Inf)
pal <- colorBin("YlOrRd", domain = ant_mun$km2, bins = bins)
labels <- mun_antioquia$NAME_2
labels## [1] "Abejorral" "Abriaquí"
## [3] "Alejandría" "Amagá"
## [5] "Amalfi" "Andes"
## [7] "Angelópolis" "Angostura"
## [9] "Anorí" "Anzá"
## [11] "Apartadó" "Arboletes"
## [13] "Argelia" "Armenia"
## [15] "Barbosa" "Bello"
## [17] "Belmira" "Betania"
## [19] "Betulia" "Bolívar"
## [21] "Briceño" "Buriticá"
## [23] "Cáceres" "Cañasgordas"
## [25] "Caicedo" "Caldas"
## [27] "Campamento" "Caracolí"
## [29] "Caramanta" "Carepa"
## [31] "Carolina del Principe" "Caucasia"
## [33] "Chigorodó" "Cisneros"
## [35] "Cocorná" "Concepción"
## [37] "Concordia" "Copacabana"
## [39] "Dabeiba" "Don Matías"
## [41] "Ebéjico" "El Bagre"
## [43] "El Carmen de Viboral" "El Santuario"
## [45] "Entrerríos" "Envigado"
## [47] "Fredonia" "Frontino"
## [49] "Gómez Plata" "Giraldo"
## [51] "Girardota" "Granada"
## [53] "Guadalupe" "Guarne"
## [55] "Guatapé" "Heliconia"
## [57] "Hispania" "Itagüí"
## [59] "Ituango" "Jardín"
## [61] "Jericó" "La Ceja"
## [63] "La Estrella" "La Unión de Sucre"
## [65] "Liborina" "Maceo"
## [67] "Marinilla" "Medellín"
## [69] "Montebello" "Murindó"
## [71] "Mutatá" "Nariño"
## [73] "Nechí" "Necoclí"
## [75] "Olaya" "Peñol"
## [77] "Pequé" "Pueblorrico"
## [79] "Puerto Berrío" "Puerto Nare"
## [81] "Puerto Triunfo" "Remedios"
## [83] "Retiro" "Rionegro"
## [85] "Sabanalarga" "Sabaneta"
## [87] "Salgar" "San Andrés de Cuerquia"
## [89] "San Carlos" "San Francisco"
## [91] "San Jerónimo" "San José de la Montaña"
## [93] "San Juan de Urabá" "San Luis"
## [95] "San Pedro de los Milagros" "San Pedro de Urabá"
## [97] "San Rafael" "San Roque"
## [99] "San Vicente" "Santa Bárbara"
## [101] "Santa Fe de Antioquia" "Santa Rosa de Osos"
## [103] "Santo Domingo" "Segovia"
## [105] "Sonsón" "Sopetrán"
## [107] "Támesis" "Tarazá"
## [109] "Tarso" "Titiribí"
## [111] "Toledo" "Turbo"
## [113] "Uramita" "Urrao"
## [115] "Valdivia" "Valparaíso"
## [117] "Vegachí" "Venecia"
## [119] "Vigía del Fuerte" "Yalí"
## [121] "Yarumal" "Yolombó"
## [123] "Yondó" "Zaragoza"Then, create the plot:
m <- leaflet(ant_mun) %>%
setView(-75.5, 7, 8) %>% addPolygons(
fillColor = ~pal(km2),
weight = 2,
opacity = 1,
color = "white",
dashArray = "3",
fillOpacity = 0.7,
highlight = highlightOptions(
weight = 5,
color = "#666",
dashArray = "",
fillOpacity = 0.7,
bringToFront = TRUE),
label = labels) %>%
addLegend(pal = pal, values = ~km2, opacity = 0.7, title = NULL,
position = "bottomright")View the plot:
mAnother way of plotting. It may be simpler:
leaflet() %>%
addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
addPolygons(data = ant_mun, popup= ant_mun$NAME_2,
stroke = TRUE, fillOpacity = 0.25, smoothFactor = 0.25
)You can try different providers to improve your map. Take advantage of the tab-completion feature to select the preferred base map just scrolling through the list of 110 providers!
Uncomment the following chunk, complete the code and create a beatiful visualization for the capital of your department. In case of any trouble, use the help of R.
#leaflet() %>%
# addProviderTiles(providers$The following output is just an example:
When your notebook looks fine, publish it on RPubs. Make sure to do it not later than on 25 March.