Trabajando con datos espaciales

En el siguiente documento se abordara un primer acercamiento al uso de datos geoestadisticos en R.

Primer Ejemplo

Para iniciar, se usaran las librerias “Tidyverse” y “Sf” por lo tanto se procede con cargarlas como se muesttra a continuacion:

library(tidyverse)
library(sf)

Una vez cargadas las librerias, se escribe un codigo con el que se busca obtener un mapa en el que se le haran pequeños ajustes para optimizar la visualizacion de la variable que se desea mostrar, en este caso la variable “BIR74”.

system.file("gpkg/nc.gpkg", package = "sf") %>% #se obtuvo un dato del package "sf"
  read_sf() %>% # se esta leyendo el dato obtenido del paquete "sf"
  st_transform(32119) %>% # se transformaron los datos
  select(BIR74) %>% # se selecciono un atributo: nacimientos en 1974
  plot(graticule= TRUE, axes= TRUE)# se coloco la rejilla y los ejes

Ahora se desea ver el contenido de toda la base de datos que contiene la libreria “Sf” respecto al estado utilizado para el mapa de ejemplo, ya que la visualizacion del mapa anterior corresponte a un filtro de datos con fines de visualizacion de una variable en especifico.

read_sf(system.file("gpkg/nc.gpkg", package = "sf"))

Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs

Pero como no vivimos en Estados Unidos, los mapas de ese pais no tienen relevancia en cuanto a lo que se desea aplicar en este documento, por tal motivo se usara la libreria “GADMTools” para obtener informacion de Colombia, por tanto se procede a cargar la libreria:

library(GADMTools)

Loading required package: sp
Loading required package: tidyverse
Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
-- Attaching packages --------------------------------------- tidyverse 1.2.1 --
v ggplot2 3.2.1     v purrr   0.3.2
v tibble  2.1.3     v dplyr   0.8.3
v tidyr   1.0.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.4.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
Loading required package: classInt
Loading required package: sf
Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
Loading required package: rgdal
rgdal: version: 1.4-6, (SVN revision 841)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/Oscar/Documents/R/win-library/3.6/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/Oscar/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.3-1 

Attaching package: 㤼㸱GADMTools㤼㸲

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

    vignette

Obtencion de datos administrativos de un pais

Ahora se procede con obtener informacion de Colombia como se muestra en el siguiente codigo:

(COL= gadm_sf_loadCountries(c("COL"), level = 0, basefile= "./"))

$basename
[1] "./"

$sf
Simple feature collection with 1 feature and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ISO   NAME_0                       geometry
1 COL Colombia MULTIPOLYGON (((-78.97958 1...

$level
[1] 0

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

Grafico de colombia

gadm_plot(COL)

Grafico de Colombia con departamentos

Se muestra el codigo con el que se muestran los departamentos (se cambia “level” de 0 a 1).

(DEPTOS= gadm_sf_loadCountries(c("COL"), level = 1, basefile= "./"))

$basename
[1] "./"

$sf
Simple feature collection with 32 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   ISO   NAME_0    NAME_1       TYPE_1   ENGTYPE_1
1  COL Colombia  Amazonas    Comisaría Commissiary
12 COL Colombia Antioquia Departamento  Department
23 COL Colombia    Arauca  Intendencia  Intendancy
27 COL Colombia Atlántico Departamento  Department
28 COL Colombia   Bolívar Departamento  Department
29 COL Colombia    Boyacá Departamento  Department
30 COL Colombia    Caldas Departamento  Department
31 COL Colombia   Caquetá  Intendencia  Intendancy
32 COL Colombia  Casanare  Intendencia  Intendancy
2  COL Colombia     Cauca Departamento  Department
                         geometry
1  MULTIPOLYGON (((-69.43138 -...
12 MULTIPOLYGON (((-76.99986 8...
23 MULTIPOLYGON (((-69.92316 6...
27 MULTIPOLYGON (((-74.8816 10...
28 MULTIPOLYGON (((-75.79681 1...
29 MULTIPOLYGON (((-74.2158 5....
30 MULTIPOLYGON (((-75.3569 4....
31 MULTIPOLYGON (((-73.88365 -...
32 MULTIPOLYGON (((-72.6363 4....
2  MULTIPOLYGON (((-77.97681 2...

$level
[1] 1

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

gadm_plot(DEPTOS)

Grafico de Colombia con municipios

Se muestra el codigo con el que se muestran los departamentos (se cambia “level” de 1 a 2).

(MUNICIP= gadm_sf_loadCountries(c("COL"), level = 2, basefile= "./"))

$basename
[1] "./"

$sf
Simple feature collection with 1065 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
      ISO   NAME_0    NAME_1           NAME_2                     TYPE_2
15629 COL Colombia  Amazonas       El Encanto Corregimiento Departamento
15634 COL Colombia  Amazonas      La Chorrera Corregimiento Departamento
15699 COL Colombia  Amazonas       La Pedrera Corregimiento Departamento
16578 COL Colombia  Amazonas          Leticia                  Municipio
15505 COL Colombia  Amazonas    Mirití-Paraná Corregimiento Departamento
16577 COL Colombia  Amazonas    Puerto Nariño                  Municipio
15494 COL Colombia  Amazonas Puerto Santander Corregimiento Departamento
16055 COL Colombia  Amazonas         Tarapacá Corregimiento Departamento
14570 COL Colombia Antioquia        Abejorral                  Municipio
14343 COL Colombia Antioquia         Abriaquí                  Municipio
                       ENGTYPE_2                       geometry
15629 Corregimiento Departamento MULTIPOLYGON (((-73.43552 -...
15634 Corregimiento Departamento MULTIPOLYGON (((-71.9802 -1...
15699 Corregimiento Departamento MULTIPOLYGON (((-72.61586 -...
16578               Municipality MULTIPOLYGON (((-70.0579 -4...
15505 Corregimiento Departamento MULTIPOLYGON (((-70.2188 -0...
16577               Municipality MULTIPOLYGON (((-70.28952 -...
15494 Corregimiento Departamento MULTIPOLYGON (((-71.4002 -0...
16055 Corregimiento Departamento MULTIPOLYGON (((-70.1296 -3...
14570               Municipality MULTIPOLYGON (((-75.5753 5....
14343               Municipality MULTIPOLYGON (((-76.0697 6....

$level
[1] 2

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

gadm_plot(MUNICIP)

Extraccion del departamento de estudios posteriores

(boyaca = gadm_subset(DEPTOS, level = 1, regions = "Boyacá"))

$basename
[1] "./"

$sf
Simple feature collection with 1 feature and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -74.6743 ymin: 4.642001 xmax: -71.98309 ymax: 7.0275
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
   ISO   NAME_0 NAME_1       TYPE_1  ENGTYPE_1                       geometry
29 COL Colombia Boyacá Departamento Department MULTIPOLYGON (((-74.2158 5....

$level
[1] 1

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

gadm_plot(boyaca) %>% gadm_showNorth("tl") %>% gadm_showScale("bl")

Si se cambia el dataframe de DEPTOS a MUNICIP salen los municipios

(boyaca = gadm_subset(MUNICIP, level = 1, regions = "Boyacá"))

$basename
[1] "./"

$sf
Simple feature collection with 123 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -74.6743 ymin: 4.642001 xmax: -71.98309 ymax: 7.0275
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
      ISO   NAME_0 NAME_1     NAME_2    TYPE_2    ENGTYPE_2
14856 COL Colombia Boyacá    Almeida Municipio Municipality
14675 COL Colombia Boyacá  Aquitania Municipio Municipality
14607 COL Colombia Boyacá   Arcabuco Municipio Municipality
14528 COL Colombia Boyacá      Belén      <NA>         <NA>
14784 COL Colombia Boyacá     Berbeo Municipio Municipality
14561 COL Colombia Boyacá  Beteitiva Municipio Municipality
14419 COL Colombia Boyacá    Boavita Municipio Municipality
14730 COL Colombia Boyacá     Boyacá Municipio Municipality
14655 COL Colombia Boyacá    Briceño Municipio Municipality
14711 COL Colombia Boyacá Buenavista Municipio Municipality
                            geometry
14856 MULTIPOLYGON (((-73.3883 4....
14675 MULTIPOLYGON (((-73.0553 5....
14607 MULTIPOLYGON (((-73.4648 5....
14528 MULTIPOLYGON (((-72.8194 5....
14784 MULTIPOLYGON (((-73.0639 5....
14561 MULTIPOLYGON (((-72.8499 5....
14419 MULTIPOLYGON (((-72.6702 6....
14730 MULTIPOLYGON (((-73.4572 5....
14655 MULTIPOLYGON (((-73.9081 5....
14711 MULTIPOLYGON (((-73.9409 5....

$level
[1] 2

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

gadm_plot(boyaca) %>% gadm_showNorth("tl") %>% gadm_showScale("bl")

Para saber como se llaman los departamentos

Si se desea saber el nombre de los departamentos o monucipios de la base de datos nacional con la que se esta trabajando, se usa el siguiente codigo (en este caso se usara para mostrar el nombre de los departamentos de Colombia).

listNames(DEPTOS, level = 1)

 [1] "Amazonas"                 "Antioquia"               
 [3] "Arauca"                   "Atlántico"               
 [5] "Bolívar"                  "Boyacá"                  
 [7] "Caldas"                   "Caquetá"                 
 [9] "Casanare"                 "Cauca"                   
[11] "Cesar"                    "Chocó"                   
[13] "Córdoba"                  "Cundinamarca"            
[15] "Guainía"                  "Guaviare"                
[17] "Huila"                    "La Guajira"              
[19] "Magdalena"                "Meta"                    
[21] "Nariño"                   "Norte de Santander"      
[23] "Putumayo"                 "Quindío"                 
[25] "Risaralda"                "San Andrés y Providencia"
[27] "Santander"                "Sucre"                   
[29] "Tolima"                   "Valle del Cauca"         
[31] "Vaupés"                   "Vichada"                 

Recortar un area o saber los limites

Si se desea saber los limites de un territorio que se esta manejando en la base de datos, se usa el siguietne Codigo (en este acaso se usa de ejemplo el departamento de Boyacá).

sbox=gadm_getBbox(boyaca)
sbox

      xmin       ymin       xmax       ymax 
-74.674301   4.642001 -71.983086   7.027500 

Para conocer la estructura de la caja obtenida con el codigo anterior se escribe el siguiente codigo:

str(sbox)

 'bbox' Named num [1:4] -74.67 4.64 -71.98 7.03
 - attr(*, "names")= chr [1:4] "xmin" "ymin" "xmax" "ymax"
 - attr(*, "crs")=List of 2
  ..$ epsg       : int 4326
  ..$ proj4string: chr "+proj=longlat +datum=WGS84 +no_defs"
  ..- attr(*, "class")= chr "crs"

Corte o enfoque de un mapa a partir de sus coordenadas

Si se desea utilizar las coordenadas de un territorio en terminos de ejes “x” y “y” se usa el siguiente codigo:

area_estudio = gadm_crop(DEPTOS, xmin=-75, ymin=4.5, xmax= -71.5, ymax= 7.3)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
attribute variables are assumed to be spatially constant throughout all geometries

area_estudio

$basename
[1] "./"

$sf
Simple feature collection with 11 features and 5 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -75 ymin: 4.5 xmax: -71.5 ymax: 7.3
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   ISO   NAME_0             NAME_1       TYPE_1  ENGTYPE_1
12 COL Colombia          Antioquia Departamento Department
23 COL Colombia             Arauca  Intendencia Intendancy
28 COL Colombia            Bolívar Departamento Department
29 COL Colombia             Boyacá Departamento Department
30 COL Colombia             Caldas Departamento Department
32 COL Colombia           Casanare  Intendencia Intendancy
6  COL Colombia       Cundinamarca Departamento Department
13 COL Colombia               Meta Departamento Department
15 COL Colombia Norte de Santander Departamento Department
20 COL Colombia          Santander Departamento Department
                         geometry
12 POLYGON ((-74.42666 7.3, -7...
23 POLYGON ((-71.5 6.176146, -...
28 POLYGON ((-74.2617 6.985301...
29 POLYGON ((-74.2158 5.4613, ...
30 POLYGON ((-75 5.710153, -74...
32 POLYGON ((-72.9527 4.5, -72...
6  POLYGON ((-73.1357 4.647601...
13 MULTIPOLYGON (((-73.73243 4...
15 POLYGON ((-72.5813 6.8853, ...
20 POLYGON ((-72.6433 6.419401...

$level
[1] 1

$hasBGND
[1] FALSE

attr(,"class")
[1] "gadm_sf"

En la consola aparecen los datos de los departamentos presentes en el recuadro obtenido a partir de las coordenadas escritas en terminos de “x” y “y”, departamentos que se pueden observar si se grafican esas coordenadas como se muestra a continuacion: ### Graficando

gadm_plot(area_estudio)

Si se desea saber el contenido presente en la variable “sf” del conjunto de datos “DEPTOS” se escribe lo siguiente (notese que se puede obtener cada departamento con sus limites en terminos de “x” y “y”):

DEPTOS$sf

Simple feature collection with 32 features and 5 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.84153 ymin: -4.228429 xmax: -66.87033 ymax: 15.91247
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   ISO   NAME_0    NAME_1       TYPE_1   ENGTYPE_1
1  COL Colombia  Amazonas    Comisaría Commissiary
12 COL Colombia Antioquia Departamento  Department
23 COL Colombia    Arauca  Intendencia  Intendancy
27 COL Colombia Atlántico Departamento  Department
28 COL Colombia   Bolívar Departamento  Department
29 COL Colombia    Boyacá Departamento  Department
30 COL Colombia    Caldas Departamento  Department
31 COL Colombia   Caquetá  Intendencia  Intendancy
32 COL Colombia  Casanare  Intendencia  Intendancy
2  COL Colombia     Cauca Departamento  Department
                         geometry
1  MULTIPOLYGON (((-69.43138 -...
12 MULTIPOLYGON (((-76.99986 8...
23 MULTIPOLYGON (((-69.92316 6...
27 MULTIPOLYGON (((-74.8816 10...
28 MULTIPOLYGON (((-75.79681 1...
29 MULTIPOLYGON (((-74.2158 5....
30 MULTIPOLYGON (((-75.3569 4....
31 MULTIPOLYGON (((-73.88365 -...
32 MULTIPOLYGON (((-72.6363 4....
2  MULTIPOLYGON (((-77.97681 2...

Paquete mapview

A continuacion se muestra la instalacion de un paquete util para obtener mapas con otra forma de visualizacion.

install.packages("mapview")

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/
Installing package into 㤼㸱C:/Users/Oscar/Documents/R/win-library/3.6㤼㸲
(as 㤼㸱lib㤼㸲 is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.6/mapview_2.7.0.zip'
Content type 'application/zip' length 2947829 bytes (2.8 MB)
downloaded 2.8 MB

package ‘mapview’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\Oscar\AppData\Local\Temp\RtmpgV327p\downloaded_packages

Usando mapview

Utilizando el conjunto de datos “DEPTOS” se puede obtener e siguiente mapa:

Mapa de Colombia con sus departamentos

suppressPackageStartupMessages(library(mapview))
DEPTOS$sf %>% mapview(zcol="NAME_1",legend=TRUE, col.regions= sf.colors)

Una vez visualizado el mapa de Colombia con departamentos a partir de la libreria mapview, se procedió con la elaboracion del mapa del departamento de Boyacá y sus municipios como se muestra a continuacion (se cambio el conjunto de datos “DEPTOS” por “boyaca” y “NAME_1” por “NAME_2”):

Mapa de Boyacá con sus municipios

suppressPackageStartupMessages(library(mapview))
boyaca$sf %>% mapview(zcol= "NAME_2",legend=TRUE, col.regions= sf.colors)

Conclusiones finales

Para concluir, se pudo observar como se pueden extraaer datos geoespacioales de diferentes bancos de informacion, ademas se mostraron varias formas de graficar la informcion obtenida de estas bancos de datos.