import pandas as pd
import geopandas as gpd
from preproceso import preprocesa
pd.options.mode.chained_assignment = None

denue = gpd.read_file("datos/DENUE_INEGI_09_.shp")
agebs = gpd.read_file("datos/ageb_urb.shp")
usos_suelo = preprocesa(denue, agebs)
usos_suelo.reset_index(drop=True, inplace=True)
usos_suelo.head()

"""
Explanation: Parte 2. Cálculos en vecindades
En la primera parte del taller trabajamos los datos para calcular las variables que nos interesan y agruparlas por AGEB. Ahora, en esta segunda parte, vamos a partir de los datos ya procesados y vamos a realizar algunos cálculos a partir de los valores en los polígonos vecinos.
El primer paso es encapsular  el preproceso de los datos en una función, de esta manera podemos ejecutar todo el flujo en una celda en este Notebook. En el archivo preproceso.py de esta carpeta pueden ver la función que realiza todo el trabajo.
Para importar la función aquí, hacemos uso de la estructura de paquetes (packages) de Python: siempre que una carpeta contenga un archivo con el nombre __init__.py, esa carpeta se convierte automáticamente en un paquete y permite importar directamente todas las funciones (o clases) definidas en los archivos que contiene.
En la siguiente celda realizamos todo el trabajo de preprocesamiento.
Nota 1: la declaración pd.options.mode.chained_assignment = None sirve para suprimir una advertencia de Pandas sobre la forma en la que asignamos valores a nuevas columnas, en este caso (y en muchos otros), esa advertencia es espúrea, sin embargo, no está de más que lean esto.
Nota 2: si se fijan, en el preproceso removimos algunas filas, por lo que el índice del DataFrame no es consecutivo, esto va a ser importante más adelante, por lo que nos conviene resetear el índice para que se vuelva a calcular como un consecutivo.
End of explanation
"""


import pysal

w = pysal.weights.Queen.from_dataframe(usos_suelo)
print(w.n)
print(w.weights[0])
print(w.neighbors[0])
print(w.neighbors[5])
print(w.histogram)

"""
Explanation: Como pueden ver, en la variable usos_suelo tenemos ya calculadas todas nuestras variables de interés, ahora lo que necesitamos es, para cada fila de nuestro GeoDataFrame, saber cuáles son los polígnos vecinos.
Para esto, vamos a utilizar la librería PySal, que provee un conjunto de métodos de análisis espacial. En particular, nos interesa la funcionalidad de crear matrices de pesos espaciales.
PySal está desarrollado para trabajar en conjunto con GeoPandas, de modo que podemos pedir la matriz de pesos directamente del GeoDataFrame y examinar el objeto que nos regresa:
End of explanation
"""


freqs = pd.DataFrame(w.histogram, columns=['vecinos', 'cuenta'])
freqs.head()

"""
Explanation: Lo primero que hicimos fue importar la librería PySal. A continuación, claculamos la matriz de pesos w usando vecindades de tipo Reina (en la documentación de PySal pueden consultar los diferentes tipos de vecindades y las fuentes de datos que pueden usar).

w.n nos dice la cantidad de renglones de la matriz
w.weights[0] nos dice los pesos que corresponden a los vecinos del elemento 0
w.neighbors[0] nos da la lista de vecinos del elemento 0
w.histogram nos da el histograma de la matriz de adyacencia, es decir, cuántos elementos tienen x número de vecinos

Como un ejercicio rápido vamos a graficar el histograma, sólo que esta vez, en lugar de usar matplotlib directamente, vamos a usar seaborn, que es una librería para graficar datos estadísticos. Además de producir, de manera sencilla, graficas más bonitas que matplotlib, seaborn tiene una sintaxis similar a la de ggplot2 de R.
Primero convertimos el histograma que nos da PySal en un DataFrame:
End of explanation
"""


%matplotlib inline
import seaborn as sns

sns.barplot(x='vecinos', y='cuenta', data=freqs)

"""
Explanation: Y luego lo graficamos:
End of explanation
"""


usos_suelo = usos_suelo.drop(usos_suelo.index[[1224]])
usos_suelo.reset_index(drop=True, inplace=True)
w = pysal.weights.Queen.from_dataframe(usos_suelo)

