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

# # Exploration of tessellation
# 
# In search of an optimal spatial unit, this notebook explores several options of tessellation.

# In[1]:


import os

import geopandas as gpd

from sqlalchemy import create_engine

user = os.environ.get('DB_USER')
pwd = os.environ.get('DB_PWD')
host = os.environ.get('DB_HOST')
port = os.environ.get('DB_PORT')


# In[2]:


db_connection_url = f"postgres+psycopg2://{user}:{pwd}@{host}:{port}/built_env"
engine = create_engine(db_connection_url)
engine.begin()


# In[10]:


x, y = 352125.32, 492802.86  # coordinates in epsg 27700
buffer = 5000  # radius in [m]


# In[11]:


sql = f'SELECT * FROM openroads_200803_topological WHERE ST_DWithin(geometry, ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer})'

roads = gpd.read_postgis(sql, engine, geom_col='geometry')


# In[12]:


roads.plot()


# In[13]:


sql = f'SELECT * FROM openmap_buildings_200814 WHERE ST_DWithin(geometry, ST_SetSRID(ST_Point({x}, {y}), 27700), {buffer})'

buildings = gpd.read_postgis(sql, engine, geom_col='geometry')
buildings.plot()


# In[14]:


buildings.shape


# In[15]:


from shapely.geometry import Point


# In[16]:


limit = Point(x, y).buffer(buffer)


# In[17]:


import momepy as mm


# In[19]:


buildings


# In[28]:


roads['uID'] = range(len(roads))


# In[29]:


tes_rd = mm.Tessellation(roads, 'uID', limit, segment=10)


# In[30]:


roads_tessellation = tes_rd.tessellation


# In[38]:


ax = roads_tessellation.plot(figsize=(12, 12))
roads.plot(ax=ax, color='k', linewidth=.2)


# In[33]:


buildings['uID'] = range(len(buildings))


# In[34]:


tes_blg = mm.Tessellation(buildings, 'uID', limit, segment=2, shrink=0)


# In[35]:


blg_tessellation = tes_blg.tessellation


# In[37]:


ax = blg_tessellation.plot(figsize=(16, 16))
buildings.plot(ax=ax, color='k')


# In[45]:


union = buildings.buffer(100).unary_union


# In[55]:


ind = roads.sindex.query(union, predicate='contains')


# In[65]:


roads[~roads.index.isin(ind)].plot()


# In[66]:


country_roads = roads[~roads.index.isin(ind)].copy()


# In[67]:


country_roads.uID = country_roads.uID + 10000


# In[69]:


mixed = country_roads[['uID', 'geometry']].append(buildings[['uID', 'geometry']])


# In[70]:


tes_mix = mm.Tessellation(mixed, 'uID', limit, segment=2, shrink=0)


# In[72]:


mixed_tessellation = tes_mix.tessellation
ax = mixed_tessellation.plot(figsize=(16, 16))
mixed.plot(ax=ax, color='k')


# In[75]:


total = buildings.total_bounds


# In[76]:


import numpy as np


# In[78]:


pts = []

for x in np.linspace(total[0], total[2], 100):
    for y in np.linspace(total[1], total[3], 100):
        pts.append(Point(x, y))


# In[79]:


pts = gpd.GeoSeries(pts)


# In[84]:


grid = pts.drop(pts.sindex.query(union, predicate='intersects'))


# In[88]:


ax = grid.plot(markersize=.1, figsize=(12, 12))
buildings.plot(ax=ax)


# In[94]:


grid = gpd.GeoDataFrame(geometry=grid, crs=buildings.crs)


# In[95]:


grid['uID'] = range(len(grid))
grid['uID'] = grid['uID'] + 100000


# In[96]:


# fix for momepy
grid.geometry = grid.buffer(.1)


# In[103]:


buildings_grid = buildings[['uID', 'geometry']].append(grid).reset_index(drop=True)


# In[104]:


buildings_grid


# In[107]:


tes_grid = mm.Tessellation(gpd.clip(buildings_grid, limit), 'uID', limit, segment=2, shrink=0)


# In[113]:


grid_tessellation = tes_grid.tessellation
ax = grid_tessellation.plot(figsize=(20, 20), edgecolor='w', linewidth=.2)
buildings_grid.plot(ax=ax, color='k')


# In[115]:


buildings_grid.to_file('units.gpkg', layer='buildings_grid', driver='GPKG')
grid_tessellation.to_file('units.gpkg', layer='grid_tessellation', driver='GPKG')
buildings.to_file('units.gpkg', layer='buildings', driver='GPKG')
roads.to_file('units.gpkg', layer='roads', driver='GPKG')
mixed.to_file('units.gpkg', layer='mixed', driver='GPKG')
blg_tessellation.to_file('units.gpkg', layer='blg_tessellation', driver='GPKG')
roads_tessellation.to_file('units.gpkg', layer='roads_tessellation', driver='GPKG')
mixed_tessellation.to_file('units.gpkg', layer='mixed_tessellation', driver='GPKG')


# In[ ]: