Cálculos en vecindades I

En este taller vamos a utilizar Python para calcular un par de variables, asociadas a cada polígono del espacio, en función de las propiedades de los polígonos vecinos.

Si pensamos en términos de rasters el ejercicio es relativamente simple ya que la función de vecindad está implícita en la estructura de los datos, por ejemplo, si tenemos el elemento $(i,j)$ de una matriz, entonces su vecindad de torre $V_{t}(i,j)$ está dada por:

De foma análoga es posible encontrar la vecindad de reina de cada elemento de la matriz.

Ahora, en el caso de coberturas (en este contexto coberturas quiere decir que los polígonos cubren totalmente el espacio de estudio), con polígonos arbitrarios la cuestión es más compleja: no tenemos un índice que defina de forma implícita las relaciones de vecindad, tenemos que encontrar a los vecinos de cada polígono a partir de operaciones geométricas. Una de las ventajas de usar Python es que, al ser muy usado por la comunidad científica en general, existen muchas librerías que nos pueden ayudar a resolver nuestros problemas, en este caso pysal provee una manera fácil de acceder a los polígonos vecinos.

El taller está organizado de la siguiente forma:

1. Leer los datos de dos fuentes diferentes y procesarlos para calcular las variables que nos interesan
2. Encontrar los vecinos de cada polígono
3. Calcular un par de variables en función de una vecindad de primer orden
4. Representar el resultado que obtengamos

1. Leyendo y transformando los datos

Lo que vamos a calcular al final del taller son dos variables que representan la mezcla de usos de suelo y la intensidad de las actividades en una ciudad, estas son dos variables muy usdas en estudios urbanos para, por ejemplo, modelar el uso del automóvil como función de los entornos urbanos locales.

Para calcular estas variables, lo primero que necesitamos es obtener un conteo de los diferentes tipos de actividades urbanas que ocurren en una zona determinada. Para esto, vamos a utilizar la base de datos del DENUE publicada por el INEGI. Duarante el taller trabajaremos con datos para la Ciudad de México, sin embargo, para la tarea ustedes trabajarán con una zona metropolitana de su elección. Para agregar los datos usaremos las Áreas Geoestadísticas Básicas (AGEBS) del INEGI.

Entonces, el primer paso es leer los datos del DENUE, para eso vamos a utilizar la librería GeoPandas que nos provee una estructura de datos similar a los DataFrames de R y que, además, soporta objetos geográficos.

Para facilitar el trabajo, descompriman los archivos en el directorio datos de este repositorio.

Entonces, para leer los datos, lo primero que necesitamos hacer es importar la librería GeoPandas y después crear un GeoDataFrame a partir del shapefile con los datos:

Como pueden ver, el archivo tiene 42 columnas, las que nos interesan en este momento son las que describen la actividad de cada unidad:

(Noten como para pedir las columnas que queremos, simplemente pasamos la lista de nombres al operador de selección [ ])

El código de la actividad corresponde a la clasificación del SCIAN y sirve para clasificar las actividades económicas. Para este taller, vamos a agregar diferentes actividades económicas en tres grupos de usos de suelo:

• Comercio: actividades que se dedican a la venta al menudeo
• Oficinas: lugares en los que se trabaja en una oficina
• Ocio: sitios de esparcimiento. Bares, teatros, etcétera

Cada uno de estos usos es una agregación de diferentes actividades del SCIAN, para este taller vamos a considerar las siguientes agregaciones:

• Comercio: empiezan con 461, 462, 463, 464, 465 o 466
• Oficinas: empiezan con 51, 521, 523, 524, 5312, 5313, 541 o 55
• Ocio: empiezan con 711121, 71212, 7132, 7139, 7211, 7224 o 7225

Para etiquetar cada renglón en la base de datos con el uso al que pertenece vamos a usar el método apply de GeoPandas (en realidad es un método de la librería Pandas que es heredado por su extensión espacial).

Lo que hace apply es aplicar una función a cada elemento de una Serie. Una Serie es una columna de un DataFrame. Para entender cómo funciona, vamos a revisar el ejemplo de la documentación:

Lo primero que hicimos fue importar dos librerías: Pandas y Numpy (en este caso numpy sólo la vamos a utilizar para generar un arreglo de números aleatorios) y después construir un DataFrame con 3 renglones y 3 columnas.

Ahora vamos a aplicar una función a cada elemento de la primera columna del DataFrame:

Como pueden ver, el resultado es una serie (una columna) con el resultado de aplicar nuestra función a cada renglón de la serie 0. Esta nueva columna la podemos agregar fácilmente a nuestros datos originales:

Ejercicio rápido: repite lo que acabamos de hacer pero con una función que calcule la mitad de cada elemento.

Ahora sí podemos regresar a nuestros propios datos. Lo primero que tenemos que hacer es escribir una función que, a partir del código del SCIAN, nos regrese nuestra propia clasificación:

La función parece bastante complicada, pero si lo pensamos con clama en realidad es bastante sencilla:

• Primero hacemos una lista con los caractéres que definen cada actividad
• Luego hacemos un diccionario que relaciona el nombre de la actividad con lla listya de caracteres que la definen
• Iteramos sobre el diccionario, es decir, para cada elemento del diccionario obtenemos el par llave:valor
• Como en este caso el valor es una lista, entonces para cada elemento de las lista de caracteres buscamos si el código de la actividad empieza con los caracteres que nos interesan y si es cierto, regresamos el nombre de la actividad

