EMPEZANDO CON DATOS ESPACIALES

Este es un cuaderno elaborado con RMarkdown. Cuando ejecutas algun codigo en este cuaderno, los resultados aparecen debajo del codigo.

Intenta ejecutar este bloque de codigo dando click en el bot?n Run en el bloque de codigo o poniendo el cursor sobre este y presionando Ctrl+Shift+Enter.

library(GADMTools)
library(tidyverse)
library(sf)
system.file("gpkg/nc.gpkg", package = "sf") %>%
  read_sf()%>%
  st_transform(32119) %>%
  select(BIR74)%>%
  plot(graticule=T, axes = T)

Con este codigo descargarmos las librerias GADMTools, tidyverse y sf las cuales se usan para obtener la informacion administrativa de los paises, uso y organizacion de listas y uso de “spatial data frames”, repectivamente. Al emplear estas librerias y seleccionar las funciones necesarias logramos visualizar el mapa de “North Caroline”o Carolina del norte, uno de los estados de estados unidos, pero este mapa tambien representamos la informacion del numero de nacimientos por condado que se dieron en 1974 en este estado.

Ahora, para leer un “simple feature” que muestre las caracteristicas asociadas al lugar ademas de la geomatria utilizamos la siguiente linea de codigo:

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

Como se observa, visualizamos la tabla de atributos del spatial data frame cargado con anterioridad, en donde nos muestra inclusive la geometria del mapa elaborado (multipoligono en este caso).

Librerias Geoespaciales

Existen diversas librerias para acceder a diversos conjuntos de datos geoespaciales, como en este caso la libreria “GADMTools” la cual se utliza para obterner los datos administrativos de un pais.

library(GADMTools)

Cargando los datos de un pa?s a nivel administrativos

Ahora que ya tenemos activada la libreria GADMTools, podemos obtener la data geoespacial de Colombia a nivel general (Todo el pais) usando la funcion “gadm_sf_loadCountries” como se observa a continuacion:

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

$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"

Al correr este codigo obtenemos la informacion referente a los datos geoespaciales de Colombia a nivel general (level=0), en donde muestra los extremos norte (ymax), sur(ymin), oriente(xmin) y occidente (xmax), el tipo de geometria entre otros atributos.

Graficando a Colombia en distintos niveles

Ahora que ya tenemos la nformacion del pais que queremos visualizar podemos graficarla con la funcion “gadm_plot” como se muestra a continuacion:

gadm_plot(COL)

De igual manera, podemos visulizar Colombia a niveles mas detallados como por Departamentos y Municipios como se muestra a continuacion:

###Departamentos

