In [1]:
import os
import requests
import numpy as np
import pandas as pd
import geopandas as gpd
import h3pandas
import googlemaps
import time
import fiona
import urllib
import folium

from tqdm import tqdm
from glob import glob
from datetime import datetime
from class_hexagons import hexagonsClass
from dotenv import load_dotenv,find_dotenv

c = hexagonsClass()
In [2]:
load_dotenv(find_dotenv())
Out[2]:
True

Carico il token di Mapbox

In [3]:
access_token = os.getenv('MAPBOX_API_KEY')

Imposto alcuni parametri

In [4]:
resolution = 8
selection = 'Sardegna'

Carico il layer dei comuni ISTAT che mi servirà come riferimento territoriale

In [5]:
comuni = gpd.read_parquet(os.path.join(c.work_path,'ISTAT','DATA','comuni_pop.parquet'))
comuni.to_crs(4326,inplace=True)

Produco la griglia H3 alla risoluzione indicata all'interno della regione (Sardegna), con un buffer di 500 metri.

In [6]:
# h3_index = f'h3_{resolution:02}'
# selection_type = 'regione'
# start = datetime.now()

# selection_df = comuni[comuni[selection_type]==selection].copy()
# selection_df = selection_df.dissolve()
# selection_df['geometry'] = selection_df['geometry'].to_crs(3857).buffer(500).to_crs(4326)

# grid = c.get_grid_of_selection(selection_df,resolution)
In [7]:
# lat = grid[grid[f'h3_{resolution:02}']=='881e90c5e5fffff']['geometry'].to_crs(3857).centroid.to_crs(4326).y.values[0]
# lon = grid[grid[f'h3_{resolution:02}']=='881e90c5e5fffff']['geometry'].to_crs(3857).centroid.to_crs(4326).x.values[0]
In [8]:
#grid = grid.sjoin(selection_df, how='inner', predicate='within' )
In [9]:
# grid.to_parquet(os.path.join(c.data_path,f'grid_{selection.lower()}_H3_{resolution:02}.parquet'))
# print(f'Esecuzione: {str(datetime.now()-start).split(".")[0]}')

Carico la griglia prodotta

In [10]:
grid = gpd.read_parquet(os.path.join(c.task_path, c.config['data_dir'], f'grid_{selection.lower()}_H3_{resolution:02}.parquet'))
In [11]:
grid.shape[0]
Out[11]:
31602

Imposto i parametri delle isocrone da produrre (7 e 15 minuti)
Attenzione il parametro "depart_at" sembra essere a pagamento, o comunque ci vuole la carta caricata nell'account

In [12]:
profiles = [
    'mapbox/driving-traffic',
    'mapbox/driving',
    'mapbox/walking',
    'mapbox/cycling'
]
api_url = 'https://api.mapbox.com'
service = 'isochrone'
version = 'v1'
profile = profiles[0]
minutes = '7,15'
depart_at = 'depart_at=2024-10-03T08:00:00Z&'

Seleziono un comune da analizzare

In [13]:
#comuni = comuni[comuni['comune']=='Roma'].copy()

Identifico le celle contenute nel comune

In [14]:
#grid = gpd.sjoin(grid, comuni, how='inner', predicate='within')

Estraggo i centroidi delle celle selezionate

In [15]:
grid_p = grid.to_crs(3857)
grid_p['geometry'] = grid_p['geometry'].centroid
grid_p.to_crs(4326,inplace=True)

Visualizzazione per controlli delle celle in mezzo al mare (ora risolto)

In [16]:
# lat = grid[grid[f'h3_{resolution:02}']=='881e90c5e5fffff']['geometry'].to_crs(3857).centroid.to_crs(4326).y.values[0]
# lon = grid[grid[f'h3_{resolution:02}']=='881e90c5e5fffff']['geometry'].to_crs(3857).centroid.to_crs(4326).x.values[0]

# m = grid[grid[f'h3_{resolution:02}']=='881e90c5e5fffff'].explore()
# m = grid_p[grid_p[f'h3_{resolution:02}']=='881e90c5e5fffff'].explore(m=m)
# url =f'{api_url}/{service}/{version}/{profile}/{lon},{lat}?contours_minutes={minutes}&polygons=true&access_token={access_token}'
# gdf = gpd.read_file(url)
# m = gdf.explore(m=m) 
# m

Produco le isocrone per ogni centroide estratto

