Geoinformática - Práctica 5¶

Autocorrelación Espacial y Análisis Exploratorio de Datos Espaciales¶

La Autocorrelación Espacial se relaciona con el grado en el que la similitud de valores entre observaciones de un conjunto de datos se conecta con la similitud en la localización de dichas observaciones. De forma muy semejante a como ocurre con la Correlación Tradicional (que determina cómo los valores de una variable cambian en función de los valores de otra) y de forma análoga a su homólogo en una serie temporal (que relaciona el valor de una variable en algún punto del tiempo con los valores en periodos anteriores), la Autocorrelación Espacial relaciona los valores de la variable de interés en un lugar dado con los valores de ella misma en sus alrededores.

Una idea clave dentro de este contexto es el de Aleatoriedad Espacial, esto es, una situación en la cual la ubicación de una observación no guarda relación alguna con el valor de una variable; en otras palabras, una variable es aleatoria espacialmente si la distribución que sigue en el espacio parece no tener algún patrón discernible. Como tal, la Autocorrelación Espacial puede ser definida formalmente como la "absencia de Aleatoriedad Espacial", y, al igual que con la correlación tradicional, puede presentarse en dos variaciones:

  • Correlación Espacial Positiva - Donde los valores similares tienden a agruparse entre sí en ubicaciones similares.
  • Correlación Espacial Negativa - Donde los valores similares tienden a estar dispersos y alejados entre sí.

En esta práctica, se aprenderá a explorar la Autocorrelación Espacial de un conjunto de datos determinado, analizando la presencia, naturaleza y robustez de los datos utilizados. Para esto, se utilizará un conjunto de herramientas colectivamente conocidos como Análisis Exploratorio de Datos Espaciales (ESDA, por sus siglas en inglés), diseñado específicamente para este propósito. El alcance de los Métodos ESDA es bastante amplio y abarca desde aproximaciones sencilllas, como Coropletas (Choropleths) y Consultas a Bases de Datos, hasta metodologías más avanzadas y robustas, como la Inferencia Estadística y la Identificación de la Dimensión Geográfica de los Datos. El propósito de esta sesión es explorar más a detalle estas últimas técnicas.

Los Métodos ESDA comúnmente se dividen en dos grupos principales: herramientas para analizar a la autocorrelación espacial a nivel global, y otras a nivel local. Las primeras consideran el patrón general que sigue la localización de los valores, y hacen posible realizar conclusiones sobre el grado de Agrupamiento (Clustering) que siguen los datos; pueden responderse preguntas como ¿Los valores en general siguen un patrón en particular en su distribución geográfica? o ¿Los valores similares se encuentran más cerca de lo que se encontrarían en una distribución aleatoria?. Se analizará la Autocorrelación Espacial Global utilizando el estadístico de I de Moran.

Los métodos utilizados para la Autocorrelación Espacial Local se enfocan en la inestabilidad espacial, esto es, que tanto se alejan ciertas ubicaciones en particular del patrón general; la idea detrás de esto es que, aún cuando existe un patrón general en términos de la naturaleza y robustez de la asociación espacial, algunas áreas en particular pueden divergir sustancialmente de éste mismo. Sin importar el grado general de concentración de los valores, existen focos donde valores inusualmente altos (o bajos) se encuentran cercanos a otros valores altos (o bajos), zonas que son conocidas como Hotspots (o Coldspots, para valores bajos); además, también es posible observar valores muy altos rodeados de valores muy bajos (o viceversa), los cuales reciben el nombre de Outliers Espaciales (Valores Atípicos). La ténica principal que se utilizará en esta sesión para explorar la Autocorrelación Espacial Local serán los Indicadores Locales de Asociación Espacial (LISA, por sus siglas en inglés).

Nota: Para instalar las dependencias en Colab, copien las siguientes líneas en una celda y ejecútenla.

