Datos de vectores geoespaciales en R

Angie Alejandra Juyo Buitrago Marzo 24-2020

Este es un cuaderno de R asignado para el curso de geomática básica de la Universidad Nacional de Colombia su objetivo es ilustra las funcionalidades geoespaciales de R con ayuda de las librerias (sf) y (tidyverse)

Lo primero que haremos ser? instalar las librerias mencionadas

#install.packages("tidyverse")
#install.packages("sf")

De forma posterior , debemos cargarlas

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- 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.4
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()

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

Luego usamos la función “read_sf” esta, nos permite leer un archivo que ha sido previamente descargado de DIVA-GIS de forma que representa a los departamentos colombianos

deptos <-  read_sf("C:/Users/LUISA CARRION/Documents/Geomática/COL_adm1.shp")

Podemos saber cuál es el contenido de las primeras filas del objeto que fue importado

head(deptos)

Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -77.149 ymin: -4.228429 xmax: -69.36835 ymax: 11.10792
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

La función st_crs es muy ?til para saber cuál es el sistema de referencia de coordenadas de los datos vectoriales almacenados en el objeto deptos:

st_crs(deptos)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

3. Usando el ggplot para visualizar los datos geoespaciales:

Podemos usar las funcionalidades del ggplot para trazar los datos del objeto deptos

ggplot() + geom_sf(data = deptos) 

Aunque podemos utilizar cualquier sistema de referencia de coordenadas para trazar los datos , no deberiamos utilizar un sistema de referencia de coordenadas que haya sido definido explícitamente para otro país o región

#El CRS 3968 es usado en Virginia , Estados Unidos 
ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3968))

En caso de que necesitemos usar tal SRC, es necesario convertir el objeto espacial de EPSG4326 en EPSG:3968

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

deptos_utm

Simple feature collection with 32 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -280425.6 ymin: -4811689 xmax: 1697588 ymax: -2284762
epsg (SRID):    3968
proj4string:    +proj=lcc +lat_1=37 +lat_2=39.5 +lat_0=36 +lon_0=-79.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

ggplot() + geom_sf(data = deptos_utm)

ggplot() + geom_sf(data = deptos_utm)

4. Filtrar datos geoespaciales en función de los atributos

Esta vez sólo trabajaremos con un departamento,del cual podemos filtrar los datos.

Vichada <-  deptos %>%   filter(NAME_1 == "Vichada")

Y trazamos el nuevo objeto

ggplot() + geom_sf(data = Vichada) 

Podemos repetir los pasos anteriores ,pero esta vez con los municipios correspondientes al departamento de Vichada

munic <-  read_sf("C:/Users/LUISA CARRION/Documents/Geomática/COL_adm2.shp")

mun_Vichada <- munic %>% filter(NAME_1 == "Vichada")
ggplot() + geom_sf(data = mun_Vichada) 

mun_Vichada

Simple feature collection with 6 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -71.0978 ymin: 2.7081 xmax: -67.42111 ymax: 6.301989
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Vichada_points<- st_centroid(mun_Vichada)

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

Vichada_points <- cbind(mun_Vichada, st_coordinates(st_centroid(mun_Vichada$geometry)))

st_centroid does not give correct centroids for longitude/latitude data

Podemos mejorar el aspecto de nuestro gráfico a traves del uso de diferentes funciones, como se muestra en el siguiente bloque de código

ggplot(Vichada) +
    geom_sf() +
    geom_sf(data = Vichada_points, fill = "Turquoise", color = "black") + 
    geom_text(data = Vichada_points, aes(x=X, y=Y,label = ID_2), size = 2) +
    coord_sf(xlim = c(-71.5, -67), ylim = c(2.5, 6.5), expand = FALSE)+ggtitle("Municipios de Vichada")

library(scales)

Attaching package: 㤼㸱scales㤼㸲

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

    discard

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

    col_factor

También,podemos hacer que cada municipio de nuestro departamento se caracterize por un color, en función del código o ID

ggplot(Vichada) + 
  geom_sf(data=Vichada_points, aes(x=X, y=Y, fill =
                                       ID_2), color = "black", size = 0.25) +
  geom_text(data = Vichada_points, aes(x=X, y=Y,label = ID_2), size = 2) +
  theme(aspect.ratio=1)+
  scale_fill_distiller(name="ID_2", palette = "PuBuGn", breaks = pretty_breaks(n = 5))+
  labs(title="Municipios de Vichada")

Ignoring unknown aesthetics: x, y

Aunque, estas representaciones no son mapas reales, podemos exportarlas como PDF o en formato PNG

ggsave("Vichada_municipios.pdf")

Saving 7 x 7 in image

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

Uso de Leaflet para la visualización de datos

Primero,instalamos el programa requerido para posteriormente cargarlo

library(leaflet)

Luego,para usar la libreria debemos convertir nuestros datos en puntos espaciales

Vichada_points <- as(Vichada_points, 'Spatial')

#head(Vichada_points)

install.packages("lwgeom")

library(lwgeom)

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

Podemos obtener el ?rea por municipios de la siguiente manera

mun_Vichada$area <- st_area(mun_Vichada)

Calculamos el área de cada municipio (en metros cuadrados):

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

Creamos un nuevo campo para el area en Kilometros Cuadrados

mun_Vichada$km2

Units: [m^2]
[1] 21531.911 20596.228  9857.561 25241.785 19206.969  4320.376

Finalmente , comprobamos la salida

mun_Vichada$km2

Units: [m^2]
[1] 21531.911 20596.228  9857.561 25241.785 19206.969  4320.376

De nuevo convertimos a poligonos espaciales

Vichada_mun <- as(mun_Vichada, 'Spatial')

Y lo comprobamos

#head(Vichada_mun)

Preparamos nuestra gráfica

bins <- c(0, 50, 100, 200, 300, 500, 1000, 2000, Inf)
pal <- colorBin("YlOrRd", domain = Vichada_mun$km2, bins = bins)


labels <- mun_Vichada$NAME_2

labels

[1] "Cumaribo"          "La Primavera"      "Puerto Carreño"    "San Jose de Ocune" "Santa Rita"       
[6] "Santa Rosalía"    

m <- leaflet(Vichada_mun) %>%
  setView(-70.8, 4, 5)  %>% 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")

m

Otra forma de representar nuestros datos puede ser

leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = Vichada_mun, popup= Vichada_mun$NAME_2,
    stroke = TRUE, fillOpacity = 0.25, smoothFactor = 0.25
  )

Ahora usaremos SetView para conocer Puerto Carreño, la capital de nuestro departamento

v <- leaflet() %>% setView(lng = -69.0, lat = 6.15, zoom = 9)
v %>% addTiles()

NA