Nonostante i numerosi vantaggi, l’accesso alla natura e agli spazi verdi sta diventando sempre più difficile nelle aree altamente urbanizzate. Alcuni temono che le comunità svantaggiate siano più esposte a questi problemi. Qui propongo un modo basato sui dati per esplorare questo aspetto.
In particolare, pongo una questione di sviluppo urbano che ultimamente ha guadagnato interesse tra gli ambienti professionali e i governi locali – ora come uguaglianza verde. Questo concetto si riferisce alle disparità tra le persone che accedono agli spazi verdi in diverse parti di una particolare città. Qui esploro la sua dimensione finanziaria e vedo se ci sono relazioni chiare tra l’area verde disponibile pro capite e il livello economico di quella stessa unità urbana.
Esplorerò due diverse risoluzioni spaziali della città: i distretti e i distretti censimentari, utilizzando gli Shapefile Esri forniti dal portale Open Data del governo austriaco. Incorporerò anche dati statistici tabellari (popolazione e reddito) nelle aree amministrative georeferenziate. Quindi, ho sovrapposto le aree amministrative con un set di dati ufficiale sulle aree verdi, registrando la posizione di ogni spazio verde in un formato geospaziale. Quindi, combino queste informazioni e quantifico la dimensione pro capite dello spazio verde totale di ciascun distretto urbano. Infine, metto in relazione lo stato finanziario di ciascuna area, catturato dal reddito netto annuo, al rapporto pro capite dell’area verde per vedere se emergono eventuali modelli.
Diamo un’occhiata al portale Open Data del governo austriaco Qui.
Mentre stavo scrivendo questo articolo, la traduzione in inglese del sito web non funzionava davvero, quindi, invece di fare affidamento sui miei 12 anni di lezioni di tedesco da tempo dimenticati, ho usato DeepL per navigare tra le sottopagine e migliaia di set di dati.
Successivamente, ho raccolto un paio di file di dati, sia georeferenziati (shapefile Esri) che semplici dati tabulari, che utilizzerò per l’analisi successiva. I dati che ho raccolto:
Confini: i confini amministrativi delle seguenti unità spaziali a Vienna:
Uso del territorio: informazioni sull’ubicazione degli spazi verdi e delle aree edificate:
Statistiche: dati sulla popolazione e sul reddito corrispondenti al livello socioeconomico di un’areaUN:
- Popolazione per distrettoregistrato annualmente dal 2002 e memorizzato suddiviso in base a gruppi di età di 5 anni, sesso e nazionalità originaria
- Popolazione per distretto di censimentoregistrato annualmente dal 2008 e archiviato suddiviso in base a tre gruppi irregolari di età, sesso e origine
- Reddito netto medio dal 2002 nei distretti di Vienna, espressi in euro per dipendente annuo
Inoltre, ho archiviato i file di dati scaricati in una cartella locale denominata data.
2.1 Confini amministrativi
Per prima cosa, leggi e visualizza i diversi shape file contenenti ciascun livello di confine amministrativo per avere una presa più stretta sulla città a portata di mano:
folder = 'data'
admin_city = gpd.read_file(folder + '/LANDESGRENZEOGD')
admin_district = gpd.read_file(folder + '/BEZIRKSGRENZEOGD')
admin_census = gpd.read_file(folder + '/ZAEHLBEZIRKOGD')display(admin_city.head(1))
display(admin_district.head(1))
display(admin_census.head(1))
Qui prendiamo nota che i nomi delle colonne BEZNR e ZBEZ corrispondono rispettivamente all’ID del distretto e all’ID del distretto del censimento. Inaspettatamente, vengono archiviati/analizzati in diversi formati, numpy.float64 e str:
print(type(admin_district.BEZNR.iloc(0)))
print(type(admin_census.ZBEZ.iloc(0)))pyth
Assicurandoci di avere effettivamente 23 distretti e 250 distretti di censimento come affermato nella documentazione dei file di dati:
print(len(set(admin_district.BEZNR)))
print(len(set(admin_census.ZBEZ)))
Ora visualizza i confini: prima la città, poi i suoi distretti e poi i distretti censuari ancora più piccoli.
f, ax = plt.subplots(1,3,figsize=(15,5))admin_city.plot(ax=ax(0),
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Reds')
admin_district.plot(ax=ax(1),
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Blues')
admin_census.plot(ax=ax(2),
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Purples')
ax(0).set_title('City boundaries')
ax(1).set_title('District boundaries')
ax(2).set_title('Census distrcit boundaries')
Questo codice restituisce le seguenti immagini di Vienna:
2.2 Aree verdi
Ora diamo un’occhiata anche alla distribuzione degli spazi verdi:
gdf_green = gpd.read_file(folder + '/GRUENFREIFLOGD_GRUENGEWOGD')
display(gdf_green.head(3))
Qui, si può notare che non esiste un modo diretto per collegare le aree verdi (ad esempio, nessun ID di distretto aggiunto) ai quartieri – quindi in seguito lo faremo manipolando le geometrie per trovare sovrapposizioni.
Ora visualizza questo:
f, ax = plt.subplots(1,1,figsize=(7,5))gdf_green.plot(ax=ax,
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Greens')
ax.set_title('Green areas in Vienna')
Questo codice mostra dove si trovano le aree verdi a Vienna:
Possiamo notare che i segmenti forestali sono ancora all’interno dei confini amministrativi, il che implica che non tutte le parti della città sono urbanizzate e popolate in modo significativo. Torneremo più avanti su questo aspetto valutando l’area verde pro capite.
2.3 Dati statistici: popolazione, reddito
Infine, diamo uno sguardo ai file di dati statistici. La prima grande differenza è che queste non sono georeferenziate ma semplici tabelle csv:
df_pop_distr = pd.read_csv('vie-bez-pop-sex-age5-stk-ori-geo4-2002f.csv',
sep = ';',
encoding='unicode_escape',
skiprows = 1)df_pop_cens = pd.read_csv('vie-zbz-pop-sex-agr3-stk-ori-geo2-2008f.csv',
sep = ';',
encoding='unicode_escape',
skiprows = 1)
df_inc_distr = pd.read_csv('vie-bez-biz-ecn-inc-sex-2002f.csv',
sep = ';',
encoding='unicode_escape',
skiprows = 1)
display(df_pop_distr.head(1))
display(df_pop_cens.head(1))
display(df_inc_distr.head(1))
3.1. Preparazione dei file di dati statistici
La sottosezione precedente mostra che le tabelle dei dati statistici utilizzano convenzioni di denominazione diverse: hanno identificatori DISTRICT_CODE e SUB_DISTRICT_CODE invece di cose come BEZNR e ZBEZ. Tuttavia, dopo aver letto la documentazione di ciascun set di dati, diventa chiaro che è facile trasformarlo dall’uno all’altro, per cui presento due brevi funzioni nella cella successiva. Elaborerò contemporaneamente i dati a livello di distretti e distretti censimentari.
Inoltre, mi interesseranno solo i (più recenti) valori aggregati e i punti dati delle informazioni statistiche, come la popolazione totale nell’istantanea più recente. Quindi, ripuliamo questi file di dati e conserviamo le colonne che utilizzerò in seguito.
# these functions convert the district and census district ids to be compatbile with the ones found in the shapefiles
def transform_district_id(x):
return int(str(x)(1:3))def transform_census_district_id(x):
return int(str(x)(1:5))
# select the latest year of the data set
df_pop_distr_2 = df_pop_distr(df_pop_distr.REF_YEAR \
==max(df_pop_distr.REF_YEAR))
df_pop_cens_2 = df_pop_cens(df_pop_cens.REF_YEAR \
==max(df_pop_cens.REF_YEAR))
df_inc_distr_2 = df_inc_distr(df_inc_distr.REF_YEAR \
==max(df_inc_distr.REF_YEAR))
# convert district ids
df_pop_distr_2('district_id') = \
df_pop_distr_2.DISTRICT_CODE.apply(transform_district_id)
df_pop_cens_2('census_district_id') = \
df_pop_cens_2.SUB_DISTRICT_CODE.apply(transform_census_district_id)
df_inc_distr_2('district_id') = \
df_inc_distr_2.DISTRICT_CODE.apply(transform_district_id)
# aggregate population values
df_pop_distr_2 = df_pop_distr_2.groupby(by = 'district_id').sum()
df_pop_distr_2('district_population') = df_pop_distr_2.AUT + \
df_pop_distr_2.EEA + df_pop_distr_2.REU + df_pop_distr_2.TCN
df_pop_distr_2 = df_pop_distr_2(('district_population'))
df_pop_cens_2 = df_pop_cens_2.groupby(by = 'census_district_id').sum()
df_pop_cens_2('census_district_population') = df_pop_cens_2.AUT \
+ df_pop_cens_2.FOR
df_pop_cens_2 = df_pop_cens_2(('census_district_population'))
df_inc_distr_2('district_average_income') = \
1000*df_inc_distr_2(('INC_TOT_VALUE'))
df_inc_distr_2 = \
df_inc_distr_2.set_index('district_id')(('district_average_income'))
# display the finalized tables
display(df_pop_distr_2.head(3))
display(df_pop_cens_2.head(3))
display(df_inc_distr_2.head(3))
# and unifying the naming conventions
admin_district('district_id') = admin_district.BEZNR.astype(int)
admin_census('census_district_id') = admin_census.ZBEZ.astype(int)
print(len(set(admin_census.ZBEZ)))
Ricontrolla i valori della popolazione totale calcolati ai due livelli di aggregazione:
print(sum(df_pop_distr_2.district_population))
print(sum(df_pop_cens_2.census_district_population))
Questi due dovrebbero fornire entrambi lo stesso risultato: 1931593 persone.
3.1. Preparazione dei file di dati geospaziali
Ora che abbiamo terminato la preparazione dei dati essenziali dei file statistici, è il momento di abbinare i poligoni dell’area verde ai poligoni dell’area amministrativa. Quindi, calcoliamo la copertura totale dell’area verde di ciascuna area amministrativa. Inoltre, per curiosità, aggiungerò la relativa copertura dell’area verde di ciascuna area amministrativa.
Per ottenere aree espresse in unità SI occorre passare ad un cosiddetto CRS locale, che nel caso di Vienna è EPSG:31282. Per saperne di più leggi di più su questo argomento, proiezione cartografica e sistemi di riferimento di coordinate Qui E Qui.
# converting all GeoDataFrames into the loca crs
admin_district_2 = \
admin_district(('district_id', 'geometry')).to_crs(31282)admin_census_2 = \
admin_census(('census_district_id', 'geometry')).to_crs(31282)
gdf_green_2 = gdf_green.to_crs(31282)
Calcola l’area dell’unità amministrativa misurata in unità SI:
admin_district_2('admin_area') = \
admin_district_2.geometry.apply(lambda g: g.area)admin_census_2('admin_area') = \
admin_census_2.geometry.apply(lambda g: g.area)
display(admin_district_2.head(1))
display(admin_census_2.head(1))
4.1 Calcolare la copertura dell’area verde in ciascuna unità amministrativa
Utilizzerò la funzione di sovrapposizione di GeoPanda per sovrapporre questi due GeoDataFrame di confine amministrativo con GeoDataFrame contenente i poligoni dell’area verde. Quindi, calcolo l’area di ciascuna sezione di area verde che rientra in diverse regioni amministrative. Successivamente, riassumo queste aree al livello di ciascuna area amministrativa, sia distrettuali che distretti censimentari. Nella fase finale, in ciascuna unità di risoluzione, aggiungo le aree amministrative dell’unità ufficiale precedentemente calcolate e calcolo il rapporto tra area totale e area verde per ciascun distretto e distretto censuario.
gdf_green_mapped_distr = gpd.overlay(gdf_green_2, admin_district_2)gdf_green_mapped_distr('green_area') = \
gdf_green_mapped_distr.geometry.apply(lambda g: g.area)
gdf_green_mapped_distr = \
gdf_green_mapped_distr.groupby(by = 'district_id').sum()(('green_area'))
gdf_green_mapped_distr = \
gpd.GeoDataFrame(admin_district_2.merge(gdf_green_mapped_distr, left_on = 'district_id', right_index = True))
gdf_green_mapped_distr('green_ratio') = \
gdf_green_mapped_distr.green_area / gdf_green_mapped_distr.admin_area
gdf_green_mapped_distr.head(3)
gdf_green_mapped_cens = gpd.overlay(gdf_green_2, admin_census_2)
gdf_green_mapped_cens('green_area') = \
gdf_green_mapped_cens.geometry.apply(lambda g: g.area)gdf_green_mapped_cens = \
gdf_green_mapped_cens.groupby(by = 'census_district_id').sum()(('green_area'))
gdf_green_mapped_cens = \
gpd.GeoDataFrame(admin_census_2.merge(gdf_green_mapped_cens, left_on = 'census_district_id', right_index = True))
gdf_green_mapped_cens('green_ratio') = gdf_green_mapped_cens.green_area / gdf_green_mapped_cens.admin_area
gdf_green_mapped_cens.head(3)
Infine, visualizza il rapporto verde per distretto e distretto di censimento! I risultati sembrano avere molto senso, con un livello di verde elevato nelle parti esterne e molto più basso nelle aree centrali. Inoltre, i 250 distretti censiti mostrano chiaramente un quadro più dettagliato e dettagliato delle diverse caratteristiche dei quartieri, offrendo spunti più profondi e più localizzati per gli urbanisti. D’altro canto, le informazioni a livello distrettuale, con unità spaziali dieci volte inferiori, mostrano invece grandi medie.
f, ax = plt.subplots(1,2,figsize=(17,5))gdf_green_mapped_distr.plot(ax = ax(0),
column = 'green_ratio',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
legend = True,
cmap = 'Greens')
gdf_green_mapped_cens.plot(ax = ax(1),
column = 'green_ratio',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
legend = True,
cmap = 'Greens')
Questo blocco di codice restituisce le seguenti mappe:
4.2 Aggiungere informazioni sulla popolazione e sul reddito per ciascuna unità amministrativa
Nella fase finale di questa sezione, mapperemo i dati statistici in aree amministrative. Promemoria: disponiamo di dati sulla popolazione sia a livello di distretti che a livello di distretti censiti. Tuttavia, ho potuto trovare il reddito (indicatore del livello socioeconomico) solo a livello dei distretti. Questo è un normale compromesso nella scienza dei dati geospaziali. Mentre una dimensione (il verde) è molto più approfondita con la risoluzione più alta (distretti censiti), i vincoli sui dati potrebbero obbligarci a utilizzare comunque la risoluzione più bassa.
display(admin_census_2.head(2))
display(df_pop_cens_2.head(2))
gdf_pop_mapped_distr = admin_district_2.merge(df_pop_distr_2, \
left_on = 'district_id', right_index = True)gdf_pop_mapped_cens = admin_census_2.merge(df_pop_cens_2, \
left_on = 'census_district_id', right_index = True)
gdf_inc_mapped_distr = admin_district_2.merge(df_inc_distr_2, \
left_on = 'district_id', right_index = True)
f, ax = plt.subplots(1,3,figsize=(15,5))
gdf_pop_mapped_distr.plot(column = 'district_population', ax=ax(0), \
edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Blues')
gdf_pop_mapped_cens.plot(column = 'census_district_population', ax=ax(1), \
edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Blues')
gdf_inc_mapped_distr.plot(column = 'district_average_income', ax=ax(2), \
edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Purples')
ax(0).set_title('district_population')
ax(1).set_title('census_district_population')
ax(2).set_title('district_average_incomee')
Questo blocco di codici risulta nella seguente figura:
4.3. Calcolo dell’area verde pro capite
Riassumiamo quello che abbiamo ora, tutto integrato in shapefile decenti corrispondenti ai distretti e ai distretti censimentari di Vienna:
A livello di distretti, abbiamo dati sul rapporto tra aree verdi, popolazione e reddito
A livello dei distretti censiti, abbiamo un rapporto tra aree verdi e dati sulla popolazione
Per catturare in modo semplice l’uguaglianza verde, unisco le informazioni sulla dimensione assoluta dell’area verde e sulla popolazione nei distretti e nei distretti censiti e calcolo la quantità totale di area verde pro capite.
Diamo un’occhiata al nostro input: copertura verde e popolazione:
# a plot for the disticts
f, ax = plt.subplots(1,2,figsize=(10,5))gdf_green_mapped_distr.plot(
ax = ax(0),
column = 'green_ratio',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Greens')
gdf_pop_mapped_distr.plot(
ax = ax(1),
column = 'district_population',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Reds')
ax(0).set_title('green_ratio')
ax(1).set_title('district_population')
# a plot for the census disticts
f, ax = plt.subplots(1,2,figsize=(10,5))
gdf_green_mapped_cens.plot(
ax = ax(0),
column = 'green_ratio',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Greens')
gdf_pop_mapped_cens.plot(
ax = ax(1),
column = 'census_district_population',
edgecolor = 'k',
linewidth = 0.5,
alpha = 0.9,
cmap = 'Reds')
ax(0).set_title('green_ratio')
ax(1).set_title('district_population')
Questo blocco di codici risulta nella seguente figura:
Per calcolare l’area verde pro capite, nei passaggi successivi unirò innanzitutto i frame di dati relativi al verde e alla popolazione. Lo farò utilizzando l’esempio dei distretti di censimento perché la loro maggiore risoluzione spaziale ci consente di osservare modelli migliori (se presenti) emergenti. Assicurati di non dividere per zero e segui anche il buon senso; lasciamo perdere quelle aree che non sono popolate.
gdf_green_pop_cens = \
gdf_green_mapped_cens.merge(gdf_pop_mapped_cens.drop( \
columns = ('geometry', 'admin_area')), left_on = 'census_district_id',\
right_on = 'census_district_id')(('census_district_id', \
'green_area', 'census_district_population', 'geometry'))gdf_green_pop_cens('green_area_per_capita') = \
gdf_green_pop_cens('green_area') / \
gdf_green_pop_cens('census_district_population')
gdf_green_pop_cens = \
gdf_green_pop_cens(gdf_green_pop_cens('census_district_population')>0)
f, ax = plt.subplots(1,1,figsize=(10,7))
gdf_green_pop_cens.plot(
column = 'green_area_per_capita',
ax=ax,
cmap = 'RdYlGn',
edgecolor = 'k',
linewidth = 0.5)
admin_district.to_crs(31282).plot(\
ax=ax, color = 'none', edgecolor = 'k', linewidth = 2.5)
Questo blocco di codici risulta nella seguente figura:
Modifichiamo leggermente la visualizzazione:
f, ax = plt.subplots(1,1,figsize=(11,7))ax.set_title('Per-capita green area in\nthe census districts of Vienna',
fontsize = 18, pad = 30)
gdf_green_pop_cens.plot(
column = 'green_area_per_capita',
ax=ax,
cmap = 'RdYlGn',
edgecolor = 'k',
linewidth = 0.5,
legend=True,
norm=matplotlib.colors.LogNorm(\
vmin=gdf_green_pop_cens.green_area_per_capita.min(), \
vmax=gdf_green_pop_cens.green_area_per_capita.max()), )
admin_district.to_crs(31282).plot(
ax=ax, color = 'none', edgecolor = 'k', linewidth = 2.5)
Questo blocco di codici risulta nella seguente figura:
E lo stesso per i distretti:
# compute the per-capita green area scores
gdf_green_pop_distr = \
gdf_green_mapped_distr.merge(gdf_pop_mapped_distr.drop(columns = \
('geometry', 'admin_area')), left_on = 'district_id', right_on = \
'district_id')(('district_id', 'green_area', 'district_population', \
'geometry'))gdf_green_popdistr = \
gdf_green_pop_distr(gdf_green_pop_distr.district_population>0)
gdf_green_pop_distr('green_area_per_capita') = \
gdf_green_pop_distr('green_area') / \
gdf_green_pop_distr('district_population')
# visualize the district-level map
f, ax = plt.subplots(1,1,figsize=(10,8))
ax.set_title('Per-capita green area in the districts of Vienna', \
fontsize = 18, pad = 26)
gdf_green_pop_distr.plot(column = 'green_area_per_capita', ax=ax, \
cmap = 'RdYlGn', edgecolor = 'k', linewidth = 0.5, legend=True, \
norm=matplotlib.colors.LogNorm(vmin=\
gdf_green_pop_cens.green_area_per_capita.min(), \
vmax=gdf_green_pop_cens.green_area_per_capita.max()), )
admin_city.to_crs(31282).plot(ax=ax, \
color = 'none', edgecolor = 'k', linewidth = 2.5)
Questo blocco di codici risulta nella seguente figura:
Mentre le tendenze significative sono chiare: perimetro esterno, più spazio verde per tutti, centro integrato, invertito. Tuttavia, questi due grafici, soprattutto quello più dettagliato a livello di distretti censiti, mostrano chiaramente una variazione nella quantità di spazio verde di cui le persone godono nelle diverse aree. Ulteriori ricerche e l’integrazione di ulteriori fonti di dati, ad esempio, sull’uso del territorio, potrebbero aiutare a spiegare meglio perché quelle aree sono più ricche di verde o di popolazione. Per ora godiamoci questa mappa e speriamo che tutti trovino la giusta quantità di verde nella propria casa!
# merging the greenery, population and financial data
gdf_district_green_pip_inc = \
gdf_green_pop_distr.merge(gdf_inc_mapped_distr.drop(columns = \
('geometry')))
Visualizza la relazione tra la dimensione finanziaria e quella del verde:
f, ax = plt.subplots(1,1,figsize=(6,4))ax.plot(gdf_district_green_pip_inc.district_average_income, \
gdf_district_green_pip_inc.green_area_per_capita, 'o')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('district_average_income')
ax.set_ylabel('green_area_per_capita')
Il risultato di questo blocco di codice è il seguente grafico a dispersione:
A prima vista, il grafico a dispersione non fornisce una forte argomentazione a favore del fatto che i fattori finanziari determinino l’accesso delle persone agli spazi verdi. Onestamente, sono un po’ sorpreso da questi risultati; tuttavia, alla luce degli sforzi consapevoli e di lunga data di Vienna per rendere più verde la propria città, potrebbe essere il motivo per cui non vediamo alcuna tendenza importante qui. Per conferma ho controllato anche le correlazioni tra queste due variabili:
print(spearmanr(gdf_district_green_pip_inc.district_average_income, gdf_district_green_pip_inc.green_area_per_capita))print(pearsonr(gdf_district_green_pip_inc.district_average_income, gdf_district_green_pip_inc.green_area_per_capita))
A causa della distribuzione pesante dei dati finanziari, qui prenderei più seriamente la correlazione di Spearman (0,13), ma anche la correlazione di Pearson (0,30) implica una tendenza relativamente debole, in linea con le mie osservazioni precedenti.
Fonte: towardsdatascience.com