# Important library for many geopython libraries
!apt install gdal-bin python-gdal python3-gdal 
# Install rtree - Geopandas requirment
!apt install python3-rtree 
# Install Geopandas
!pip install git+git://github.com/geopandas/geopandas.git
# Install descartes - Geopandas requirment
!pip install descartes 
# Install Folium for Geographic data visualization
!pip install folium
# Install plotlyExpress
!pip install mapclassify
!pip install palettable
!pip install pysal
!pip install esda
In [20]:
%matplotlib inline

import seaborn as sns
import pandas as pd
from libpysal.weights import W,Queen, Rook, KNN, DistanceBand, min_threshold_distance, block_weights, lag_spatial
import esda
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import mapclassify

np.random.seed(123)

Datos¶

En esta práctica, se continuarán utilizando las AGEB's como unidad espacial, así como la variable de Personas sin Derechohabiencia a Servicios Públic de Salud colectada por el Consejo Nacional de Evaluación de Política del Desarrollo Social (CONEVAL) en 2010.

NOTA: Bajar datos para trabajar en Colab

url = "https://www.dropbox.com/s/idewfb8v9stbt7b/pob_sinderechohab.zip?dl=1"
r = requests.get(url, allow_redirects=True)
open('/content/pob_sinderechohab.zip', 'wb').write(r.content)

Como es de costumbre, se creará una variable que contenga la ruta donde se encuentran los datos a utilizar, si van a trabajar en Colab, ajusten el valor:

In [2]:
f = 'data/'

Importando y Analizando los Datos¶

Como se ha hecho siempre, se utiliza GeoPandas para importar e inspeccionar los datos espaciales de la práctica:

In [4]:
# Importar el ShapeFile
agebs = gpd.read_file(f + 'pob_sinderechohab.shp')

# Graficar el GeoDataFrame
agebs.plot(cmap = 'Set1', figsize = (10,10))
Out[4]:
<AxesSubplot:>

Antes de continuar con la parte analítica, resulta pertinente crear un Mapa de Coropletas para, de entrada, tener una idea visual de cómo es que se distribuye la variable en el espacio. A continuación, se muestra el código necesario para construir este mapa, en el cual cabe resaltar el uso del argumento scheme dentro de la función .plot() para indicar cómo deben de ser clasificados los datos, para después asignarles el color correspondiente:

In [5]:
# Preparación de la figura y sus filas
fig , filas = plt.subplots(1, figsize = (10,10))

# Gráfica del GeoDataFrame, generando un Mapa de Coropletas a través de 'scheme'
agebs.plot(ax = filas, column = 'p_sderech', scheme = 'fisherjenks' , cmap = 'YlGn', legend = True)

# Adición de título a la gráfica
fig.suptitle('Población sin Derechohabiencia a Servicios de Salud Públicos en la CDMX', size = 20)

# Mostrar la gráfica
plt.show()

Ejercicio Opcional¶

Crea un mapa similar al anterior, pero utilizando como Métodos de Clasificación los Cuantiles (quantiles) e Intervalor Iguales (equalinterval), y responde las preguntas: ¿En qué difieren los mapas generados? ¿Cómo varía la interpretación de lo observado en cada uno de los casos? ¿Cuál es la distribución de los datos? Confirma esta última respuesta generando un Diagrama de Densidad de Kernal o un Histograma.

Matriz de Pesos Espaciales¶

Como se ha discutido anteriormente, una Matriz de Pesos Espaciales es la forma en la que el espacio geográfico es codificado de forma numérica, a modo de que sea fácil para una computadora (o un método estadístico) el interpretarlo. Se han estudiado algunas de las varias formas conceptuales en las que puede definirse una Matriz de Pesos Espaciales, como los criterios de Contigüidad, Basados en Distancia, o de Bloques.

En este caso, se utilizará uno de los criterios más comúnes en este tipo de análisis: la Matriz de Contigüidad de Reina:

In [7]:
m_reina = Queen.from_dataframe(agebs, idVariable='cvegeo')
m_reina.transform = 'R'

('WARNING: ', '090090015012A', ' is an island (no neighbors)')