In [17]:
# # 6min - Roma
# # 2:47:22 - Sardegna
# start = datetime.now()
# iso_list = list()
# errors = list()
# grid_dicts = grid_p.to_dict('records')
# for cell in tqdm(grid_dicts, total=len(grid_dicts)):
#     for attempt in range(10):
#         try:
#             lon = cell['geometry'].x
#             lat = cell['geometry'].y
#             url =f'{api_url}/{service}/{version}/{profile}/{lon},{lat}?contours_minutes={minutes}&polygons=true&access_token={access_token}'
#             gdf = gpd.read_file(url)    
#         except urllib.error.HTTPError as err:
#             print(f'tentativo = {attempt} - {err}')
#             time.sleep(20)
#         except fiona.errors.DriverError as err:
#             print(f'Location out of road - {err}')
#             errors.append(cell['geometry'])
#             df = pd.DataFrame([{'geometry': None}])
#             gdf = gpd.GeoDataFrame(df, geometry=df['geometry'])
#             break
#         else:
#             break
#     else:
#         # entra qui se ha fallito tutti i tentativi
#         print('Errore non risolvibile')
#         break
        
#     gdf[f'h3_{resolution:02}'] = cell[f'h3_{resolution:02}']
#     iso_list.append(gdf)


# iso = pd.concat(iso_list)
# iso_et = str(datetime.now()-start).split(".")[0]
# print(f'Esecuzione: {iso_et}')

Salvo le isocrone prodotte

In [18]:
# iso.to_parquet(os.path.join(c.data_path,f'iso_{selection.lower()}_H3_{resolution:02}.parquet'))

Carico le isocrone prodotte

In [19]:
iso = gpd.read_parquet(os.path.join(c.data_path,f'iso_{selection.lower()}_H3_{resolution:02}.parquet'))
In [20]:
#iso[iso[f'h3_{resolution:02}']==cell_sel].head()

Eseguo il join con la griglia RISTRETTA per eliminare le celle nel mare, a regime la griglia sarà direttamente clippata sulla regione

In [21]:
iso = iso.merge(grid[[f'h3_{resolution:02}']], on=f'h3_{resolution:02}')

Eseguo l'intersect tra le isocrone prodotte e i servizi

In [22]:
osm_pois = gpd.read_parquet(os.path.join(c.work_path,'OSM','GEOP','OSM_POIS.parquet'))
osm_pois.drop_duplicates(inplace=True)
In [23]:
iso = iso[['contour','metric',f'h3_{resolution:02}','geometry']].copy()
osm_pois = osm_pois[['osm_id','fclass','geometry']].copy()
In [24]:
start = datetime.now()
i = osm_pois.sjoin(iso, predicate='intersects')
int_et = str(datetime.now()-start).split(".")[0]
print(f'Esecuzione: {int_et}')

Esecuzione: 0:02:34

Salvo lo spatial join tra pois e iso

In [25]:
#i.to_parquet(os.path.join(c.data_path,f'iso_pois_{selection.lower()}_H3_{resolution:02}.parquet'))

Rileggo l'intersect

In [26]:
i_int = gpd.read_parquet(os.path.join(c.data_path,f'iso_pois_{selection.lower()}_H3_{resolution:02}.parquet'))

Stats

In [27]:
errors = iso[iso['geometry'].isnull()].shape[0]
In [28]:
print(f'Numero di esagoni alla risoluzione {resolution} per la regione {selection}: {grid.shape[0]}')
print(f'Profilo: {profiles[0]}')
print(f'Tempi di percorrenza isocrone (minuti): {" - ".join(minutes.split(","))}')
print(f'Totale richieste: {grid.shape[0] * len(minutes.split(","))}')
print(f'Numero di isocrone prodotte ({iso.shape[0]}) + Location-out-of-road ({errors}) = {iso.shape[0] + errors}')
print('Tempo di esecuzione per le isocrone: 2:47:22')
print(f'Tempo di esecuzione per l\'intersect (poi vs iso): 0:02:14')

Numero di esagoni alla risoluzione 8 per la regione Sardegna: 31602
Profilo: mapbox/driving-traffic
Tempi di percorrenza isocrone (minuti): 7 - 15
Totale richieste: 63204
Numero di isocrone prodotte (63198) + Location-out-of-road (0) = 63198
Tempo di esecuzione per le isocrone: 2:47:22
Tempo di esecuzione per l'intersect (poi vs iso): 0:02:14

Calcolo i conteggi per le categorie di POI

In [29]:
g = i.groupby([f'h3_{resolution:02}','contour','fclass'], as_index=False).size()
In [30]:
#136
len(g['fclass'].unique())
Out[30]:
136

Inserisco un indirizzo da cui estrarre le coordinate per visualizzare un esempio del processo

In [31]:
#search_text = 'via filippi 16 roma rome, italy'
#search_text = 'via cornelio celso 11 roma rome, italy'
search_text = 'Via Cino da Pistoia, 20, 09128 Cagliari CA'

Geocoding con GoogleMaps

