Binning spaziale con Google BigQuery

Il binning dei dati è un'utile pratica comune nella scienza dei dati e nell'analisi dei dati in diversi modi: la discretizzazione di una variabile continua in Machine Learning o semplicemente la creazione di un istogramma per facilitare la visualizzazione sono un paio di esempi familiari. Spesso i valori sono raggruppati in bin di uguale dimensione, ad esempio intervalli di lunghezza costante o rettangoli di area costante rispettivamente per variabili unidimensionali e bidimensionali, ma finché ogni valore cade in uno e un solo bin siamo completamente liberi di definire loro come ci piace. Un banale esempio di binning potrebbe essere realizzato troncando un numero decimale per conservare solo la sua parte integrale o qualsiasi numero di cifre decimali, a seconda della dimensione dell'intervallo necessario.

Vogliamo qui concentrarci su un caso particolare di variabili continue, ovvero longitudine e latitudine relative a dati georeferenziati. Come nel caso generale, abbiamo modi virtualmente infiniti per definire contenitori bidimensionali non sovrapposti che ricoprono la superficie della Terra: una scelta semplice potrebbe essere quella di utilizzare qualsiasi divisione amministrativa esistente in stati, province, città e distretti per citarne alcuni. L'unica limitazione qui è avere accesso a tali dati geografici, che fortunatamente sono spesso pubblici: BigQuery ci consente di accedere a diversi set di dati geografici pubblici come il set di dati geo us boundaries , che contiene tabelle che hanno poligoni che rappresentano i confini degli Stati Uniti.

Contenitori geometricamente definiti

Questa soluzione non può essere sempre soddisfacente in quanto vorremmo definire i bin anche geometricamente in termini di funzioni matematiche come normalmente facciamo con altre variabili continue. Sarebbe bello suddividere una data area geografica in diversi poligoni più piccoli che sono regolari in un certo senso. Qui arriva la parte difficile: la superficie terrestre (presumibilmente ) non è piatta, quindi la sua rappresentazione planare causerà inevitabilmente alcune distorsioni nelle forme o nelle lunghezze. Si sarebbe tentati di dimenticare per un momento la natura angolare delle coordinate che usiamo abitualmente per localizzare un punto, cioè longitudine e latitudine, e trattarle come coordinate cartesiane. Questo ci permetterebbe di definire facilmente bin geografici nello stesso modo suggerito prima per il caso generale di variabili bidimensionali: possiamo troncare sia la longitudine che la latitudine e mantenere solo un dato numero di cifre decimali. Mostriamo cosa otterremmo, usando Google BigQuery per creare i poligoni e GoogleGeo Viz per visualizzarli (Figura 1):

with starting_point as(
select
  0 sp_lat,
  18 sp_long,
  1 digits
)
select
  st_makepolygon(st_makeline(array[
  st_geogpoint(pol_long,pol_lat), 
  st_geogpoint(pol_long,pol_lat + pow(10,-digits)),
  st_geogpoint(pol_long + pow(10,-digits),pol_lat + pow(10,-digits)),
  st_geogpoint(pol_long + pow(10,-digits),pol_lat) 
  ])) polygon
from(
  select
    digits,
    trunc(sp_long,digits) + x*pow(10,-digits) pol_long,
    trunc(sp_lat,digits) + y*pow(10,-digits) pol_lat,
  from
    starting_point,
    unnest(generate_array(-10,10,1)) x,
    unnest(generate_array(-10,10,1)) y)

Sembra carino!
Cambiamo il punto di partenza:

with starting_point as(
select
  59 sp_lat,
  18 sp_long,
  1 digits
)

Questa griglia ha un aspetto leggermente diverso dall'ultima. Come mai?
Ricorda che longitudine e latitudine sono coordinate angolari, e mentre spostando 1° lungo un meridiano (longitudine costante) vicino ai poli si ottiene lo stesso spostamento misurato vicino all'equatore, lo stesso non vale per i movimenti lungo le parallele (latitudine costante): man mano che ci avviciniamo ai poli una variazione di 1° di longitudine è progressivamente minore ed esattamente 0 ai poli. Possiamo facilmente convincerci di una tale caratteristica misurando la distanza tra due vertici qualsiasi dei poligoni variando il punto di partenza:

st_length(st_makeline(st_geogpoint(pol_long,pol_lat),st_geogpoint(pol_long,pol_lat + pow(10,-digits)))) length_constant_longitude,
st_length(st_makeline(st_geogpoint(pol_long,pol_lat) ,st_geogpoint(pol_long + pow(10,-digits),pol_lat))) length_constant_latitude

Hanno la stessa larghezza ma diverse altezze. Ciò è causato esclusivamente dalle inevitabili distorsioni della proiezione utilizzata per mappare la superficie curva della Terra su un piano. Ecco perché la Groenlandia sembra grande quasi quanto l'Africa!

Realizzazione di contenitori quadrati

La tecnica del binning appena proposta è perfettamente valida in quanto ogni punto geografico cade in uno ed un solo poligono (che d'ora in poi chiameremo piastrella ) di questa griglia (una volta estesa a tutta la superficie), ma manca della soddisfacente regolarità visiva di una griglia composta da quadrati perfettamente allineati (almeno lontano dall'equatore!). A un algoritmo di apprendimento automatico potrebbe non interessare e, in tal caso, ci sono probabilmente soluzioni migliori in cui le lunghezze o le aree sono mantenute costanti in tutta la griglia, ma se siamo più interessati a un output visivo ordinato, questa soluzione si adatta meglio. Seguendo la struttura dell'esempio precedente abbiamo:

with starting_point as( -- step 1
select
  0 sp_lat,
  18 sp_long,
  12 zoom
)
select
 zoom,
 tile_x,
 tile_y,
 st_makepolygon(st_makeline([
  st_geogpoint(ne_long,ne_lat),
  st_geogpoint(se_long,se_lat),
  st_geogpoint(sw_long,sw_lat),
  st_geogpoint(nw_long,nw_lat)])) polygon,
from( --step 3
 select
  zoom,
  tile_x,
  tile_y,
  180/acos(-1)*(2*acos(-1)*tile_x/pow(2,zoom)-acos(-1)) nw_long,
  360/acos(-1)*(atan(exp(acos(-1)-2*acos(-1)/pow(2,zoom)*tile_y))-acos(-1)/4) nw_lat,
  180/acos(-1)*(2*acos(-1)*(tile_x+1)/pow(2,zoom)-acos(-1)) ne_long,
  360/acos(-1)*(atan(exp(acos(-1)-2*acos(-1)/pow(2,zoom)*tile_y))-acos(-1)/4) ne_lat,
  180/acos(-1)*(2*acos(-1)*(tile_x+1)/pow(2,zoom)-acos(-1)) se_long,
  360/acos(-1)*(atan(exp(acos(-1)-2*acos(-1)/pow(2,zoom)*(tile_y+1)))-acos(-1)/4) se_lat,
  180/acos(-1)*(2*acos(-1)*(tile_x)/pow(2,zoom)-acos(-1)) sw_long,
  360/acos(-1)*(atan(exp(acos(-1)-2*acos(-1)/pow(2,zoom)*(tile_y+1)))-acos(-1)/4) sw_lat        
 from( --step 2
  select
   zoom,
   trunc((sp_long*acos(-1)/180+acos(-1))/(2*acos(-1))*pow(2,zoom)) + x tile_x,
   trunc(pow(2,zoom)/(2*acos(-1))*(acos(-1)-safe.ln(tan(acos(-1)/4 + (sp_lat/180*acos(-1))/2)))) + y tile_y,
  from
   starting_point,
   unnest(generate_array(-10,10,1)) x,
   unnest(generate_array(-10,10,1)) y
     )
    )

1. Come nell'esempio precedente, iniziamo sia definendo il punto di partenza attorno al quale generare ciascuna tessera della griglia, sia scegliendo la dimensione desiderata tramite il parametro zoom . Ciò che conta è che il livello di zoom più basso è 0 (cioè c'è solo una tessera che copre l'intera mappa) e ogni volta che saliamo di un livello il numero di tessere quadruplica. Qui scegliamo arbitrariamente un livello di 12.
2. Determiniamo le " coordinate delle tessere " di ogni quadrato nella griglia dai due array di numeri interi. Queste coordinate e il livello di zoom formeranno insieme un identificatore univoco per ogni riquadro. Nota che usiamo acos(-1) per ottenere il valore di π.
3. Successivamente, applichiamo le funzioni inverse appena utilizzate per ottenere le coordinate dei quattro vertici di ciascun quadrato.
4. Infine, generiamo il riquadro come tipo di dati geografico utilizzando le consuete funzioni di geografia. Nota che questi due ultimi passaggi sono necessari solo se vogliamo visualizzare la griglia sulla mappa.

E se cambiassimo il punto di partenza? Diamo un'occhiata (Figura 5):

È chiaro che ora abbiamo due griglie identiche composte da quadrati regolari della stessa dimensione apparente, indipendentemente da dove si trovino sulla mappa. Vogliamo sottolineare ancora una volta che sembrano solo uguali, quando in realtà non lo sono. La valutazione della loro area rivela questa caratteristica: l'area media delle tegole vicino all'equatore è di circa 95,5 km², mentre l'area media delle tegole più vicine al Polo Nord è di circa 25,3 km²!

Siamo finalmente pronti per utilizzare questa tecnica per creare molte fantastiche visualizzazioni sia su Geo Viz che su Data Studio!

Divertiti!