2. Escala, agregación y distancia

2.1. Escala y resolución

Dentro del análisis de datos espaciales se utiliza frecuentemente el término escala. Este se refiere a la relación entre la distancia en un mapa (normalmente de papel) y la distancia real. Por ejemplo, si 1 cm de distancia representa 100m (100m x 100cm/m=10.000cm) en el mundo real, la relación es 1/10.000, 1:10.000 o 10-4. Esta escala nos permitirá comparar mapas para saber si uno de ellos es relativamente grande en comparación a otro. Dado que hoy en día, el mapeo se hace mayoritariamente en una computadora, la escala es aún más traicionera, debido a que podemos utilizar los mismos datos para hacer mapas de distintos tamaños con una escala diferente. Es así que, en la visualización de mapas, nos solemos preocupar más de la resolución que de la escala.

Recordemos cuando revisamos los datos en ráster, en ellos la resolución está bien definida y viene dada por el tamaño de las celdas. Para el caso de datos vectoriales la resolución puede variar en gran medida dentro de un conjunto de datos, pero se suele pensar en ella como la distancia promedio entre los nodos de las líneas o polígonos. En el caso de los datos de puntos no existe el concepto de resolución, a no ser que varios puntos se combinen como uno solo en varias ocasiones. Dada esta definición es importante tener cuidado, dado que un raster se puede dividir en más celdas de menor tamaño o a un polígono se le puede agregar más nodos para crear una resolución falsa.

En esta imagen, supongamos que el ráster original es el mostrado en el tercer cuadro (4x4 celdas), y es así como se han recogido los datos. Si tomamos las celdas de este y las dividimos a cada una en 4 celdas más pequeñas (segundo cuadro), parecería que hubiésemos aumentado la granularidad de nuestro mapa; sin embargo, los datos originales no se puede dividir equitativamente a dichas celdas. Así, habríamos creado una resolución falsa.

Esta temática es importante ya afecta directamente los análisis que hagamos, sea que estos busquen medir el límite de un territorio (como se muestra en la imagen anterior) o entender variables de interés (las cuales buscaremos recoger en la mayor resolución espacial o granularidad posible).

2.2. Zonificación

Cuando analizamos datos geográficos lo ideal es tenerlos al nivel más granular (y significativo) posible; sin embargo, dado distintas situaciones que deben ser consideradas en el levantamiento de información (principalmente costos), los datos geográficos suelen ser agregados por zonas. Por ejemplo, en lugar de tener datos sobre individuos con o sin empleo, tenemos tasas de desempleo agregadas a nivel provincial.

El nivel de agregación al que los datos son analizados podrían influir sobre los resultados que obtengamos. Es por ello que el problema que los investigadores enfrentan al analizar datos agregados se conoce como inferencia ecológica.

Veamos un ejemplo con datos simulados.

Para este ejemplo, cargaremos primero las librerías necesarias (las cuales corresponden principalmente al tidyverse y a aquellas dedicadas al procesamiento de datos geográficos).

# Librerías necesarias
library(tidyverse)
library(ggpubr)
library(raster)
library(ggrepel)
library(ineq)
# Declaramos la función select para que sea tomada de dplyr y no de raster
select <- dplyr::select

Luego generamos el dato de una variable simulada (ingreso personal) para 1000 personas (puntos) de manera aleatoria, bajo las siguientes premisas:

Las coordenadas x e y se encuentran dentro del rango 0-100.
El ingreso se genera de manera aleatoria, pero de tal forma que los valores más altos caigan cerca de las coordenadas (0,0), (0,100), (100,0) y (100,100); es decir, cerca de las esquinas del territorio.

# Generación de datos aleatorios
set.seed(123)
datos_ingreso <- tibble(x=runif(1000, 0, 100), y=runif(1000, 0, 100)) %>% 
  mutate(ingreso = (runif(1000) * abs((x - 50) * (y - 50))))
# Vistazo a los datos
datos_ingreso %>% glimpse()

## Rows: 1,000
## Columns: 3
## $ x       <dbl> 28.757752, 78.830514, 40.897692, 88.301740, 94.046728, 4.55565~
## $ y       <dbl> 27.362273, 59.386693, 16.018481, 85.343024, 84.773916, 47.7886~
## $ ingreso <dbl> 76.783417, 39.109340, 46.143023, 696.389311, 754.852361, 61.93~

Una vez que tengamos nuestros datos simulados, realizaremos el primer paso necesario en cualquier proceso de análisis de datos: exploración. En este caso particular dicha exploración debe ser geográfica y no geográfica. Para ello veremos los valores del ingreso ordenados de menor a mayor, la distribución del ingreso con un histograma y su distribución geográfica en un gráfico bidimensional.

# Ingreso ordenado
ingreso_ordenado <- datos_ingreso %>% 
  ggplot(aes(x=1:nrow(.),y=sort(ingreso),color=terrain.colors(nrow(.))))+
  geom_point(show.legend = F)+
  labs(x="Índice", y="Ingreso ordenado")
# Distribución del ingreso
ingreso_histograma <- datos_ingreso %>% 
  ggplot(aes(x=ingreso))+
  geom_histogram(bins = 30,
                 fill=rev(terrain.colors(30)),
                 color="black")+
  labs(x="Ingreso", y="Frecuencia")
# Distribución geográfica del ingreso
ingreso_localizado <- datos_ingreso %>% 
  ggplot(aes(x=x, y=y, size=ingreso, color=ingreso))+
  geom_point()+
  scale_color_gradientn(colors=terrain.colors(10))
# Combinación de gráficos
ggarrange(ingreso_ordenado, ingreso_histograma, ingreso_localizado, ncol = 3)

Como podemos notar, son mucho más frecuentes los valores de ingreso bajo. Así mismo, los valores de ingreso alto se ubican en las esquinas del territorio.

Ahora, en base a estos datos, calculemos el índice de Gini. Este índice, creado por Gini (1921) es una medida de económica que sirve para calcular la desigualdad (generalmente de ingresos) que existe entre los ciudadanos de un territorio. El valor del índice de Gini se encuentra entre 0 y 1, siendo cero la máxima igualdad (todos los ciudadanos reciben el mismo ingreso) y uno, la máxima desigualdad (un solo ciudadano recibe todo el ingreso). Formalmente:

\[G = 1-\sum_{k=1}^{n-1}{(X_{k+1}-X_k)(Y_{k+1}+Y_k)}\]

Donde \(G\) es el índice de Gini, \(X\) es la proporción acumulada de población e \(Y\) es la proporción acumulada de ingresos.

Si lo calculamos a detalle, lo podemos hacer con la siguiente tabla:

datos_gini <- datos_ingreso %>% 
  arrange(ingreso) %>% 
  mutate(poblacion=1) %>% 
  mutate(ingreso_acumulado = cumsum(ingreso)/sum(ingreso),
         poblacion_acumulada = cumsum(poblacion)/1000,
         lambda = (lead(poblacion_acumulada,1)-poblacion_acumulada)*(lead(ingreso_acumulado)+ingreso_acumulado))
datos_gini %>% tail()

## # A tibble: 6 x 7
##       x      y ingreso poblacion ingreso_acumulado poblacion_acumulada   lambda
##   <dbl>  <dbl>   <dbl>     <dbl>             <dbl>               <dbl>    <dbl>
## 1 90.3   0.607   1814.         1             0.969               0.995  0.00194
## 2 89.3   2.93    1832.         1             0.975               0.996  0.00196
## 3 10.3   1.92    1894.         1             0.981               0.997  0.00197
## 4 96.7   0.484   1907.         1             0.987               0.998  0.00198
## 5 93.3  99.1     1962.         1             0.993               0.999  0.00199
## 6  5.28 95.4     2016.         1             1                   1     NA

Si calculamos la sumatoria final obtendremos:

datos_gini %>% 
  summarise(G=1-sum(lambda, na.rm = T))

## # A tibble: 1 x 1
##       G
##   <dbl>
## 1 0.588

Otra forma de calcularlo en R viene dada por:

n <- nrow(datos_ingreso)
Gini <- (2 * sum(sort(datos_ingreso$ingreso) * 1:n)/sum(datos_ingreso$ingreso) - (n + 1)) / n
Gini

## [1] 0.5877254

Y si queremos usar una librería:

ineq(datos_ingreso$ingreso, type = "Gini")

## [1] 0.5877254

Este dato nos servirá más adelante para notar los efectos de la agregación de los datos geográficos.

Supongamos entonces que los hogares fueron agrupados por algún tipo de zona censal. Por objetivos de comparación, crearemos varias combinaciones de estos distritos con rásters y calcularemos nuevamente el índice de Gini, basado en datos agregados (ingreso medio).