In [32]:
gmaps = googlemaps.Client(key=os.getenv('GOOGLE_API_KEY'))
resp = gmaps.geocode(search_text,language='it')
print(resp[0]['formatted_address'])
lat = resp[0]['geometry']['location']['lat']
lon = resp[0]['geometry']['location']['lng']
point = [{'id':1,'lat':lat,'lon':lon}]
point = pd.DataFrame(point)
point = gpd.GeoDataFrame(point, geometry=gpd.points_from_xy(point['lon'],point['lat']), crs=4326)

Via Cino da Pistoia, 20, 09128 Cagliari CA, Italia

Geocoding con Mapbox

In [33]:
# search_text = search_text.replace(' ','+')
# url = f'https://api.mapbox.com/search/geocode/v6/forward?q={search_text}&access_token={access_token}'

# point = gpd.read_file(url)
# point = point.head(1)

Identifico la cella di appartenenza del punto

In [34]:
cell_sel = gpd.sjoin(grid[[f'h3_{resolution:02}','geometry']], point, how='inner', predicate='intersects')
cell_sel = cell_sel[f'h3_{resolution:02}'].values[0]
In [35]:
g[g[f'h3_{resolution:02}']==cell_sel]
Out[35]:
h3_08 contour fclass size
61506 881e92992bfffff 7.0 archaeological 2
61507 881e92992bfffff 7.0 arts_centre 1
61508 881e92992bfffff 7.0 artwork 3
61509 881e92992bfffff 7.0 atm 2
61510 881e92992bfffff 7.0 attraction 4
... ... ... ... ...
61683 881e92992bfffff 15.0 wastewater_plant 9
61684 881e92992bfffff 15.0 water_tower 8
61685 881e92992bfffff 15.0 water_works 1
61686 881e92992bfffff 15.0 wayside_cross 4
61687 881e92992bfffff 15.0 wayside_shrine 2

182 rows × 4 columns

Visualizzo il risultato del processo per la cella che contiene il punto (iso a 7 minuti)

In [36]:
#cell_sel='881e90c5e5fffff'
In [37]:
m = iso[(iso[f'h3_{resolution:02}']==cell_sel)&(iso['contour']==15)].explore(color='color')
# m = point.explore(m=m, color='red')
m = grid[grid[f'h3_{resolution:02}']==cell_sel].explore(m=m, color='black',style_kwds={'fill':False})
m = grid_p[grid_p[f'h3_{resolution:02}']==cell_sel].explore(m=m, color='black')
m = i[(i[f'h3_{resolution:02}']==cell_sel)&(i['contour']==15)].explore(m=m, column='fclass')
folium.LayerControl().add_to(m)
m
Out[37]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [38]:
#m.save(os.path.join(c.work_path,'HTML','iso_poi_sel.html'))

Sezioni di censimento¶

Leggo il file con la selezione delle variabili di interesse

In [39]:
dfVar = pd.read_excel('../CSV/metadati_ISTAT_censimento.xlsx')
dfVar['NOME_CAMPO'] = dfVar['NOME_CAMPO'].str.lower()
dfVar = dfVar[dfVar['Sardegna'].notnull()].copy()
dfVar.shape
Out[39]:
(35, 3)
In [40]:
etichetteVariabili = {}
for codice, etichetta in zip(dfVar['NOME_CAMPO'].tolist(), dfVar['Sardegna'].tolist()):
    etichetteVariabili[codice] = etichetta

Leggo gli indicatori, seleziono solo la regione di interesse, rinomino le colonne degli indicatori con l'etichetta parlante

In [41]:
sezioni_21_ind = pd.read_parquet(os.path.join(c.work_path, 'ISTAT', 'DATA', 'sezioni_21_ind.parquet'))
sezioni_21_ind = sezioni_21_ind[sezioni_21_ind['regione']==selection].copy()
cols = ['sez21_id'] + list(etichetteVariabili)
sezioni_21_ind = sezioni_21_ind[cols].copy()
sezioni_21_ind = sezioni_21_ind.rename(columns=etichetteVariabili)
In [42]:
sezioni_21_ind.head()
Out[42]:
sez21_id POP_TOT POP_M POP_F M_5 M_5_9 M_10_14 M_15_19 M_20_24 M_25_29 ... F_30_34 F_35_39 F_40_44 F_45_49 F_50_54 F_55_59 F_60_64 F_65_69 F_70_74 F_74
439302 900010000001 293 137 156 3 1 4 12 5 3 ... 7 7 12 12 14 10 8 8 12 30
439303 900010000002 360 163 197 5 8 3 4 4 9 ... 10 9 12 11 17 18 16 16 12 34
439304 900010000003 453 219 234 4 8 9 17 7 8 ... 12 12 17 19 17 13 16 19 14 47
439305 900010000004 47 28 19 0 1 2 0 3 1 ... 0 1 2 2 2 4 2 2 0 0
439306 900010000010 5 1 4 0 1 0 0 0 0 ... 0 0 1 0 0 0 0 0 0 1

