Trabajando con datos georreferenciados en R

Trabajando con datos georreferenciados en R

La idea de este post es mostrar de qué forma R (e imagino que mcuhos otros lenguajes, como Python o Julia, por ejemplo) pueden ayudar a trabajar con datos georreferenciados. No vamos a hacer un tutorial del tema porque hay muchos y muy buenos (por ejemplo, acá, acá o acá). Pero sí queríamos ilustrar la utilidad de R para estos menesteres a partir de un ejemplo bien práctico.

La semana pasada tuve que realizar las siguientes tareas (a partir de una muestra de radios censales proporcionada por un muestrista):

1) seleccionar todos los radios adyacentes

2) calcular una probabilidad de selección de los radios adyacentes en función de la densidad poblacional

3) ordenar aleatoriamente los radios en función de dicha probabilidad

4) plotear cada conjunto de radios colindantes con el orden como etiqueta

Si esto hubiera que hacerlo "manualmente" con algún sistema de GIS podría llevar tiempo -Hay excepciones: QGIS permite trabajar con scripts escritos en Python-. Así que traté de armar una función en R para generar este output. La verdad que resulta muy sencillo editar los datos de un .SHP en R.

Mientras esperaba que el muestrista enviara los radios finales, empecé a preparar el script. La idea era trabajar con una muestra artificial de radios para ir probando cómo hacer.

# IMPORTA LAS LIBRERIAS A USAR #
library(MASS)
library(maptools)
library(rgeos)
library(spdep)
library(igraph)
library(ggmap)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Linking to sp version: 1.2-3
options(stringsAsFactors = F)
setwd("C:/Users/asdf/Documents/asdf/Cartografias/Radios INDEC/Z_TOTAL PAIS")

# IMPORTA EL SHAPE #
radios<-readOGR("Radios_INDEC_total_pais.shp", layer="Radios_INDEC_total_pais")
## OGR data source with driver: ESRI Shapefile 
## Source: "Radios_INDEC_total_pais.shp", layer: "Radios_INDEC_total_pais"
## with 52408 features
## It has 8 fields
# GENERA CAMPO DE DENSIDAD DE POBLACION POR RADIO #
radios$totalpobl[is.na(radios$totalpobl)==TRUE | radios$totalpobl==0]<-0.00001
radios$dens<-radios$totalpobl/gArea(radios,byid = TRUE)*1000

# GENERA MUESTRA ALEATORIA DE RADIOS #
set.seed(45)
index<-sample(1:nrow(radios),8,replace=FALSE)
s_radios<-radios[index,]

La cosa fue relativamente simple: tengo un subset del archivo .SHP de radios censales. A partir de la función gTouches() del paquete rgeos se arma una matriz de adyacencias de radios (todos los que comparten al menos un punto). Luego, se transforma eso a una lista y se "trae" del archivo .SHP original los radios (identificados por su ID) que están en la matriz.

***Esta última parte fue un toque densa y creo que puede mejorarse un poco (al menos estéticamente)***

# FUNCION PARA EL SUBSETEO DE RADIOS #
shp_ady<-function(shp_s,shp){
        library(MASS)
        library(maptools)
        library(rgeos)
        library(spdep)
        library(igraph)
        library(ggmap)
        library(rgdal)

        # GENERA MATRIZ DE ADYACENCIA DE RADIOS #
        ad_radios<-gTouches(shp_s,shp, byid=TRUE)
        row.names(ad_radios)<-shp$link
        colnames(ad_radios)<-shp_s$link

        # GENERA LISTA DE ADYACENCIA DE RADIOS #
        ad_r<-apply(ad_radios,2,function(x){names(x[x==TRUE])})
        ad_r1<-list()
        for (i in 1:length(ad_r)){
                ad_r1[[i]]<-c(names(ad_r)[[i]],unlist(ad_r[[i]]))
        }

        ### LOOP QUE SUBSETEA DEL SHAPE TOTAL LOS RADIOS SELECCIONADOS Y SUS COLINDANTES ###
        shp_col<-list()
        for (i in 1:length(ad_r1)){
                index<-ad_r1[[i]]
                shp_col[[i]]<-radios[radios$link %in% index,]
                index2<-(shp_col[[i]]$link==index[1])
                shp_col[[i]]$p<-0
                shp_col[[i]]$p[index2]<-1
                shp_col[[i]]$p[!index2]<-(shp_col[[i]]$dens[!index2]/sum((shp_col[[i]]$dens[!index2])))
                shp_col[[i]]$orden<-0
                x<-2:length(shp_col[[i]])
                size<-length(shp_col[[i]])-1
                p<-shp_col[[i]]$p[!index2]
                        if (size==1){
                                shp_col[[i]]$orden[!index2]<-2
                                shp_col[[i]]$orden[index2]<-1
                        }
                        else {
                                shp_col[[i]]$orden[!index2]<-sample(x,size,replace=FALSE,prob=p)
                                shp_col[[i]]$orden[index2]<-1
                        }

        }
        return(shp_col)
}
# APLICA LA FUNCION Y TOMA EL TIEMPO #
ptm<-proc.time()
q<-shp_ady(s_radios,radios)
proc.time() - ptm
##    user  system elapsed 
##    0.97    0.00    0.97

Y ahora solamente queda armar una función para poder graficar todo y meterla en un lapply:

# FUNCION QUE GRAFICA LOS RADIOS Y ETIQUETA #
pplot<-function(x){
        plot(x)
        text(coordinates(gCentroid(x,byid=TRUE)),labels=x$orden)
}
# Gráfico de radios seleccionados con sus colindantes #
par(mfrow=c(2,2))
lapply(q,pplot)

Gráfico de radios seleccionados con sus colindantes

Así que espero que con esto se pueda notar lo simple que es trabajar con datos georreferenciados en R.

Comentarios

Entradas populares de este blog

Sobre la multicausalidad

Número efectivo de partidos en elecciones presidenciales 1983-2011