"""
Explanation: Intensidad
Después de este intermedio, ahora sí vamos a hacer nuestro primer cómputo en vecindades. Vamos a comenzar por la intensidad.
La intensidad es simplemente la cantidad de actividades en un área determinada. En nuestro caso, vamos a calcular el total de actividades (de cualquier tipo) que hay en la vecindad inmediata de cada AGEB (si lo piensan un poco, esto se parece bastante a los filtros tipo blur en procesamiento de imágenes).
Para calcular la intensidad, lo que necesitamos hacer es recorrer la lista de elementos del GeoDataFrame y, para cada elemento, obtener la lista de vecinos, sacar sus variables y sumarlas.
Antes de calcular, vamos a eliminar el elemento que no tiene ningún vecino, reindexar los datos y volver a calcular los pesos (para que los índices de la matriz de pesos y del DataFrame coincidan):
End of explanation
"""


usos_suelo.iloc[[0]][['clase_comercio', 'clase_ocio', 'clase_oficinas']].as_matrix()

import numpy as np

intensidad =[]
for i in range(0, w.n):
    vecinos = w.neighbors[i]
    total = 0.0
    suma = np.zeros((3),dtype=np.float)
    valores = usos_suelo.iloc[[i]][['clase_comercio', 'clase_ocio', 'clase_oficinas']].as_matrix()
    for j in range(0,len(vecinos)):    
        data = usos_suelo.iloc[[j]][['clase_comercio', 'clase_ocio', 'clase_oficinas']].as_matrix()
        suma = suma + data
        total += sum(data)
    intensidad.append((i, sum(total)))

print(intensidad[0:10])

"""
Explanation: Ahora sí, recorremos la lista de vecinos y calculamos la intensidad para cada elemento:
End of explanation
"""


intensidad_df = pd.DataFrame(intensidad, columns=['gid', 'intensidad'])
datos_intensidad = usos_suelo.merge(intensidad_df, left_index=True, right_on='gid', how='inner')
datos_intensidad.head()

"""
Explanation: Al parecer lo que estamos haciendo es muy complicado, sin embargo, una vez más, si lo vemos con detenimiento es relativamente simple:

Primero estamos definiendo una lista vacía intensidad que nos va a servir para guardar los resultados
Luego, en el for externo, estamos recorriendo la matriz de adyacencia, entonces el índice del for, es el identificador de cada polígono
Inicializamos un array de 3 entradas con ceros (esto nos va a servir para guardar la suma para cada uso de suelo)
Con iloc tomamos la fila correspondiente en el DataFrame y as_matrix() convierte los valores de las columnas en un array
Recorremos en el for interno los vecinos de cada elemento y tomamos, como array, sus valores
Sumamos los arrays entrada por entrada (esto realmente no es necesario aquí, pero va a ser útil más adelante cuando hagamos un cálculo más complejo)
A la salida de los dos for, agregamos a la lista intensidad una tupla con el índice y el valor de la intensidad

Entonces, podemos convertir la lista intensidad en un DataFrame para después unirlo con nuestros datos:
End of explanation
"""


entropia =[]
for i in range(0, w.n):
    vecinos = w.neighbors[i]
    total = 0.0
    suma = np.zeros((3),dtype=np.float)
    valores = usos_suelo.iloc[[i]][['clase_comercio', 'clase_ocio', 'clase_oficinas']].as_matrix()
    for j in range(0,len(vecinos)):    
        data = usos_suelo.iloc[[j]][['clase_comercio', 'clase_ocio', 'clase_oficinas']].as_matrix()
        suma = suma + data
        total += np.sum(data)
    p = np.nan_to_num(suma/total)
    lp = np.select([p == 0,p > 0],[p, np.log(p)])
    entropia.append((i, np.sum(p*lp)))

print(entropia[0:10])

"""
Explanation: Ejercicio
Hagan un mapa que destaque las diferencias en intensidad.
Entropía
La entropía es una medida de la mezcla de usos de suelo, está basada en la forma en la que se calcula la entropía en mecánica estadística:
$$ E = \sum\limits_{j}{\frac{p_{j}*ln(p_{j})}{ln(J)}} $$
Donde $p_{j}$ representa la proporción del $j-ésimo$ uso de suelo con respecto al total y $J$ es el número de usos de suelo considerados. Valores cercanos a 0 indican poca mezcla de usos de suelo y valores cercanos a -1 indican una mezcla balanceada.
Entonces, para calcular la entropía, basta con modificar un poco el for que usamos para calcular la intensidad:
End of explanation
"""


