#!/usr/bin/env python
# coding: utf-8

# # Los 40,000$m^2$ mas llenos de edificios en España

# In[14]:


get_ipython().run_line_magic('matplotlib', 'inline')

import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

google = "http://www.google.cn/maps/vt?lyrs=s@189&gl=cn&x=tileX&y=tileY&z=tileZ"


# Primero, vamos a cargar el grid con celdas donde hay *algún* edificio:

# In[2]:


db = gpd.read_file('../../01 AUDES/grid/grid.shp')


# In[3]:


db.head()


# La celda con más edificios puede ser identificada seleccionando la fila con el mayor valor en la columna `N. Buildings`:

# In[4]:


densest = db.loc[db['N. Buildin'] == db['N. Buildin'].max(), ]
densest


# In[76]:


def plot_cell(densest, figsize=(9, 9), ax=None):
    if ax is None:
        f, ax = plt.subplots(1, figsize=figsize)
    densest.to_crs(epsg=3857)\
           .plot(facecolor='none', edgecolor='red', linewidth=5, ax=ax)
    ctx.add_basemap(ax, zoom=18, url=google, 
                    attribution='Imagery by Google Tile Map Service')
    minX, maxX, minY, maxY = ax.axis()
    p = densest.to_crs(epsg=3857).iloc[0].geometry
    shade = Polygon(shell=((minX, minY), (maxX, minY), 
                           (maxX, maxY), (minX, maxY), 
                           (minX, minY)),
                    holes=[tuple(p.exterior.coords)])
    gpd.GeoSeries(shade, crs={'init': 'epsg:3857'}).plot(ax=ax, color='k', alpha=0.5)
    ax.set_axis_off()
    return None

plot_cell(densest)


# Usando la misma lógica, podemos tambien explorar el top 9 (encaja mejor en una matriz 3x3):

# In[42]:


top9 = db.sort_values('N. Buildin', ascending=False)\
         .head(10)
top9


# Y las imagenes en un grid:

# In[91]:


f, axs = plt.subplots(3, 3, figsize=(12, 12))
axs = axs.flatten()
for i in range(axs.shape[0]):
    plot_cell(top9.iloc[[i], :], ax=axs[i])
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('top9.png', dpi=300)
plt.show()


# In[92]:


get_ipython().system(' du -h top9.png')