Lavorare con i dati geospaziali a livello di codice postale

'Utilizzare i dati geospaziali a livello di CAP'

Come associare i codici postali “punto” con i dati del censimento dell’area

In alcuni paesi, i codici postali sono punti o percorsi e non aree. Ad esempio, le ultime tre cifre di un codice postale canadese corrispondono all’unità di consegna locale, che può corrispondere alle case su un lato di una strada o a un percorso rurale. Allo stesso modo, i codici postali del Regno Unito hanno un formato del tipo “YO8 9UR”. Questo potrebbe essere anche un singolo edificio a Londra. In un codice postale statunitense 5+4, gli ultimi quattro numeri determinano un percorso di consegna postale (quindi, un insieme di indirizzi) e non un’area. Contrariamente alla credenza comune, i codici ZIP statunitensi a 5 cifre non sono nemmeno aree: sono semplicemente una raccolta di percorsi postali a 5+4, tipicamente serviti da un’unica posta.

La Francia, come si addice all’ideatore del sistema metrico, è molto logica. In Francia, i codici postali corrispondono a un’area: le ultime due cifre corrispondono all’arrondissement, quindi 75008 corrisponde all’ottavo arrondissement di Parigi ed è veramente un’area. Le rotte di consegna della posta sono probabilmente subottimali, però.

Poiché le persone e i negozi hanno indirizzi, che hanno codici postali associati, la maggior parte dei dati dei consumatori viene riportata a livello di codice postale. Per effettuare calcoli come la copertura territoriale, la quota di mercato, ecc. è necessario determinare l’estensione territoriale di un codice postale. Questo è facile in Francia, ma sarà difficile in qualsiasi paese in cui i codici postali siano percorsi postali e non aree.

I codici postali del Regno Unito sono indirizzi di consegna postale reale, non aree. Questo è vero anche per il Canada e gli Stati Uniti. Foto di Monty Allen su Unsplash

Poiché i loro codici postali sono indirizzi di consegna postale, è possibile disegnare infiniti poligoni per suddividere il Regno Unito/Canada/USA in “regioni” di codici postali validi. Questo è il motivo per cui i dati demografici del Regno Unito sono pubblicati dall’Ufficio di Statistica Nazionale (ONS) in aree amministrative (come le contee) e non nei codici postali. Il censimento degli Stati Uniti pubblica dati a livello di “area di tabulazione del codice postale” (ZCTA), e i dati di voto degli Stati Uniti sono pubblicati a livello di contea. Quando si lavora con dati del Regno Unito/Canada/USA, spesso si ha una miscela simile di indirizzi (che sono punti) e dati spaziali raccolti su un’area. Come si associano insieme?

A titolo illustrativo, metterò insieme i dati del codice postale del Regno Unito e i dati del censimento in questo articolo.

Se siete di fretta, potete scaricare i risultati della mia analisi da https://github.com/lakshmanok/lakblogs/tree/main/uk_postcode — ci sono un paio di file CSV lì e contengono i dati di cui potreste aver bisogno.

ukpopulation.csv.gz ha le seguenti colonne:

codice_postale,latitudine,longitudine,codice_area,nome_area,tutte_le_persone,donne,uomini

ukpostcodes.csv.gz ha una colonna in più: il poligono per ogni codice postale nel formato WKT:

codice_postale,latitudine,longitudine,codice_area,nome_area,tutte_le_persone,donne,uomini,geometria_wkt

Si prega di notare che l’uso dei dati o del codice è a proprio rischio: viene distribuito “COME È”, SENZA GARANZIE O CONDIZIONI DI ALCUN TIPO, esplicite o implicite.

In questo articolo, passerò attraverso come ho creato il dataset in quel repository GitHub. Potete seguirmi usando il notebook uk_postcodes.ipynb .

Dati grezzi

Partiamo da tre fonti di dati grezzi rilasciate con licenza Open Government del Regno Unito:

  1. Free Map Tools ha un file scaricabile che contiene la latitudine e la longitudine del centroide per ogni codice postale. Questo non è sufficiente per l’analisi spaziale, poiché non è possibile fare cose come ST_INTERSECTS con un solo punto. Ma è un buon punto di partenza.
  2. L’ONS ha pubblicato dati di censimento per aree come “County Durham”. Questi non sono codici postali, però. Sono tipicamente a livello di ward, contea o regione.
  3. L’ufficio di statistica del Regno Unito ha gentilmente identificato le aree associate a ogni codice postale . Ogni codice postale appartiene a diverse aree definite per diversi scopi e risoluzioni. Ciò include, ma non si limita a, la parrocchia, il ward, la contea e la regione in cui si trova la regione. A titolo di completezza, ecco tutte le altre associazioni disponibili (a seconda del vostro dataset spaziale, potreste avere bisogno delle altre colonne):
pcd,pcd2,pcds,dointr,doterm,oscty,ced,oslaua,osward,parish,usertype,oseast1m,osnrth1m,osgrdind,oshlthau,nhser,ctry,rgn,streg,pcon,eer,teclec,ttwa,pct,itl,statsward,oa01,casward,park,lsoa01,msoa01,ur01ind,oac01,oa11,lsoa11,msoa11,wz11,sicbl,bua11,buasd11,ru11ind,oac11,lat,long,lep1,lep2,pfa,imd,calncv,icb,oa21,lsoa21,msoa21

Il mio notebook scarica i file di dati utilizzando wget:

mkdir -p indatacd indataif [ ! -f census2021.xlsx ]; then    wget -O census2021.xlsx https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/populationandhouseholdestimatesenglandandwalescensus2021/census2021/census2021firstresultsenglandwales1.xlsxfi

Lettura dei dati

La lettura diretta di un file CSV in Pandas è semplice:

import pandas as pdpd.read_csv(POSTCODES_CSV)

Questo mi fornisce il centroide lat-lon di ogni codice postale:

Centroide lat-lon per ogni codice postale

Esistono molti pacchetti che consentono di leggere file Excel in Pandas, ma ho deciso di usare DuckDB perché lo utilizzerò in seguito nel notebook per unire i tre dataset utilizzando SQL:

import duckdbconn = duckdb.connect()ukpop = conn.execute(f"""install spatial;load spatial;SELECT   * FROMst_read('{POPULATION_XLS}', layer='P01')""").df()

Questo file Excel ha 7 righe di informazioni di intestazione che posso eliminare. Rinomino anche le colonne in variabili significative:

ukpop = ukpop.drop(range(0,7), axis=0).rename(columns={    'Field1': 'area_code',     'Field2': 'area_name',     'Field3': 'all_persons',     'Field4': 'females',     'Field5': 'males'})

Quella era la scheda denominata P01. Nota che la scheda P04 ha informazioni sulla densità di popolazione, ma non è utile perché la popolazione non è distribuita uniformemente sul codice area. Ricaveremo la popolazione di ogni codice postale.

Salvo queste informazioni in un file CSV in modo da poterlo leggere facilmente da DuckDB.

ukpop.to_csv("temp/ukpop.csv", index=False)

Allo stesso modo, estraggo le colonne necessarie dal file dell’ufficio di statistica del Regno Unito e lo scrivo in un file CSV:

onspd = pd.read_csv(ONSPD_CSV)onspd = onspd[['pcd', 'oslaua', 'osward', 'parish']]onspd.to_csv("temp/onspd.csv", index=False)

Associazione dei dati

Ora, possiamo usare DuckDB per unire i tre dataset preparati per ottenere la densità di popolazione di ogni codice postale. Perché DuckDB? Anche se posso fare l’unione in Pandas, trovo SQL molto più leggibile. Inoltre, mi ha dato una scusa per utilizzare la nuova tendenza.

Unisco i dataset leggendoli prima in DuckB usando read_csv_auto. Quindi, cerco l’area (parrocchia, distretto, o contea) in cui si trova il codice postale e trovo l’area (parrocchia, distretto o contea) in cui vengono riportati i dati sulla densità di popolazione:

/* pcd,oslaua,osward,parish */WITH onspd AS (    SELECT       *     FROM    read_csv_auto('temp/onspd.csv', header=True)),/* area_code,area_name,all_persons,females,males */ukpop AS (    SELECT       *     FROM    read_csv_auto('temp/ukpop.csv', header=True)),/* id,postcode,latitude,longitude */postcodes AS (    SELECT       *     FROM    read_csv_auto('indata/ukpostcodes.csv', header=True)),/* postcode, area_code */postcode_to_areacode AS (  SELECT     pcd AS postcode,    ANY_VALUE(area_code) as area_code  FROM onspd  JOIN ukpop   ON (area_code = oslaua OR area_code = osward OR area_code = parish)  GROUP BY pcd)SELECT  postcode, latitude, longitude, /* from postcodes */  area_code, area_name, /* from ukpop */  all_persons,females,males /* from ukpop, but has to be spatially corrected */FROM postcode_to_areacodeJOIN postcodes USING (postcode)JOIN ukpop USING (area_code)

