Datos espaciales con R

Este documento es un Notebook de R Markdown. En principio, cuando se ejecuta un trozo de código los resultados se ven debajo del mismo.

Está basado en un documento desarrollado por Robin Lovelace (R.Lovelace@leeds.ac.uk), James Cheshire, Rachel Oldroyd y otros dentro del contexto de la vista espacial del repositorio CRAN https://cran.r-project.org/doc/contrib/intro-spatial-rl.pdf

Comenzamos cargando las librerías necesarias

x <- c("ggmap", "rgdal", "rgeos", "maptools", "dplyr", "tidyr", "tmap","sf")
x

[1] "ggmap"    "rgdal"    "rgeos"    "maptools" "dplyr"    "tidyr"    "tmap"     "sf"      

lapply(x, library, character.only = TRUE) # load the required packages

Google Maps API Terms of Service: http://developers.google.com/maps/terms.
Please cite ggmap if you use it: see citation('ggmap') for details.
package ‘rgdal’ was built under R version 3.4.3rgdal: version: 1.2-16, (SVN revision 701)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
 Path to GDAL shared files: /Users/agus/Library/R/3.4/library/rgdal/gdal
 GDAL binary built with GEOS: FALSE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: /Users/agus/Library/R/3.4/library/rgdal/proj
 Linking to sp version: 1.2-5 
package ‘rgeos’ was built under R version 3.4.2rgeos version: 0.3-26, (SVN revision 560)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
 Linking to sp version: 1.2-5 
 Polygon checking: TRUE 

Checking rgeos availability: TRUE
package ‘dplyr’ was built under R version 3.4.2
Attaching package: ‘dplyr’

The following objects are masked from ‘package:rgeos’:

    intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

package ‘tidyr’ was built under R version 3.4.3package ‘tmap’ was built under R version 3.4.3package ‘sf’ was built under R version 3.4.3Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3

[[1]]
 [1] "ggmap"     "ggrepel"   "ggplot2"   "knitr"     "RJSONIO"   "sp"        "stats"    
 [8] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     

[[2]]
 [1] "rgdal"     "ggmap"     "ggrepel"   "ggplot2"   "knitr"     "RJSONIO"   "sp"       
 [8] "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     

[[3]]
 [1] "rgeos"     "rgdal"     "ggmap"     "ggrepel"   "ggplot2"   "knitr"     "RJSONIO"  
 [8] "sp"        "stats"     "graphics"  "grDevices" "utils"     "datasets"  "methods"  
[15] "base"     

[[4]]
 [1] "maptools"  "rgeos"     "rgdal"     "ggmap"     "ggrepel"   "ggplot2"   "knitr"    
 [8] "RJSONIO"   "sp"        "stats"     "graphics"  "grDevices" "utils"     "datasets" 
[15] "methods"   "base"     

[[5]]
 [1] "dplyr"     "maptools"  "rgeos"     "rgdal"     "ggmap"     "ggrepel"   "ggplot2"  
 [8] "knitr"     "RJSONIO"   "sp"        "stats"     "graphics"  "grDevices" "utils"    
[15] "datasets"  "methods"   "base"     

[[6]]
 [1] "tidyr"     "dplyr"     "maptools"  "rgeos"     "rgdal"     "ggmap"     "ggrepel"  
 [8] "ggplot2"   "knitr"     "RJSONIO"   "sp"        "stats"     "graphics"  "grDevices"
[15] "utils"     "datasets"  "methods"   "base"     

[[7]]
 [1] "tmap"      "tidyr"     "dplyr"     "maptools"  "rgeos"     "rgdal"     "ggmap"    
 [8] "ggrepel"   "ggplot2"   "knitr"     "RJSONIO"   "sp"        "stats"     "graphics" 
[15] "grDevices" "utils"     "datasets"  "methods"   "base"     

[[8]]
 [1] "sf"        "tmap"      "tidyr"     "dplyr"     "maptools"  "rgeos"     "rgdal"    
 [8] "ggmap"     "ggrepel"   "ggplot2"   "knitr"     "RJSONIO"   "sp"        "stats"    
[15] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     

Una vez que tenemos los paquetes vamos a cargar algunos datos, de momento los cogeremos de:

