<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  Como pueden ver, el archivo tiene 42 columnas, las que nos interesan en este momento son las que describen la actividad de cada unidad
Step2:  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
Step3:  Lo primero que hicimos fue importar dos librerías
Step4:  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
Step5:  Ejercicio rápido
Step6:  La función parece bastante complicada, pero si lo pensamos con clama en realidad es bastante sencilla
Step7:  Ejercicio rápido
Step8:  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
Step9:  Ejercicio rápido
Step10:  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
Step11:  Entonces, ya sabemos como calcular nuestro identificador único de AGEB, ahora agreguémoslo como columna
Step12:  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
Step13:  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
Step14:  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
Step15:  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
Step16:  Como sólo vamos a utilizar las coumnas de geometría y el identificador de AGEB, guardemos el GeoDataFrame con sólo esas columnas
Step17:  Ahora vamos a unir las dos bases de datos
Step18:  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
Step19:  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)
<ASSISTANT_TASK:>
Python Code:
import geopandas as gpd
denue = gpd.read_file("datos/DENUE_INEGI_09_.shp")
denue.head()
Explanation: 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:
$$V_{t}(i,j) = [(i-1,j), (i+1,j), (i,j+1), (i,j-1)]$$
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:
Leer los datos de dos fuentes diferentes y procesarlos para calcular las variables que nos interesan
Encontrar los vecinos de cada polígono
Calcular un par de variables en función de una vecindad de primer orden
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:
End of explanation
denue[["codigo_act", "nom_estab"]].head()
Explanation: 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 [ ])
End of explanation
import pandas as pd
import numpy as np
df = pd.DataFrame(np.random.randn(3, 3))
df
Explanation: 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:
End of explanation
def square(x):
    return x**2
df[[0]].apply(square)
Explanation: 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:
End of explanation
df['squared'] = df[[0]].apply(square)
df
Explanation: 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:
End of explanation
def clasifica(codigo):
    comercio = ['461', '462', '463', '464', '465', '466']
    oficinas = ['51', '521', '523', '524', '5312', '5313', '541', '55']
    ocio = ['711121', '71212', '7132', '7139', '7211', '7224', '7225']
    usos = {'comercio': comercio, 'oficinas':oficinas, 'ocio': ocio}
    for actividad, claves in usos.items():
        for c in claves:
            if str(codigo).startswith(c):
                return actividad
            
Explanation: 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:
End of explanation
clases = denue['codigo_act'].apply(clasifica)
denue['clase'] = clases
denue[['codigo_act', 'clase']][denue['clase']== 'oficinas'].head()
Explanation: 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:
End of explanation
denue.loc[denue['clase'].isnull()].head()
Explanation: 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:
End of explanation
denue = denue.loc[denue['clase'].notnull()]
denue.head()
Explanation: 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:
End of explanation
def concatena_claves(x):
    return '{}{}{}{}'.format (x['cve_ent'], x['cve_mun'], x['cve_loc'], x['ageb'])
denue.apply(concatena_claves, axis=1).head()
Explanation: 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:
End of explanation
print('{} {}      {}'.format('a', 'b', 'c'))
Explanation: 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:
End of explanation
denue['cve_geo'] = denue.apply(concatena_claves, axis=1)
denue[['cve_geo','clase']].head()
Explanation: Entonces, ya sabemos como calcular nuestro identificador único de AGEB, ahora agreguémoslo como columna:
End of explanation
variables = denue[['cve_geo','clase']]
variables.head()
Explanation: 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):
End of explanation
pd.get_dummies(variables, columns=['clase']).head(20)
Explanation: 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:
End of explanation
variables = pd.get_dummies(variables, columns=['clase'])
por_ageb = variables.groupby(['cve_geo']).sum()
por_ageb.head()
Explanation: 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:
End of explanation
agebs = gpd.read_file("datos/ageb_urb.shp")
agebs.head()
Explanation: 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:
End of explanation
agebs = agebs[['CVEGEO', 'geometry']]
agebs.head()
Explanation: Como sólo vamos a utilizar las coumnas de geometría y el identificador de AGEB, guardemos el GeoDataFrame con sólo esas columnas:
End of explanation
usos_suelo = agebs.merge(por_ageb, left_on='CVEGEO', right_index=True, how='inner')
usos_suelo.head()
Explanation: 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:
End of explanation
%matplotlib inline
usos_suelo.plot(column='clase_comercio', figsize=(20,10))
Explanation: 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
End of explanation
usos_suelo.plot(column='clase_comercio', figsize=(20,10), cmap='OrRd',scheme='quantiles')
usos_suelo.plot(column='clase_comercio', figsize=(20,10), cmap='OrRd',scheme='fisher_jenks')
Explanation: 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):
End of explanation
<END_TASK>