Si noti che le quantità spaziali sono scalari che corrispondono all’intera area e non al codice postale. Devono essere suddivisi tra i codici postali.

Dividere le quantità di area tra i codici postali

Tutti gli individui, le femmine e i maschi corrispondono all’intera area, non a un codice postale specifico. Potremmo farlo proporzionalmente in base all’area del codice postale, ma ci sono infiniti poligoni che possono adattarsi ai codici postali, e come vedremo più avanti, l’estensione areale dei codici postali vicino ai parchi e ai laghi è un po’ dubbia. Quindi faremo qualcosa di semplice che ci dia una risposta singola e unica: divideremo il valore scalare in modo uniforme tra tutti i codici postali nell’area! Questo non è strano come sembra: nei quartieri ad alta densità, ci sono più codici postali, quindi dividere in modo equo tra i codici postali equivale a distribuire la quantità scalare proporzionalmente alla densità di popolazione.

npostcodes = ukpop.groupby('area_code')['postcode'].count()for col in ['females', 'males', 'all_persons']:    ukpop[col] = ukpop.apply(lambda row:  row[col]/npostcodes[row['area_code']], axis=1)

A questo punto, abbiamo la quantità per ogni codice postale – questa è l’associazione di cui avevamo bisogno:

Quindi, scrivilo:

ukpop.to_csv("ukpopulation.csv", index=False)

Estensione areale dei codici postali

Per molte analisi, vogliamo che i codici postali siano aree e non punti. Anche se ci sono infiniti poligoni che possiamo usare per dividere il Regno Unito in modo che ci sia solo un centroide di codice postale in ciascun poligono, esiste una nozione del “migliore” poligono. Questo è la partizione di Voronoi, che divide l’area in modo che ogni punto appartenga al codice postale più vicino ad esso:

Analisi di Voronoi di 20 punti, da Wikipedia

Per calcolarlo, possiamo utilizzare scipy:

import numpy as npfrom scipy.spatial import Voronoi, voronoi_plot_2dpoints = df[['latitude', 'longitude']].to_numpy()vor = Voronoi(points)

Assumo qui che l’area sia abbastanza piccola da non avere molta differenza tra la distanza geodesica e la distanza euclidea calcolata dalla latitudine e longitudine. I codici postali del Regno Unito sono abbastanza piccoli da soddisfare questa condizione.

Il risultato è organizzato in modo tale che, per ogni punto, ci sia una regione composta da un insieme di vertici. Possiamo creare una stringa poligono WKT per ogni punto usando:

def make_polygon(point_no):    region_no = vor.point_region[point_no]    region = vor.regions[region_no]    if len(region) >= 3:        # chiudi l'anello        closed_region = region.copy()        closed_region.append(closed_region[0])        # crea un WKT dei punti        polygon = "POLYGON ((" + ','.join([ f"{vor.vertices[v][1]} {vor.vertices[v][0]}" for v in closed_region]) + "))"        return polygon    else:        return None

Ecco un esempio di risultato:

POLYGON ((-0.32491691953979235 51.7393550489536,-0.32527234008402217 51.73948967705648,-0.32515738641624575 51.73987124225542,-0.3241646650618929 51.74087626616231,-0.3215663358407994 51.742660660928614,-0.32145633473723817 51.742228570262824,-0.32491691953979235 51.7393550489536))

Possiamo creare un GeoDataFrame e tracciare un sottoinsieme di codici postali:

import geopandas as gpdfrom shapely import wktdf['geometry'] = gpd.GeoSeries.from_wkt(df['geometry_wkt'])gdf = gpd.GeoDataFrame(df, geometry='geometry')gdf[gdf['area_name'] == 'St Albans'].plot()
Estensione spaziale del codice postale per St. Albans. Immagine generata dall'autore.

Ecco Birmingham:

gdf[gdf['area_name'] == 'Birmingham'].plot()
Estensione spaziale del codice postale per Birmingham. Immagine generata dall'autore.

Aree senza popolazione

Notare la zona a forma di corno in alto e la grande area blu al centro. Che succede? Diamo un’occhiata a Birmingham su Google Maps:

Birmingham, screenshot di Google Maps di autore.

Notate le aree di parchi? La Royal Mail non è tenuta a consegnare a nessuno lì. Quindi non ci sono codici postali in quelle zone. Di conseguenza, i codici postali vicini vengono “estesi” in quelle aree. Ciò causerà problemi nei calcoli spaziali poiché quei codici postali appariranno molto più grandi di quanto siano effettivamente.