5 rows × 36 columns

Leggo le geometrie delle sezioni 2021

In [43]:
sezioni_21_geom = gpd.read_parquet(os.path.join(c.work_path, 'ISTAT', 'DATA', 'sezioni_21_geom.parquet'))
sezioni_21_geom = sezioni_21_geom[['sez21_id', 'geometry']].copy()
sezioni_21_geom['sez_a_m2'] = sezioni_21_geom['geometry'].area
In [44]:
sezioni_21_geom.head()
Out[44]:
sez21_id geometry sez_a_m2
0 10010000001 POLYGON ((404026.314 5024216.152, 403976.239 5... 269456.721061
1 10010000002 POLYGON ((403621.103 5024831.964, 403622.510 5... 466547.650192
2 10010000006 POLYGON ((403424.939 5025678.139, 403430.532 5... 17608.572199
3 10010000007 POLYGON ((403859.447 5026020.664, 403869.348 5... 191053.288816
4 10010000008 POLYGON ((404912.809 5025930.662, 404931.933 5... 45837.207412

Facciamo overlay con gli esagoni!

In [45]:
start = datetime.now()
o = grid.overlay(sezioni_21_geom.to_crs(4326), how='intersection')
o['int_a_m2'] = o['geometry'].to_crs(3857).area
o['k'] = o['int_a_m2']/o['sez_a_m2']
o = o[[f'h3_{resolution:02}', 'sez21_id', 'geometry', 'k']].copy()
int_et = str(datetime.now()-start).split(".")[0]
print(f'Esecuzione: {int_et}')

Esecuzione: 0:01:20
In [46]:
o.head()
Out[46]:
h3_08 sez21_id geometry k
0 883865a699fffff 1110490000014 POLYGON ((8.76850 39.15124, 8.76872 39.15190, ... 0.001786
1 883865a691fffff 1110490000014 POLYGON ((8.77147 39.14303, 8.77251 39.14624, ... 0.025561
2 883865a6d1fffff 1110490000014 POLYGON ((8.75941 39.15792, 8.75486 39.16126, ... 0.082375
3 883865a6d5fffff 1110490000014 POLYGON ((8.74873 39.15974, 8.75486 39.16126, ... 0.274187
4 883865a695fffff 1110490000014 POLYGON ((8.76534 39.14151, 8.77147 39.14303, ... 0.091525

Merge tra indicatori e geometrie

In [47]:
sezioni_21 = sezioni_21_ind.merge(o, on='sez21_id', how='right')
sezioni_21 = gpd.GeoDataFrame(sezioni_21)
sezioni_21.fillna(0, inplace=True)
sezioni_21.shape
Out[47]:
(105962, 39)
In [48]:
sezioni_21[sezioni_21[f'h3_{resolution:02}']==cell_sel].explore()
Out[48]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [49]:
indicatori = [etichetteVariabili[key] for key in etichetteVariabili]
In [50]:
for col in indicatori:
    sezioni_21[col] = sezioni_21[col] * sezioni_21['k']

Ricalcolo gli indicatori aggregati sull'esagono (perdo le geometrie)

In [51]:
grid_ind = sezioni_21.groupby(f'h3_{resolution:02}', as_index=False).agg({indicatore: 'sum' for indicatore in indicatori})

Ci ri-attacco le geometrie

In [52]:
grid_adv = grid.merge(grid_ind, on=f'h3_{resolution:02}')
In [53]:
grid_adv.head(100).explore()
Out[53]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [54]:
grid_adv.shape
Out[54]:
(31602, 38)

Mergiamo tutto: popolazione e OSM¶

Selezioniamo solo i pois con le iso a 15 min

In [55]:
g15 = g[g['contour']==15]

ristrutturiamo la tabella con i dati OSM

In [56]:
pivot = pd.pivot_table(g15, values='size', index=f'h3_{resolution:02}', columns='fclass', aggfunc='sum', fill_value=0).reset_index()

Possiamo mergiare i dati della popolazione con quelli osm

In [57]:
final = grid_adv.merge(pivot, on=f'h3_{resolution:02}', how='left').fillna(0)
In [58]:
final[final[f'h3_{resolution:02}']==cell_sel].explore()
Out[58]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]:
In [ ]:

Esempio di creazione della macrocategoria dell'indicatore

In [59]:
final['sport'] = 0
for col in ['sports_centre', 'pitch', 'stadium']:
    final['sport'] += final[col]
In [ ]: