Lectura, filtrado y visualización de datos vectoriales geoespaciales en R

Paula Juliana Virguez Gomez 23 de marzo de 2020

1.Introducción

En este notebook se mostraran los pasos a seguir para llevar a cabo la visualización de datos geoespaciales del departamento del Amazonas.

El primer paso sera instalar los paquetes necesarios

install.packages(c("tidyverse", "sf"))

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
probando la URL 'https://cran.rstudio.com/bin/windows/contrib/3.6/tidyverse_1.3.0.zip'
Content type 'application/zip' length 440114 bytes (429 KB)
downloaded 429 KB

probando la URL 'https://cran.rstudio.com/bin/windows/contrib/3.6/sf_0.8-1.zip'
Content type 'application/zip' length 39623129 bytes (37.8 MB)
downloaded 37.8 MB

package ‘tidyverse’ successfully unpacked and MD5 sums checked
package ‘sf’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\bvirguez\AppData\Local\Temp\RtmpS0SakH\downloaded_packages

Ahora hay que cargar ambas librerias:

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.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()

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

2. Lectura de datos vectoriales

En este punto se leera un archivo que representa a los departamentos de Colombia.Estos fueron descargados previamente.

deptos <- read_sf("G:/Geomatica/QGIS/Amazonas/COL_adm1.shp")

Para saber el contenido de deptos escribimos el siguiente comando:

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

Por otro lado en la libreria sf se pueden saber cuales son los sistemas de referencia de coordenadas de los datos vectoriales guardados en deptos usando st_crs

st_crs(deptos)

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

3. Using ggplot to visualize geospatial data:

Ahora usaremos ggplot para poder visualizar los datos

ggplot() + geom_sf(data = deptos)

También se pueden usar otros sistemas de referencia de coordenadas para visualizar los datos

#El CRS 6372 es usado en Mexico
ggplot() + geom_sf(data = deptos) + coord_sf(crs = st_crs(6372))

El codigo EPSG 32618 se usa para buscar las propiedades de CRS

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

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
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

ggplot() + geom_sf(data = deptos_utm)

4.Filtering geospatial data based on attributes

Ahora filtramos los datos solo para el departamento de interes, que en este caso es el amazonas

amazonas <- deptos %>% filter(NAME_1 == "Amazonas")

ggplot() + geom_sf(data = amazonas)

Repetimos los pasos para filtrar solo los municipios del Amazonas

munic <- read_sf("G:/Geomatica/QGIS/Amazonas/COL_adm2.shp")
mun_amazonas <- munic %>% filter(NAME_1 == "Amazonas")
ggplot() + geom_sf(data = mun_amazonas)

mun_amazonas

amazonas_points<- st_centroid(mun_amazonas)

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

amazonas_points <- cbind(mun_amazonas, st_coordinates(st_centroid(mun_amazonas$geometry)))

st_centroid does not give correct centroids for longitude/latitude data

Otro metodos para ejecutar el comando es:

ggplot(amazonas) + 
  geom_sf() +
  geom_sf(data = amazonas_points, fill = "pink") +
  geom_text(data = amazonas_points, aes(x=X, y=Y, label = ID_2), size=2) +
  coord_sf(xlim = c(-74.7, -69), ylim = c(-4.4, 0.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(amazonas) +
  geom_sf(data = amazonas_points, aes(x=X, y=Y, fill =
                                        ID_2), color= "black", size = 0.25) + 
  geom_text(data = amazonas_points,aes(x=X, y=Y, label = ID_2), color="red", size = 3) + 
  theme(aspect.ratio = 1) +
  scale_fill_distiller(name="ID_2", palette = "Blues", breaks = pretty_breaks(n = 5)) +
  labs(title = "Otro mapa del Amazonas")

Ignoring unknown aesthetics: x, y

Para guardar como pdf o png:

ggsave("amazonas_municipios.pdf")

Saving 7 x 7 in image

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

5. Uso de leaflet para visualizar los datos

Primero, intalamos el paquete

install.packages("leaflet")

library(leaflet)

Necesitamos convertir de caracteristicas simples a puntos espaciales

amz_points <- as(amazonas_points, 'Spatial')

head(amz_points)

install.packages("lwgeom")

library(lwgeom)

Calcular el area de cada municipio (m^2):

mun_amazonas$area <- st_area(mun_amazonas)

Ahora el area en km^2

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

Verificar

mun_amazonas$km2

Units: [m^2]
[1] 10276.428 22265.516 24285.214  4172.503 16557.607  1991.018  9042.967 18871.544

Conversion de caracteristicas a poligonos espaciales:

amz_mun <- as(mun_amazonas, 'Spatial')

head(amz_mun)

Ahora, preparar el plot

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


labels <- mun_amazonas$NAME_2

labels

[1] "El Encanto"       "La Chorrera"      "La Pedrera"       "Leticia"          "Mirití-Paraná"    "Puerto Nariño"    "Puerto Santander"
[8] "Tarapacá"        

Crear el plot

m <- leaflet(amz_mun) %>%
  setView(-71.5, -2.2,7) %>% addPolygons(
    fillColor = ~pal(km2),
    weight = 2,
    opacity = 1,
    color = "pink",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 5,
      color = "#111",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels) %>%
  addLegend(pal = pal, values = ~km2, opacity = 0.7, title = NULL,
            position = "bottomright")

Ver del plot

m

Otra opcion para plotear:

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

NA

leaflet(amz_mun) %>%
  setView(-69.94, -4.212,15.5)%>%
  addProviderTiles(providers$Esri.WorldImagery, options = providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = amz_mun, popup = amz_mun$NAME_2,
              stroke = TRUE, fillOpacity = 0.13, smoothFactor = 0.25)