Per risolvere questo problema, adotterò un approccio piuttosto euristico. Divido il Regno Unito in celle di risoluzione 0,01×0,01 (circa 1 km²) e trovo le celle di griglia che non contengono codici postali al loro interno:

GRIDRES = 0.01min_lat, max_lat = np.round(min(df['latitude']), 2) - GRIDRES, max(df['latitude']) + GRIDRESmin_lon, max_lon = np.round(min(df['longitude']), 2) - GRIDRES, max(df['longitude']) + GRIDRESprint(min_lat, max_lat, min_lon, max_lon)npostcodes = np.zeros([ int(1+(max_lat-min_lat)/GRIDRES), int(1+(max_lon-min_lon)/GRIDRES) ])for point in points:    latno = int((point[0] - min_lat)/GRIDRES)    lonno = int((point[1] - min_lon)/GRIDRES)    npostcodes[latno, lonno] += 1unpop = []for latno in range(len(npostcodes)):    for lonno in range(len(npostcodes[latno])):        if npostcodes[latno][lonno] == 0:            # nessuno vive qui.            # creo un codice postale per questa posizione            # codice postale latitudine longitudine area_code area_name persone_per_km²            unpop.append({                'codice_postale': f'NONPOP {latno}x{lonno}',                'latitudine': min_lat + latno * 0.01,                'longitudine': min_lon + lonno * 0.01,                'tutte_le_persone': 0            })               

Creeremo codici postali falsi nel centro di queste celle di griglia prive di popolazione e assegneremo una densità di popolazione zero a queste. Aggiungeremo questi codici postali falsi ai codici postali effettivi e ripeteremo l’analisi di Voronoi:

df2 = pd.concat([df, pd.DataFrame.from_records(unpop)])points = df2[['latitudine', 'longitudine']].to_numpy()vor = Voronoi(points)df2['geometry_wkt'] = [make_polygon(x) for x in range(len(vor.point_region))]df2['geometry'] = gpd.GeoSeries.from_wkt(df2['geometry_wkt'])gdf = gpd.GeoDataFrame(df2, geometry='geometry')

Ora, quando disegniamo Birmingham, otteniamo qualcosa di molto più bello:

Poligoni dei codici postali di Birmingham. Immagine generata dall'autore

È questo dataframe che ho salvato come secondo file csv:

gdf.to_csv("ukpostcodes.csv", index=False)

[Opzionale] Caricamento su BigQuery

Possiamo caricare il file CSV su BigQuery e effettuare alcune analisi spaziali, ma è meglio far sì che BigQuery analizzi prima l’ultima colonna di stringhe come una geometria e avere i dati raggruppati per codice postale:

CREATE OR REPLACE TABLE uk_public_data.postcode_popgeo2CLUSTER BY postcodeASSELECT   * EXCEPT(geometry_wkt),  SAFE.ST_GEOGFROMTEXT(geometry_wkt, make_valid=>TRUE) AS geometry,FROM uk_public_data.postcode_popgeo

Ora possiamo interrogarlo facilmente. Ad esempio, possiamo utilizzare ST_AREA per i codici postali:

SELECT   COUNT(*) AS num_postcodes,  SUM(ST_AREA(geometry))/1e6 AS total_area,  SUM(all_persons) AS population,  area_name FROM uk_public_data.postcode_popgeo2 GROUP BY area_name ORDER BY population DESC

Sommario

L’analisi spaziale richiede spesso una estensione areale, non solo delle posizioni puntuali. Nei paesi in cui i codici postali sono punti/percorsi, ci sono infinite modalità per generare un’estensione spaziale poligonale per i codici postali. Un approccio ragionevole è utilizzare le regioni di Voronoi per creare poligoni che contengano tali codici postali. Tuttavia, se lo si fa, si otterranno poligoni anormalmente grandi vicino a laghi o parchi dove l’ufficio postale non recapita la posta. Per risolvere questo problema, è anche possibile suddividere il paese in griglia e creare codici postali artificiali nelle celle della griglia non popolate. In questo articolo ho dimostrato come farlo per il Regno Unito. Il notebook associato può essere adattato ad altri luoghi.

Passaggi successivi

  1. Controlla il codice completo su https://github.com/lakshmanok/lakblogs/blob/main/uk_postcode/uk_postcodes.ipynb
  2. Scarica i dati da https://github.com/lakshmanok/lakblogs/tree/main/uk_postcode