https://github.com/Robinlovelace/Creating-maps-in-R Una vez descargados y descomprimidos ya estamos listos para cargarlos en nuestro proyecto. Daros cuenta de que además de cargar los datos estamos poniendo una serie de condiciones en los mismos, la primera que los enteros de 64 bits nos los convierta en enteros de 32 bits (la opción contraria sería convertirlos a texto) y la segunda es que los atributos de texto no sean factores.

#setwd("/Users/agus/Documents/curso_r/Creating-maps-in-R-master/data")
#lnd <- readOGR(dsn = ".", layer = "london_sport", integer64="allow.loss",stringsAsFactors = FALSE)
#save(lnd, file="lnd.Rdata")
load("/Users/agus/Documents/curso_r/Creating-maps-in-R-master/vignettes/data/lnd.RData")

En la segunda línea del código anterior, la función readOGR se usa para cargar un archivo sho y asignarlo a un nuevo objeto espacial llamado “lnd”; abreviatura de Londres. readOGR es una función que acepta dos argumentos: dsn que significa “nombre de fuente de datos” y especifica el directorio en el que se almacena el archivo, y la capa que especifica el nombre del archivo (tened en cuenta que no es necesario incluir la extensión de archivo .shp ) Los argumentos están separados por una coma y el orden en que se especifican es importante. No hay que escribir explícitamente dsn = o layer = como R sabe qué orden aparecen, por lo que readOGR (“datos”, “london_sport”) funcionaría igual de bien. Para mayor claridad, es una buena práctica incluir nombres de argumento al aprender nuevas funciones, así que continuaremos haciéndolo. El archivo que le asignamos al objeto Lnd contiene la población de los distritos londinenses en 2001 y el porcentaje de la población que participa en actividades deportivas. Estos datos provienen de la Encuesta de personas activas. Los datos de límites provienen de la Ordenance Survey equivalente a nuestro IGN. Para obtener información sobre cómo cargar diferentes tipos de datos espaciales, consultad la documentación de ayuda para readOGR. Se puede llegar a ella escribiendo: ? ReadOGR Para otro ejemplo, en el que se carga una traza de GPS, vea Cheshire y Lovelace (2014).

Los objetos espaciales como ‘Lnd’ están formados por bloques diferentes, los bloques clave son @data (datos de atributos no geográficos) y @polygons (o @lines para datos de línea). El bloque de datos puede considerarse como una tabla de atributos y el de geometría son los polígonos que conforman los límites físicos. A cada uno de ellos se puede acceder usando el símbolo @. Analicemos ahora el objeto ‘lndp’ con algunos comandos básicos:

Siempre podemos conventir los atributos a datos de otra clase

lnd$Pop_2001 <- as.numeric(as.character(lnd$Pop_2001))

Ahora que ya tenemos más o menos claro lo que hay en el archivo vamos a hacer una primera impresión.

plot(lnd)

El que este familiarizado con Londres a lo mejor lo entiende. Como hemos visto en sesiones anteriores una de las cosas buenas que tiene R es la capacidad de seleccionar datos. Por ejemplo, podríamos seleccionar los distritos en los que la participación es inferior al 15%

Pero claro ya que hemos llegado hasta aquí lo que nos gustaría es hacer un mapa que cumpla determinadas condiciones. Vamos a tratar de dibujar los distritos en los que el índice de participación es entre el 20 y 25 %

head(sel) # test output of previous selection (not shown)

[1]  TRUE FALSE  TRUE FALSE  TRUE FALSE

Pero claro ahora, yo que conozco poco Londres, no tengo idea de donde están esos barrios. Por ello, lo mejor sería ponerlos encima del mapa base y pintarlos con algún color. Eso se puede hacer de diferentes maneras pero una es superponiendo el segundo gráfico al primero.

plot(lnd, col = "darkgrey") # plot the london_sport object
sel <- lnd$Partic_Per > 25
plot(lnd[ sel, ], col = "red", add = TRUE) # add selected zones to map

Ahora vamos a ver cuantos de esos distritos están cerca del centro de Londres, para ello calcularemos un radio de 10 Km y lo añadiremos a nuestro mapa. Pero vayamos por partes, primero hallamos el centro de Londres.

file.edit("data/intro-spatial.Rmd")
lat <- coordinates(gCentroid(lnd))[[1]]
long <- coordinates(gCentroid(lnd))[[2]]
longlat<-cbind(long,lat)

Vamos a tratar de quedarnos solo con los distritos que están cerca del centro, una vez hemos hallado el centro, solo tenemos que hacer un área de influencia de 10 Km a su alrededor.

SpatialPoints(coords=cbind(long,lat))

class       : SpatialPoints 
features    : 1 
extent      : 179643.4, 179643.4, 531334, 531334  (xmin, xmax, ymin, ymax)
coord. ref. : NA 

gb<-gBuffer(SpatialPoints(coords=cbind(lat,long)),width=10000,byid = FALSE)
gb@proj4string<-lnd@proj4string
lnd@proj4string

CRS arguments:
 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy
+units=m +no_defs 

plot(lnd)
plot(lnd[ sel, ], col = "turquoise", add = TRUE)

plot(gb,add=TRUE,border=2,lwd=2)

También podemos recortarlos

require(raster)
#plot(intersect(gb,lnd))
plot(gIntersection(gb, lnd))
plot(lnd[ sel, ], col = "turquoise", add = TRUE)

En este caso el resultado es el que buscábamos pero no siempre tiene porque se así. Para eso tenemos que instalar otros paquetes

require(devtools)

Loading required package: devtools
package ‘devtools’ was built under R version 3.4.3

#devtools::install_github("tidyverse/ggplot2")

Ahora podemos intentar dibujar lo que que queríamos haciendo uso de objetos ‘sf’ que puede que sean el futuro de los objetos espaciales en R. Primero dibujamos todos los porcentajes, luego los barrios ‘tocados’ por la selección

plot(st_join(st_as_sf(lnd[,3]),st_as_sf(gb), left = TRUE)) # Sale el conjunto completo

plot(st_join(st_as_sf(lnd[,3]),st_as_sf(gb), left = FALSE))

plot(st_join(st_as_sf(lnd[ sel, ][3]),st_as_sf(gb), left = TRUE))

n same as number of different finite values\neach different finite value is a separate class

plot(st_join(st_as_sf(lnd[ sel, ][3]),st_as_sf(gb), left = FALSE))

n same as number of different finite values\neach different finite value is a separate class

plot(gb,add=TRUE,border=2,lwd=2)

plot(lnd,add=TRUE)

Cuadrantes También podemos hallar los cuadrantes de forma más o menos sencilla, hay paquetes que lo hacen de otra forma, basándonos simplemente en las coordenadas del centro

east <- sapply(coordinates(lnd)[,1], function(x) x > lat)
north <- sapply(coordinates(lnd)[,2], function(x) x > long)

Añadimos esa información a la tabla de atributos

lnd@data$quadrant[east & north] <- "northeast"

Comprobamos que hemos obtenido el resultado que queríamos

lnd@data$quadrant

 [1] NA          NA          NA          "northeast" NA          NA          NA         
 [8] NA          NA          NA          NA          NA          NA          NA         
[15] NA          NA          NA          NA          NA          NA          "northeast"
[22] "northeast" "northeast" NA          "northeast" "northeast" "northeast" NA         
[29] NA          NA          NA          "northeast" "northeast"

Hacemos lo mimos con el resto de cuadrantes

west <- sapply(coordinates(lnd)[,1], function(x) x < lat)
south <- sapply(coordinates(lnd)[,2], function(x) x < long)
lnd@data$quadrant[east & south] <- "southeast"
lnd@data$quadrant[west & north] <- "northwest"
lnd@data$quadrant[south & west] <- "southwest"

Comprobamos el resultado final.

lnd@data$quadrant

 [1] "southeast" "southwest" "northwest" "northeast" "southwest" "southwest" "southwest"
 [8] "southwest" "southwest" "southeast" "southwest" "southeast" "southeast" "southeast"
[15] "northwest" "southwest" "northwest" "northwest" "northwest" "northwest" "northeast"
[22] "northeast" "northeast" "northwest" "northeast" "northeast" "northeast" "southeast"
[29] "southwest" "northwest" "northwest" "northeast" "northeast"

Ya estamos listo para dibujar los cuadrantes. Para ello, convertimos la información en factorial y usamos esos valores para escoger los colores de cada cuadrante, si usamos los mismos valores para dibujar los límites conseguimos un efecto de continuidad que no esta en los datos. Le ponemos algunos adornos.

plot(lnd,col=as.factor(lnd@data$quadrant),border=as.factor(lnd@data$quadrant))
box()

grid()
title("Cuadrantes geográficos", xlab="Este-Oeste",
   ylab="Sur-Norte")

for (i in c(1,2,4))
    {axis(i)}