/home/plablo/miniconda3/envs/covid/lib/python3.10/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected: 
 There are 8 disconnected components.
 There is 1 island with id: 090090015012A.
  warnings.warn(message)

Cabe recordar que, al modificar el parámetro transform de la matriz, los valores de los pesos asignados se alteran; por ejemplo, al modificar este parámetro por 'R', se han estandarizado los pesos de modo que la suma de éstos sea igual a uno:

In [8]:
m_reina['090150001056A']
Out[8]:
{'0901500010517': 0.25,
 '0901500010610': 0.25,
 '0901500010574': 0.25,
 '0901500010502': 0.25}

También necesitamos quitar la isla que detectamos, más adelante esta sla nos va a impedir realizar algunos cálculos

In [27]:
agebs = agebs.set_index('cvegeo').drop('090090015012A').reset_index()

Volvemos a calcular nuestra matriz de vecindad

In [28]:
m_reina = Queen.from_dataframe(agebs, idVariable='cvegeo')
m_reina.transform = 'R'

/home/plablo/miniconda3/envs/covid/lib/python3.10/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected: 
 There are 7 disconnected components.
  warnings.warn(message)

Rezago Espacial¶

Una vez que se encuentran listos los datos y la Matriz de Pesos Espaciales, es posible calcular el Rezago Espacial de la variable p_sderech. Importante recordar que el rezago espacial es el resultado del producto entre la Matriz de Pesos Espaciales y una variable dada y que, si $W$ se encuentra estandarizada en filas, los resultados equivalen al valor promedio de una variable en la vecidad de una observación dada.

En la siguiente línea, se calcula el Rezago Espacial para la variable p_sderech y es almacenada directamente en la tabla original:

In [9]:
agebs['rez_espacial'] = lag_spatial(m_reina , agebs['p_sderech'])

Puede compararse entre sí la apariencia del rezago espacial calculado con respecto a su variable original:

In [10]:
agebs[['p_sderech' , 'rez_espacial']].head()
Out[10]:
p_sderech rez_espacial
0 389.0 337.0
1 323.0 419.5
2 448.0 484.5
3 240.0 409.0
4 1017.0 500.0

La forma en la que se interpreta rez_espacial para una observación dada es la siguiente: el AGEB con Clave '0901500010235' (la segunda en la lista anterior) tiene una Población sin Derechohabiencia de 323 personas, y se encuentra rodeada por otras AGEB's que, en promedio, tienen 419.5 Pesonas sin Derechohabiencia.

Para ejemplificar lo anterior correctamente, pueden obtenerse primero los vecinos de la AGEB mencionada:

In [11]:
m_reina.neighbors['0901500010235']
Out[11]:
['090150001024A',
 '0901500010076',
 '0901500010377',
 '0901500010362',
 '0901500010381',
 '0901500010220']

Y después, verificar los valores de las variables para estos vecinos:

In [14]:
vecinos = agebs.set_index('cvegeo').loc[m_reina.neighbors['0901500010235'] , 'p_sderech']
vecinos
Out[14]:
cvegeo
090150001024A     226.0
0901500010076    1428.0
0901500010377     226.0
0901500010362     324.0
0901500010381     267.0
0901500010220     146.0
Name: p_sderech, dtype: float64

De forma sencilla, puede calcularse el valor promedio de estos valores, en cuyo caso podemos comprobar que es muy similar al Rezago Espacial calculado para dicha observación:

In [15]:
vecinos.mean()
Out[15]:
436.1666666666667

Para algunas técnicas que se estudiarán más adelante, tiene más sentido el operar con la versión estandarizada de la variable, en lugar de sus valores brutos; estandarizar significa sencillamente el restar el valor de la media de la variable y dividirlo por la desviación estándar de cada observacion de la columna. Esto se consigue rápidamente a través de los siguientes comandos:

In [16]:
agebs['p_sderech_std'] = (agebs['p_sderech'] - agebs['p_sderech'].mean()) / agebs['p_sderech'].std()
agebs.head()
Out[16]:
cvegeo alcaldia p_sderech p_nescu p_hacin p_analf geometry rez_espacial p_sderech_std
0 0900700013628 09007 389.0 475.0 0.0 9.0 POLYGON ((2810132.372 824698.172, 2810169.540 ... 337.0 -0.959260
1 0901500010235 09015 323.0 481.0 6.0 7.0 POLYGON ((2798881.634 831643.241, 2798843.076 ... 419.5 -1.034395
2 0900200010097 09002 448.0 780.0 5.0 37.0 POLYGON ((2792415.239 836846.390, 2792356.808 ... 484.5 -0.892094
3 0900200011184 09002 240.0 389.0 0.0 11.0 POLYGON ((2792260.139 836768.777, 2792333.695 ... 409.0 -1.128883
4 0900300011285 09003 1017.0 1291.0 0.0 23.0 POLYGON ((2802121.600 817466.682, 2802124.157 ... 500.0 -0.244339

Y, de forma análoga, debe de calcularse el Rezago Espacial para estos valores estandarizados:

In [17]:
agebs['rez_espacial_std'] = ps.lib.weights.lag_spatial(m_reina , agebs['p_sderech_std'])
agebs.head()
Out[17]:
cvegeo alcaldia p_sderech p_nescu p_hacin p_analf geometry rez_espacial p_sderech_std rez_espacial_std
0 0900700013628 09007 389.0 475.0 0.0 9.0 POLYGON ((2810132.372 824698.172, 2810169.540 ... 337.0 -0.959260 -1.018458
1 0901500010235 09015 323.0 481.0 6.0 7.0 POLYGON ((2798881.634 831643.241, 2798843.076 ... 419.5 -1.034395 -0.924539
2 0900200010097 09002 448.0 780.0 5.0 37.0 POLYGON ((2792415.239 836846.390, 2792356.808 ... 484.5 -0.892094 -0.850542
3 0900200011184 09002 240.0 389.0 0.0 11.0 POLYGON ((2792260.139 836768.777, 2792333.695 ... 409.0 -1.128883 -0.936492
4 0900300011285 09003 1017.0 1291.0 0.0 23.0 POLYGON ((2802121.600 817466.682, 2802124.157 ... 500.0 -0.244339 -0.832897

Autocorrelación Espacial Global¶

La Autocorrelación Espacial Global se relaciona con el patrón geográfico general presente en un conjunto de datos. Los estadísticos diseñados para medir ésta caracterizan a los datos espaciales en términos de su grado de agrupamiento (clustering) y lo resumen. Este resumen puede ser visual o numérico; en esta sección, se estudiará un ejemplo de cada uno de éstos: el Gráfico de Moran, y el Estadístico de I de Moran de Autocorrelación Espacial

Gráfico de Moran¶

El Gráfico de Moran es una forma de visualizar un conjunto de datos espaciales para explorar la naturaleza y robustez de la Autocorrelación Espacial. En escencia se trata de un Diagrama de Dispersión tradicional en el que la variable de interés se contrasta contra su Rezago Espacial. Para poder interpretar los valores en términos de la media, y sus cantidades en términos de desviaciones estándar, la variable de interés usualmente se estandariza.

Desde el punto de vista técnico, crear un Gráfico de Moran en Python sigue un proceso muy similar al de calcular cualquier otro Diagrama de Dispersión, dado que se tenga la variable estandarizada y su Rezago Espacial de antemano:

In [18]:
# Preparación de la Gráfica
fig, filas = plt.subplots(1, figsize = (10,10))

# Generación del Diagrama de Dispersión
sns.regplot(x = 'p_sderech_std' , y = 'rez_espacial_std', data = agebs)

# Adición de lineas horizontal y vertical
plt.axvline(0, c='k', alpha=0.5)
plt.axhline(0, c='k', alpha=0.5)

# Muestra del resultado
plt.show()