Reading, filtering and visualizing geospatial vector data in R

By Valentina Valle Velasco

March 29th, 2020

This R Notebook is in order to teach how to read, filter and plot vector geospatial data.

1st. Introduction.

You should install the simple features (sf) and the tidyverse libraries. You can use the next code:

install.packages(c(“tidyverse”,“sf”))

Next step: We need to load both libraries.

library(tidyverse)

library(sf)

Maybe you will find this when you load the tidyverse library.

## -- Attaching packages ------------------------------------------------------tidyverse 1.3.0  

## v ggplot2 3.3.0     v purrr   0.3.3  
## v tibble  2.1.3     v dplyr   0.8.5  
## v tidyr   1.0.2     v stringr 1.4.0  
## v readr   1.3.1     v forcats 0.5.0

## ---Conflicts--------------------------------------------tidyverse_conflicts()  
## x dplyr::filter() masks stats::filter()                        
## x dplyr::lag()    masks stats::lag()

And this when you load the sf library.

## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2 

2nd. How to read a vector data.

So we are going to read a shape file that we have downloaded before from DIVA-GIS, it is about Colombian departments. You need to use the next code:

deptos <- read_sf("C:/Users/vale_/Desktop/UNAL/6to semestre/COL_adm/COL_adm1.shp")

The code depends on where you have saved the shape file. If you want to know what is the contents of deptos, you should put this code

head(deptos)

It will show you the top lines of a table.
No, we can review basic functions provided by the sf library. For example, you can use the st_crs function to know what is the coordinate reference system of the vector data stored in deptos object.

st_crs(deptos)

3rd. Using ggplot in order to visualize geospatial data:

We are going to use ggplot for plotting the data. This is the code and it will show you a map of Colombia.

ggplot() + geom_sf(data = deptos)

You can choose any coordinate reference system, but this must be defined for your country. The next code allows change the coordinate reference sytem, for example, to CRS 3978 that is used in Canada.

ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3978))

Now, we are going to look for the properties of the CRS with EPSG code 32618. If we need that CRS, we have to convert the spatial object from EPSG4326 into EPSG:32618. For that we will use the next code

deptos_utm <- st_transform(deptos, crs = st_crs(32618))

It will appear a table like the follow, after de code:

deptos_utm

Simple feature collection with 32 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -245935.3 ymin: -469204.3 xmax: 1407491 ymax: 1763314
CRS:            EPSG:32618

ggplot() + geom_sf(data = deptos_utm)

4th. How to filter geospatial data based on attributes

We are going to filter the data because we are just interested in one departament. Note the following code and change it to match your departament:

magdalena <- deptos %>% filter(NAME_1 == "Magdalena")

Now, we must plot the new object:

ggplot() + geom_sf(data = magdalena)

You can repeat the previous steps to load the Colombian municipalities and filter the Magdalena ones:

munic <- read_sf("C:/Users/vale_/Desktop/UNAL/6to semestre/COL_adm/COL_adm2.shp")

Again, the code depend on where you saved the shape file.

mun_magdalena <- munic %>% filter(NAME_1 == "Magdalena")

ggplot() + geom_sf(data = mun_magdalena)

 mun_magdalena

Simple feature collection with 20 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -74.9569 ymin: 8.956799 xmax: -73.6021 ymax: 11.34958
CRS:            4326

magdalena_points <- st_centroid(mun_magdalena)

st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

magdalena_points <- cbind(mun_magdalena, st_coordinates(st_centroid(mun_magdalena$geometry)))

st_centroid does not give correct centroids for longitude/latitude data

We can improve our plot with the next chunk:

ggplot(magdalena) + geom_sf() + geom_sf(data = magdalena_points, fill = "antiquewhite") + geom_text(data = magdalena_points, aes(x=X, y=Y,label = ID_2), size = 2) + coord_sf(xlim = c(-75, -73.55), ylim = c(8.9, 11.4), 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_factor

ggplot(magdalena) + geom_sf(data = magdalena_points, aes(x=X, y=Y, fill = ID_2), color = "white", size = 0.25) + geom_text(data = magdalena_points, aes(x=X, y=Y,label = ID_2), size = 2) + theme(aspect.ratio = 1.5)+ scale_fill_distiller(name="ID_2", palette = "YlGn", breaks = pretty_breaks(n = 5))+ labs(title = "Another Map of Magdalena")

Ignoring unknown aesthetics: x, y

You can save it either as a pdf or as a png.

ggsave("magdalena_municipios.pdf")

Saving 7 x 7 in image

ggsave("map_magdalena.png", width = 6, height = 6, dpi = "screen")

5th. How to visualiza data using leaflet

We need to install the required library:

install.packages(“leaflet”)

Then load the library

library(leaflet)

Now, we need to convert from simple features (sf) to spatial points. Use this code:

ant_points <- as(magdalena_points, 'Spatial')

Uncomment the following command to see what is inside the object:

#head(ant_points)

Get area of municipalities:

install.packages(“lwgeom”)

Load the library:

library(lwgeom)

Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.6.1, PROJ 4.9.3

Attaching package: 㤼㸱lwgeom㤼㸲

The following object is masked from 㤼㸱package:sf㤼㸲:

    st_make_valid

Let’s calculate the area of each municipality (in square meters):

mun_magdalena$area <- st_area(mun_magdalena) #Take care of units

Now, let’s create a new field to store area in square kilometers:

mun_magdalena$km2 <- mun_magdalena$area/(1000000)

Check the output:

mun_magdalena$km2

Units: [m^2]
 [1] 1929.0519 1663.8062
 [3]  279.0380  596.2701
 [5] 1757.0900  798.8020
 [7]  593.2350 1084.1263
 [9]  555.4957  527.5458
[11] 2170.1403 2584.7184
[13]  238.2250  547.6491
[15]  184.8835  417.2063
[17]  271.4884 2211.6646
[19] 2390.0368  746.1562

Again, we need a coversion from simple features to spatial polygons:

ant_mun <- as(mun_magdalena, 'Spatial')

Uncomment the following code to see what is inside the spatial object:

#head(ant_mun)

Now, we are going to 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_magdalena$NAME_2

labels

 [1] "Aracataca"                  
 [2] "Ariguaní"                   
 [3] "Cerro de San Antonio"       
 [4] "Chivolo"                    
 [5] "Ciénaga"                    
 [6] "El Banco"                   
 [7] "El Piñón"                   
 [8] "Fundación"                  
 [9] "Guamal"                     
[10] "Pedraza"                    
[11] "Pivijay"                    
[12] "Plato"                      
[13] "Pueblo Viejo"               
[14] "Remolino"                   
[15] "Salamina"                   
[16] "San Sebastián de Buenavista"
[17] "San Zenón"                  
[18] "Santa Ana"                  
[19] "Santa Marta (Dist. Esp.)"   
[20] "Tenerife"                   

Then, we are goning to create the plot:

m <- leaflet(ant_mun) %>% setView(lng = -74.47, lat = 10.5, zoom = 6.5) %>% addPolygons(fillColor = ~pal(km2), weight = 2, opacity = 1, color = "white", dashArray = "3", fillOpacity = 0.7, highlight = highlightOptions(weight = 5, color = "#666", dashArray = "", bringToFront = TRUE), label = labels) %>% addLegend(pal = pal, values = ~km2, opacity = 0.7, title = NULL, position = "bottomleft")

View the plot:

m

There is an simpler way of plotting, yous should use this code:

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 ist of 110 providers!
Uncomment the following chunk, complete the code and CREATE a beautiful visualization for the capital of your department. In case of any trouble, use the help of R.

santa_marta <- mun_magdalena %>% filter(NAME_2 == "Santa Marta (Dist. Esp.)")

leaflet() %>% addProviderTiles(providers$Esri.WorldImagery, options = providerTileOptions(opacity = 3)) %>% addPolygons(data = santa_marta, popup = santa_marta$NAME_2, stroke = TRUE, fillOpacity = 0.1, smoothFactor = 0.15) %>% setView(lng = -74.201, lat = 11.247, zoom = 13.5)

Make sure all your codes are correctly written, with the spaces, commas and other signs; you have downloaded all the packages and you have the data and the files you need.