Apliquemos ahora la función a nuestros datos:

Ejercicio rápido: Verifica las demás clases de actividades.

Ahora, si se dan cuenta, la clasificación que hicimos no regresa valores para todas la filas. Es decir, no todas las actividades del DENUE corresponden a alguna de nuestras categorías:

Fíjense como estamos usando aquí loc para seleccionar filas que cumplan con una condición. Pandas tiene varias formas de seleccionar (ya hemos usado por lo menos dos aquí). Aquí puedes ver la documentación oficial de los métodos para seleccionar filas en Pandas, en esta otra liga hay algunos ejemplos extra.

Filtremos entonces las filas que no corresponden a alguna de nuestras clases:

Ejercicio rápido: Comprueba que, en efecto, denue ya no contenga valores nulos en la columna clase.

Recapitulemos un momento. Hasta aquí tenemos un GeoDataFrame con todos los puntos del DENUE que corresponden a una de las categorías de uso de suelo que nos interesan. Como dijimos al principio, queremos calcular las mezclas de usos de suelo al nivel AGEB, entonces, todavía necesitamos agregar nuestras variables por AGEB.

Una forma de agregar nuestros datos por AGEB es hacer una unión espacial de los puntos del DENUE con las AGEBS, de modo que a cada punto se le asigne el identificador de la AGEB en la que está. En este caso, vamos a tomar una aproximación diferente, basada en la construcción de un identificador único para cada AGEB. Si abren en QGIS el shape de las AGEBs, van a notar que tiene 3 columnas que, en conjunto, identifican de forma única a cada AGEB: cve_ent, cve_mun, cve_loc y cve_ageb. Las mismas columnas aparecen en el DENUE, de modo que, si concatenamos estas 3 columnas en ambas bases de datos, tendremos un identificador que nos permite realizar la unión sin necesidad de calcular una unión espacial.

Empecemos por calcular el identificador de AGEB en el DENUE. Intuitivamente, lo que necesitamos es algo similar al apply que utilizamos para clasificar las actividades, la diferencia está en que en aquel caso, la función se aplicaba sobre una sola columna y ahora la función debe tomar varias columnas como entrada. Pandas provee una función apply que opera sobre todo un DataFrame, ya sea sobre las columnas o sobre las filas. En este caso, vamos a utilizar la versión que opera sobre las columnas, de esta forma podremos concatenar las columnas que nos interesan:

Fíjense que aquí estamos introduciendo el operador de formato de caracteres para rellenar un string con los valores que obtenemos de las columnas. Aquí hay un ejemplo rápido para que vean cómo funciona format:

Entonces, ya sabemos como calcular nuestro identificador único de AGEB, ahora agreguémoslo como columna:

Ahora sí, ya que tenemos un identificador de AGEB para cada punto, podemos agregarlos por AGEB y calcular cuantas ocurrencias de cada una de nuestras categorías hay en cada AGEB, como lo haríamos con un group by de SQL. El problema en este caso es que las clases son valores a lo largo de una única columna, de modo que no tenemos una forma directa de agregar la cantidad de ocurrencias de cada uso de suelo. Entonces, antes de agregar por AGEB es necesario promover los valores de cada clase a columnas, es decir, agregar 3 columnas a nuestra base de datos: una por cada clase y llenar las filas con unos o ceros de acuerdo a la clase a la que corresponde el punto.

Primero vamos a seleccionar y guardar en un DataFrame las columnas que nos interesan (para mantener las cosas simples):

Fíjense que lo que queremos hacer es codificar la variable clase como tres variables dummies, es decir convertir la columna en tres columnas indicadoras. Pandas nos provee el método get_dummies para hacer esto:

Ahora tenemos una columna por cada clase que nos indica si el punto es o no de dicha clase. Entonces ya sólo es cuestión de agrupar por nuestro identificador de AGEB y sumar las columnas indicadoras:

Ahora sí, finalmente tenemos la cantidad de ocurrencias de cada uso de suelo para cada AGEB, Ahora sólo necesitamos unirlas con la geometría de las AGEBS y tendremos nuestros datos listos!

Vamos entonces a leer el archivo con las geometrías de las AGEBS:

Como sólo vamos a utilizar las coumnas de geometría y el identificador de AGEB, guardemos el GeoDataFrame con sólo esas columnas:

Ahora vamos a unir las dos bases de datos: agebs y por_ageb. En SQL esta operación es un join usando las llaves, en Pandas (o en R), esta operación se llama merge:

Fíjense que el merge funciona exáctamente igual que un join en SQL, la única peculiaridad en este caso es el uso de right_index=True. Esto es necesario porque, por la forma en la que hicimos el DataFrame por_ageb, la clave de AGEB es el índice de la tabla y no una columna (de alguna forma, el índice es el nombre de cada fila, se parece bastante a una columna pero no tienen nombre).

Para terminar esta parte del taller, vamos a hacer un mapa coloreado por la cantidad de comercios:

Nota: la línea %matplotlib inline es un comando de IPython para indicar que la gráfica debe imprimirse directamente en la página

Como se ve, el mapa no muestra grandes diferencias, eso es porque por defecto, GeoPandas clasifica en intervalos iguales. Si queremos destacar las diferencias, podemos usar otros métodos de clasificación (esta opción sólo está disponible si PySal está también instalado):