# Creación de distritos
xy <- datos_ingreso %>% select(x,y) %>% as.matrix()
r1 <- raster(ncol=1, nrow=4, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r1 <- rasterize(xy, r1, datos_ingreso$ingreso, mean)
r2 <- raster(ncol=4, nrow=1, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r2 <- rasterize(xy, r2, datos_ingreso$ingreso, mean)
r3 <- raster(ncol=2, nrow=2, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r3 <- rasterize(xy, r3, datos_ingreso$ingreso, mean)
r4 <- raster(ncol=3, nrow=3, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r4 <- rasterize(xy, r4, datos_ingreso$ingreso, mean)
r5 <- raster(ncol=5, nrow=5, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r5 <- rasterize(xy, r5, datos_ingreso$ingreso, mean)
r6 <- raster(ncol=10, nrow=10, xmn=0, xmx=100, ymn=0, ymx=100, crs=NA)
r6 <- rasterize(xy, r6, datos_ingreso$ingreso, mean)

Con los rásters creados y el ingreso agregado, creemos su visualización geográfica:

# Creamos una lista de gráficos (uno por cada ráster)
raster_plot_list <- map(list(r1,r2,r3,r4,r5,r6),
    function(rasterdata){
      rasterdata %>% 
  as.data.frame(xy=T) %>% 
  ggplot(aes(x,y,fill=layer)) +
  geom_raster()+
  scale_fill_gradientn(colours = rev(terrain.colors(10)))+
  labs(fill="Ingreso medio")
    })
# Unimos todos los gráficos
ggarrange(plotlist = raster_plot_list, common.legend = T,
          ncol = 3, nrow = 2)

Al realizar los gráficos agregados con distintos niveles de ráster se puede observar que mientras más pequeña es la agregación, mejor se captura el patrón real de los datos.

Revisemos ahora si es que capturamos la distribución de los datos reales.

# Creamos una lista de gráficos (uno por cada ráster)
raster_hist_list <- map(list(r1,r2,r3,r4,r5,r6),
    function(rasterdata){
      rasterdata %>% 
  as.data.frame(xy=T) %>% 
  ggplot(aes(layer)) +
  geom_histogram(fill=rev(terrain.colors(30)))+
  labs(x="Ingreso medio", y="Frecuencia")
    })
# Unimos todo en una sola visualización
ggarrange(plotlist = raster_hist_list, common.legend = T,
          ncol = 3, nrow = 2)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En esta gráfica podemos notar que, pese a que la frecuencia con la agregación más pequeña tiene la distribución de ingreso más similar a la original, ésta no es igual.

Finalmente, calculemos el índice de gini para cada agregación.

map(list(r1,r2,r3,r4,r5,r6),
    function(x){
     ineq(x %>% as.data.frame(xy = T) %>% pull("layer"), type = "Gini") 
    })

## [[1]]
## [1] 0.2936912
## 
## [[2]]
## [1] 0.272314
## 
## [[3]]
## [1] 0.04284505
## 
## [[4]]
## [1] 0.3772084
## 
## [[5]]
## [1] 0.4459655
## 
## [[6]]
## [1] 0.4847319

Ninguno de estos valores es igual al 0.5877254 de los datos desagregados, y todos indican un menor nivel de desigualdad al que existe realmente. Una vez que hemos entendido estos efectos, pasemos al estudio de distancias con datos geográficos.

2.3. Distancia

Cuando hablamos de datos geográficos hablamos sin duda alguna de distancias. La distancia es una descripción numérica del espacio existente entre los objetos y es uno de los conceptos fundamentales de la geografía. No por nada la Primera Ley de la Geografía de Waldo Tobler (o principio de autocorrelación espacial) establece que todas las cosas están relacionadas entre sí, pero aquellas más próximas en el espacio tienen una relación mayor que las distantes.

Ahora el problema entre manos es ¿cómo calcular la distancia?. Quizás la primera forma que se nos viene a la mente es la euclídea o distancia en línea recta, pero esta no suele ser relevante. En datos geográficos a veces es necesario considerar distancias en tiempo o distancia de viaje por carretera, por avión, o incluso distancia a pie considerando la curvatura de la tierra y la altitud de las montañas. Cada uno de estos métodos nos dará un resultado distinto, y éste debe estar alineado a nuestros propósitos.

Matriz de distancias

Cuando calculamos distancias, lo más frecuente suele ser colocarlas en una matriz de distancias, la cual contiene un número para cada par de objetos de interés. Si la distancia entre cada par de puntos es simétrica, es necesario llenar tan solo la mitad de la matriz (nótese que este punto se vuelve relevante cuando el número de puntos que debemos analizar se crece).

Para entender cómo se construye una matriz de distancias, realicemos un ejemplo con datos simulados.

# Generación de coordenadas
set.seed(123)
puntos <- tibble(etiqueta = paste("Punto",1:6),
                 x = rnorm(6, 60, 30), y = rnorm(6, 60, 30))
puntos

## # A tibble: 6 x 3
##   etiqueta     x     y
##   <chr>    <dbl> <dbl>
## 1 Punto 1   43.2  73.8
## 2 Punto 2   53.1  22.0
## 3 Punto 3  107.   39.4
## 4 Punto 4   62.1  46.6
## 5 Punto 5   63.9  96.7
## 6 Punto 6  111.   70.8

Dibujemos estos datos en un plano bidimensional.

puntos %>% 
  ggplot(aes(x, y, label=etiqueta))+
  geom_point(color="red", size=4)+
  geom_text_repel()

Con estos datos, podemos construir la matriz de distancias utilizando la función dist.

distancias <- puntos %>% select(x,y) %>% dist(method = "euclidean") %>% as.matrix()
rownames(distancias) <- puntos$etiqueta
colnames(distancias) <- puntos$etiqueta
distancias

##          Punto 1  Punto 2  Punto 3  Punto 4  Punto 5  Punto 6
## Punto 1  0.00000 52.71893 72.30133 33.13642 30.86059 68.33357
## Punto 2 52.71893  0.00000 56.40030 26.18481 75.44895 76.03794
## Punto 3 72.30133 56.40030  0.00000 45.22854 71.59206 31.74843
## Punto 4 33.13642 26.18481 45.22854  0.00000 50.12334 54.93653
## Punto 5 30.86059 75.44895 71.59206 50.12334  0.00000 54.18010
## Punto 6 68.33357 76.03794 31.74843 54.93653 54.18010  0.00000

La función dist por defecto calcula la distancia euclídea. Confirmemos entonces la distancia entre los puntos 3 y 4 utilizando el teorema de pitágoras.

print(puntos %>% filter(etiqueta=="Punto 3"))

## # A tibble: 1 x 3
##   etiqueta     x     y
##   <chr>    <dbl> <dbl>
## 1 Punto 3   107.  39.4

print(puntos %>% filter(etiqueta=="Punto 4"))

## # A tibble: 1 x 3
##   etiqueta     x     y
##   <chr>    <dbl> <dbl>
## 1 Punto 4   62.1  46.6

cateto_x <- abs((puntos %>% filter(etiqueta=="Punto 3") %>% pull(x))-
                  (puntos %>% filter(etiqueta=="Punto 4") %>% pull(x)))
cateto_y <- abs((puntos %>% filter(etiqueta=="Punto 3") %>% pull(y))-
                  (puntos %>% filter(etiqueta=="Punto 4") %>% pull(y)))
hipotenusa <- sqrt(cateto_x^2+cateto_y^2)
hipotenusa

## [1] 45.22854

La distancia calculada con la función dist fue:

distancias[rownames(distancias)=="Punto 3",colnames(distancias)=="Punto 4"]

## [1] 45.22854

Como podemos notar ambos resultados coinciden.

Al final, las matrices de distancias son indispensables para casi cualquier tipo de análisis (incluso no geográfico), dado que la distancia se considera como una medida de similaridad (se recomienda al alumno revisar el concepto de clúster).

Distancias para coordenadas

En datos geográficos podemos calcular distancias utilizando la función pointDistance, que calcula la distancia en línea recta, pero considerando la curvatura de la tierra. Revisemos nuestro último ejemplo.

Cabe notar que esta función siempre utiliza las coordenadas en longitud y latitud, y su resultado se encuentra en metros.

geo_distancia <- pointDistance(puntos %>% select(x,y), lonlat = T)
geo_distancia

##         [,1]    [,2]    [,3]    [,4]    [,5] [,6]
## [1,]       0      NA      NA      NA      NA   NA
## [2,] 5788278       0      NA      NA      NA   NA
## [3,] 5030455 5400260       0      NA      NA   NA
## [4,] 3171815 2845939 3675765       0      NA   NA
## [5,] 2556718 8313655 6390497 5586658       0   NA
## [6,] 2210407 6568843 3505774 3736953 2895177    0

¿Has entendido todos los conceptos? No dudes en hacer preguntas. ¡Nos vemos en el próximo capítulo!

Capítulo 1: Introducción, escala, agregación y distancia

Hugo Porras

2021-09-13

1. Introducción al análisis de datos espaciales