(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                       geometry
1  COL Colombia  Amazonas    Comisaría Commissiary MULTIPOLYGON (((-69.43138 -...
12 COL Colombia Antioquia Departamento  Department MULTIPOLYGON (((-76.99986 8...
23 COL Colombia    Arauca  Intendencia  Intendancy MULTIPOLYGON (((-69.92316 6...
27 COL Colombia Atlántico Departamento  Department MULTIPOLYGON (((-74.8816 10...
28 COL Colombia   Bolívar Departamento  Department MULTIPOLYGON (((-75.79681 1...
29 COL Colombia    Boyacá Departamento  Department MULTIPOLYGON (((-74.2158 5....
30 COL Colombia    Caldas Departamento  Department MULTIPOLYGON (((-75.3569 4....
31 COL Colombia   Caquetá  Intendencia  Intendancy MULTIPOLYGON (((-73.88365 -...
32 COL Colombia  Casanare  Intendencia  Intendancy MULTIPOLYGON (((-72.6363 4....
2  COL Colombia     Cauca Departamento  Department MULTIPOLYGON (((-77.97681 2...

$level
[1] 1

$hasBGND
[1] FALSE

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

DEPTOS

$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                       geometry
1  COL Colombia  Amazonas    Comisaría Commissiary MULTIPOLYGON (((-69.43138 -...
12 COL Colombia Antioquia Departamento  Department MULTIPOLYGON (((-76.99986 8...
23 COL Colombia    Arauca  Intendencia  Intendancy MULTIPOLYGON (((-69.92316 6...
27 COL Colombia Atlántico Departamento  Department MULTIPOLYGON (((-74.8816 10...
28 COL Colombia   Bolívar Departamento  Department MULTIPOLYGON (((-75.79681 1...
29 COL Colombia    Boyacá Departamento  Department MULTIPOLYGON (((-74.2158 5....
30 COL Colombia    Caldas Departamento  Department MULTIPOLYGON (((-75.3569 4....
31 COL Colombia   Caquetá  Intendencia  Intendancy MULTIPOLYGON (((-73.88365 -...
32 COL Colombia  Casanare  Intendencia  Intendancy MULTIPOLYGON (((-72.6363 4....
2  COL Colombia     Cauca Departamento  Department MULTIPOLYGON (((-77.97681 2...

$level
[1] 1

$hasBGND
[1] FALSE

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

gadm_plot(DEPTOS)

Municipios

(MUNC<-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                  ENGTYPE_2
15629 COL Colombia  Amazonas       El Encanto Corregimiento Departamento Corregimiento Departamento
15634 COL Colombia  Amazonas      La Chorrera Corregimiento Departamento Corregimiento Departamento
15699 COL Colombia  Amazonas       La Pedrera Corregimiento Departamento Corregimiento Departamento
16578 COL Colombia  Amazonas          Leticia                  Municipio               Municipality
15505 COL Colombia  Amazonas    Mirití-Paraná Corregimiento Departamento Corregimiento Departamento
16577 COL Colombia  Amazonas    Puerto Nariño                  Municipio               Municipality
15494 COL Colombia  Amazonas Puerto Santander Corregimiento Departamento Corregimiento Departamento
16055 COL Colombia  Amazonas         Tarapacá Corregimiento Departamento Corregimiento Departamento
14570 COL Colombia Antioquia        Abejorral                  Municipio               Municipality
14343 COL Colombia Antioquia         Abriaquí                  Municipio               Municipality
                            geometry
15629 MULTIPOLYGON (((-73.43552 -...
15634 MULTIPOLYGON (((-71.9802 -1...
15699 MULTIPOLYGON (((-72.61586 -...
16578 MULTIPOLYGON (((-70.0579 -4...
15505 MULTIPOLYGON (((-70.2188 -0...
16577 MULTIPOLYGON (((-70.28952 -...
15494 MULTIPOLYGON (((-71.4002 -0...
16055 MULTIPOLYGON (((-70.1296 -3...
14570 MULTIPOLYGON (((-75.5753 5....
14343 MULTIPOLYGON (((-76.0697 6....

$level
[1] 2

$hasBGND
[1] FALSE

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

MUNC

$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                  ENGTYPE_2
15629 COL Colombia  Amazonas       El Encanto Corregimiento Departamento Corregimiento Departamento
15634 COL Colombia  Amazonas      La Chorrera Corregimiento Departamento Corregimiento Departamento
15699 COL Colombia  Amazonas       La Pedrera Corregimiento Departamento Corregimiento Departamento
16578 COL Colombia  Amazonas          Leticia                  Municipio               Municipality
15505 COL Colombia  Amazonas    Mirití-Paraná Corregimiento Departamento Corregimiento Departamento
16577 COL Colombia  Amazonas    Puerto Nariño                  Municipio               Municipality
15494 COL Colombia  Amazonas Puerto Santander Corregimiento Departamento Corregimiento Departamento
16055 COL Colombia  Amazonas         Tarapacá Corregimiento Departamento Corregimiento Departamento
14570 COL Colombia Antioquia        Abejorral                  Municipio               Municipality
14343 COL Colombia Antioquia         Abriaquí                  Municipio               Municipality
                            geometry
15629 MULTIPOLYGON (((-73.43552 -...
15634 MULTIPOLYGON (((-71.9802 -1...
15699 MULTIPOLYGON (((-72.61586 -...
16578 MULTIPOLYGON (((-70.0579 -4...
15505 MULTIPOLYGON (((-70.2188 -0...
16577 MULTIPOLYGON (((-70.28952 -...
15494 MULTIPOLYGON (((-71.4002 -0...
16055 MULTIPOLYGON (((-70.1296 -3...
14570 MULTIPOLYGON (((-75.5753 5....
14343 MULTIPOLYGON (((-76.0697 6....

$level
[1] 2

$hasBGND
[1] FALSE

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

gadm_plot(MUNC)

Extraccion de un departamento

Ademas de las anteriores funciones, tambien es posible visualizar mapas un poco mas especificos. A continuacion mostramos como, mediante la funcion “gadm_subset”, se pueden pueden extraer los datos para visualizar solo el departamento de interes, que en este caso es el Huila:

HUILA<-gadm_subset(DEPTOS, level = 1, regions = "Huila")
gadm_plot(HUILA) %>% gadm_showNorth("tl") %>% gadm_showScale("bl")

Este mismo proceso lo podemos realizar para extraer y visualizar los municipios del departamento como se muestra a continuacion:

HUMUN<- gadm_subset(MUNC, level = 1, regions = "Huila");str(HUMUN)

List of 4
 $ basename: chr "./"
 $ sf      :Classes ‘sf’ and 'data.frame':  37 obs. of  7 variables:
  ..$ ISO      : chr [1:37] "COL" "COL" "COL" "COL" ...
  ..$ NAME_0   : chr [1:37] "Colombia" "Colombia" "Colombia" "Colombia" ...
  ..$ NAME_1   : chr [1:37] "Huila" "Huila" "Huila" "Huila" ...
  ..$ NAME_2   : chr [1:37] "Acevedo" "Agrado" "Aipe" "Algeciras" ...
  ..$ TYPE_2   : chr [1:37] "Municipio" "Municipio" "Municipio" "Municipio" ...
  ..$ ENGTYPE_2: chr [1:37] "Municipality" "Municipality" "Municipality" "Municipality" ...
  ..$ geometry :sfc_MULTIPOLYGON of length 37; first list element: List of 1
  .. ..$ :List of 1
  .. .. ..$ : num [1:64, 1:2] -75.8 -75.8 -75.8 -75.9 -75.9 ...
  .. ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
  ..- attr(*, "sf_column")= chr "geometry"
  ..- attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
  .. ..- attr(*, "names")= chr [1:6] "ISO" "NAME_0" "NAME_1" "NAME_2" ...
 $ level   : num 2
 $ hasBGND : logi FALSE
 - attr(*, "class")= chr "gadm_sf"

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

class(HUMUN)

[1] "gadm_sf"

Listados

Cuando usamos este tipo de datos resulta ?til conocer distintos componentes que constituyen nuestro objeto de interes, por lo que a continuacion mostramos con ejemplos como generar estas listas, ya sea usando una funcion (“listNames” en este caso) o llamanndo especificamente al componente deseado:

LDEPTOS<-listNames(DEPTOS,1);LDEPTOS

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

LMUNC<- HUMUN$sf$NAME_2;LMUNC

 [1] "Acevedo"      "Agrado"       "Aipe"         "Algeciras"    "Altamira"     "Baraya"       "Campoalegre" 
 [8] "Colombia"     "Elías"        "Garzón"       "Gigante"      "Guadalupe"    "Hobo"         "Iquira"      
[15] "Isnos"        "La Argentina" "La Plata"     "Nátaga"       "Neiva"        "Oporapa"      "Paicol"      
[22] "Palermo"      "Palestina"    "Pital"        "Pitalito"     "Rivera"       "Saladoblanco" "San Agustín" 
[29] "Santa María"  "Suaza"        "Tarquí"       "Tello"        "Teruel"       "Tesalia"      "Timaná"      
[36] "Villavieja"   "Yaguará"     

Recortando areas

Con unos cuantos codigos adicionales es posible realizar un mapa mas elaborado. El siguiente bloque de codigos se usa para visualizar el mapa del departamento con los departamentos con los que limita, todo visualizado en una cuadricula en la cual se pueden regular sus limites en funcion de los limites que presenta el mapa.

(sbox <- gadm_getBbox(HUILA))

      xmin       ymin       xmax       ymax 
-76.640198   1.493471 -74.505402   3.807901 

Este codigo sirve para exhibir los limites cardinales que presenta el departamento.

str(sbox)

 'bbox' Named num [1:4] -76.64 1.49 -74.51 3.81
 - 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"

Aqui la funcion “str” es usada para obtener la informacion sobre el tipo de objeto que es el nuevo data frame generado anteriormente.

(STUDY_AREA<- gadm_crop(DEPTOS, xmin =-76.66, ymin = 1.49, xmax = -74.49, ymax = 3.81))

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

$basename
[1] "./"

$sf
Simple feature collection with 7 features and 5 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.66 ymin: 1.49 xmax: -74.49 ymax: 3.81
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
   ISO   NAME_0          NAME_1       TYPE_1  ENGTYPE_1                       geometry
31 COL Colombia         Caquetá  Intendencia Intendancy POLYGON ((-76.20522 1.49, -...
2  COL Colombia           Cauca Departamento Department POLYGON ((-76.66 3.108991, ...
6  COL Colombia    Cundinamarca Departamento Department POLYGON ((-74.49 3.655067, ...
9  COL Colombia           Huila Departamento Department POLYGON ((-76.13251 1.49347...
13 COL Colombia            Meta Departamento Department POLYGON ((-74.49 1.817293, ...
22 COL Colombia          Tolima Departamento Department POLYGON ((-74.809 3.493, -7...
24 COL Colombia Valle del Cauca Departamento Department POLYGON ((-76.5175 3.282401...

$level
[1] 1

$hasBGND
[1] FALSE

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

Con la funcion “gadm_crop” se pueden establecer los limites que va a tener la cuadricula en la que se va a plasma el mapa a desarrollar.

gadm_plot(STUDY_AREA)

Finalmente, este es el mapa obtenido, en donde se puede observar al departamento del Huila rodeado por los departamentos vecinos.

Obtencion de un “Simple feature (sf)”

Cuando se estemos trabajando con datos geograficos, resutla util la obtencion del simple feature, ya que este contiene la informacion general del mapa desarrollado.

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                       geometry
1  COL Colombia  Amazonas    Comisaría Commissiary MULTIPOLYGON (((-69.43138 -...
12 COL Colombia Antioquia Departamento  Department MULTIPOLYGON (((-76.99986 8...
23 COL Colombia    Arauca  Intendencia  Intendancy MULTIPOLYGON (((-69.92316 6...
27 COL Colombia Atlántico Departamento  Department MULTIPOLYGON (((-74.8816 10...
28 COL Colombia   Bolívar Departamento  Department MULTIPOLYGON (((-75.79681 1...
29 COL Colombia    Boyacá Departamento  Department MULTIPOLYGON (((-74.2158 5....
30 COL Colombia    Caldas Departamento  Department MULTIPOLYGON (((-75.3569 4....
31 COL Colombia   Caquetá  Intendencia  Intendancy MULTIPOLYGON (((-73.88365 -...
32 COL Colombia  Casanare  Intendencia  Intendancy MULTIPOLYGON (((-72.6363 4....
2  COL Colombia     Cauca Departamento  Department MULTIPOLYGON (((-77.97681 2...

Paquete “mapview”

El paquete “mapview” es una herramienta util que permite una visualizacion mas dinamica de los mapas elaborados, permitiendo asignar colores a las divisiones que tenga el mapa o simplemente alejar o acercar el mapa. Como ejemplo se muestran los siguientes mapas:

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/usuagro/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 2947906 bytes (2.8 MB)
downloaded 2.8 MB

package ‘mapview’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\usuagro\AppData\Local\Temp\Rtmpcvc5IS\downloaded_packages

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

HUMUN$sf %>% mapview(zcol="NAME_2", legend= T, col.regions=sf.colors)