Il progetto di ricerca è un contributo metodologico all’attività trasversale di elaborazione statistica in PoliS-Lombardia, istituto regionale per il supporto alle politiche di Regione Lombardia con sede a Milano.

Si presenta una forma innovativa di reporting e programmazione statistica con software R, proponendo come esempio i dati regionali del trasporto intermodale. Lo stesso file .docx del progetto di ricerca è fabbricato in questo ambiente di lavoro.

La programmazione in un linguaggio specialistico come R è propedeutica al campo dei Big Data e dell’Intelligenza Artificiale, di interesse comune al ricercatore e all’amministrazione pubblica.

Ambiente, anzi ecosistema R

Software

R è il più popolare linguaggio di programmazione per la statistica, e rimane tra i primi linguaggi generici, con Python e C, a 20 anni dalla versione inaugurale, uscita il 29 febbraio 2000. L’ultima major release è la 4.0 del 2020.

Diversamente dagli omologhi SPSS, SAS e Stata, R aderisce all’importante cyber-filosofia dell’open source. Ross Ihaka e Robert Gentleman derivarono R dal linguaggio proprietario S, prefigurando un codice sorgente gratuito, come l’installazione del software base, delle librerie di estensione, o pacchetti, e dei c.d. «IDE».

IDE (≠ GUI)

Integrated Development Environment. Più che una semplice GUI o Graphic User Interface, l’IDE è una piattaforma di incubazione del processo che parte dall’import dei dati originali e termina col confezionamento ed export dei risultati informativi.

Coniugato con l’IDE, R non è un semplice linguaggio, ma un vero e proprio ambiente di progetto statistico, perché autoconsistente nel richiamo di librerie ad hoc (library()), di routine eseguibili da altri script (source()) e di c.d. «vignette» (vignette()), cioè risorse conoscitive e dimostrazioni di cosa l’utilizzo combinato delle funzioni importate da una libreria può fare.

Il plotting di anteprime tabellari e grafici come istogrammi, curve e nuvole di punti è parte integrante di questo continuo rapporto dialogico ed ellittico con R.

RStudio

RStudio è l’IDE principe di R, il più funzionale e aperto all’interazione agile con altri linguaggi, tra cui il già citato Python e quello SQL dei database relazionali.

L’IDE non sostituisce il software base: per usare RStudio occorre aver installato R. Esiste anche RStudio Cloud, una versione browser che gira su server Linux e non richiede di installare nulla.

Meta-programmazione

In particolare, col sourcing (source()) supportato dall’IDE si possono cogliere i benefici di una buona meta-programmazione del lavoro, specie se di gruppo.

Negli script troppo lunghi è facile confondere l’ordine delle operazioni. Anziché eccedere nel numero di comandi per script, un progetto statistico complesso può essere gestito simultaneamente attraverso script più brevi concepiti secondo l’ordine procedurale di scomposizione del problema.

Questa strategia, comune alla programmazione in altri linguaggi, permette di localizzare più speditamente gli errori.

Tra gli errori, che spesso dipendono da software non aggiornato, si riscontrano quelli sistematici e accidentali. I primi, oltre ai refusi, solitamente riguardano argomenti obbligatori che non vengono adeguatamente specificati all’interno di una funzione. I secondi possono consistere in comandi di cui viene misinterpretato il senso per disattenzione dell’analista.

Community

La comunità online di entusiasti intrattiene rapporti costanti con gli sviluppatori di R base e dei pacchetti di estensione, di fatto tutelando la sfera di influenza del linguaggio nel lavoro degli utenti.

La conseguenza è che questi possono elaborare un dato dall’inizio alla fine nello stesso software senza falsarne la piena riproducibilità con passaggi estranei di tipo «point & click» in Microsoft Excel o altra piattaforma.

Tidyverse

Con la riforma linguistica Tidyverse (Wickham et al., 2019), R da ambiente diventa progressivamente un ecosistema dove la combinazione di funzioni statistiche da parte dell’utente è allo stesso tempo modulare e uniforme. Il senso di Tidyverse sta nella complementarità dei pacchetti che lo compongono.

Non è una riforma esclusiva: è possibile servirsi di R base e Tidyverse nello stesso script.

Il contesto Tidyverse propone una curva di apprendimento molto meno ripida di R base, che si lavori con dati tabellari o spaziali.

Per una padronanza estensiva dei formati, specie quelli meno comuni come le matrici a n-dimensioni (array), è raccomandabile partire con gli operatori e la razionalità di R base, attraverso il quale sono opportunamente compresi i segreti insiti nella struttura di immagazzinamento alla radice dei dati.

Sapere, ad esempio, che una tabella (data.frame/tibble) è un caso particolare di lista (list) i cui vettori colonna hanno eguale lunghezza torna utile nel momento in cui si programma con lo stesso Tidyverse un ciclo di trasformazione su più variabili.

Grammatica

Formalmente, Tidyverse è un ensemble di pacchetti che condividono il medesimo quadro logico e sintattico. La sostanza è che questo rende più intellegibile la relazione tra i comandi di più pacchetti nella loro sequenzialità. Si tratta di un beneficio tanto più apprezzabile quanto più sono numerose le operazioni concorrenti.

Nell’esempio in basso, l’obiettivo è calcolare il coefficiente di variazione (deviazione standard su media) di un vettore numerico con sette valori, di cui due «Not Available» o NA (oggetto vettore). A scopo dimostrativo, si crea una funzione personalizzata (function()) per restituire il coefficiente di variazione (oggetto coefficiente).

Si riconosce il senso di usare Tidyverse in contrapposizione al nesting (annidamento), tipico del software base, preferendo piuttosto il piping (conduttura) (operatore %>%), ovvero l’ideale trasmissione dell’output ottenuto dalla funzione precedente all’input della funzione successiva.

In questo modo, si disbroglia il filo delle operazioni in cui altrimenti tenderebbe a perdersi il soggetto/oggetto del trattamento dati (vettore), e si abilita anche il c.d. «data masking», cioè ogni funzione può fare a meno di ripetere il dato in entrata come primo argomento.

# installo pacchetto tidyverse
install.packages(“tidyverse”)

# carico pacchetto tidyverse
library(tidyverse)

# creo vettore numerico con NA
vettore <- c(5.6, 2.1, NA, 8.3, 6.4, NA, 9.2)

# creo funzione coefficiente variazione
coefficiente <- function(.) sd(.) / mean(.)

# calcolo coefficiente variazione

## nesting
coefficiente(keep(vettore, ~ !is.na(.)))

## piping linea unica
vettore %>% keep(~ !is.na(.)) %>% coefficiente()

## piping indentato
vettore %>% 
  keep(~ !is.na(.)) %>% 
  coefficiente()

Nel piping risiedono vantaggi molteplici: si riducono al minimo gli errori di compitazione (spelling) degli oggetti, e quelli di posizione degli argomenti e delle parentesi; il soggetto/oggetto del trattamento dati è immediatamente riconoscibile.

A mano a mano che aumentano le operazioni da svolgere, si evita di creare oggetti temporanei o di applicare trasformazioni indesiderate a oggetti preesistenti, per il cui ripristino occorrerebbe rieseguire, dopo averle individuate, tutte le routine che li coinvolgono nell’ordine corretto.

A monte di questi vantaggi c’è una comprensione evoluta tra operatore e macchina, più vicina al modo di pensiero sequenziale dell’essere umano. La centralità è del verbo, cioè della funzione da eseguire incrementalmente sui dati originali.

Non a caso un effetto collaterale e voluto di Tidyverse è la sostituzione degli enunciati o «statement» di R base, come for, espressivo di un ciclo o loop, e i condizionali if, else if ed else, con equivalenti funzioni, come map() e possibly(). I cicli programmati con for sono tasselli di codice a sé stanti, mentre map() li rimpiazza senza interrompere la cascata del piping.

Un’anticipazione di questo slittamento linguistico esiste nella famiglia delle c.d. «apply functions», funzioni di automatizzazione con le quali è possibile fare a meno di espliciti e ineleganti enunciati for.

L’omogeneizzazione di R sul modello f(x) è evidentemente una duplice volontà estetica e strumentale: da un linguaggio così normalizzato non ci si aspettano sorprese, ovvero le eccezioni proprie di una qualsiasi lingua naturale. Ogni operazione è eseguita da un comando invocato allo stesso modo, le cui specifiche, obbligatorie e non, corrispondono ai suoi argomenti.

Inoltre, con gli ultimi aggiornamenti, il functional programming e lo stesso piping di Tidyverse sono scritti in C, un linguaggio dalla computazione rapida, quindi riescono a manipolare una grande mole di dati più efficientemente di quanto faccia la scrittura in R tradizionale.

Organizzazione

All’interno dello script, Tidyverse si installa e carica allo stesso modo di un pacchetto ordinario, rispettivamente con install.packages() e library(), ma da questa operazione hanno origine tre costellazioni di pacchetti tra loro simbiotici nel flusso di lavoro:

Tidymodels

Lungamente Tidyverse ha lasciato scoperto il più raffinato dominio della modellazione statistica, sale del data science, di fatto costringendo a ibridare funzioni innovative con la morfologia di R base e convenzioni di nomenclatura dissimili (fai_questo() vs. faiQuesto()).

Con l’ascesa recente di Tidymodels, ensemble sulla falsariga «Core»/«Non-Core»/«Tidy-Compatible», l’ecosistema di R Tidyverse può vantare robustezza e un’esatta traducibilità di quanto era possibile svolgere nella sola grammatica tradizionale.

Passi in avanti notevoli riguardano lo sviluppo di sofisticati algoritmi di machine learning sia per la regressione (variabili continue) che per la classificazione (variabili discrete).

Markdown

La confezione del progetto statistico in R riflette una bellezza ontologica autentica, perché la reportistica prodotta è anch’essa immagine di uno script riproducibile.

Il Markdown, scrittura in «rich text format», non è un semplice scritto, ma è connaturato al codice nella sua essenza. Questo mette in condizione di elaborare documenti dinamici e integrabili anche con dashboard interattive (shiny).

La redazione in Markdown è procedurale e multipiattaforma come un qualunque script R: può contemplare l’inserimento dei c.d. «code chunks» (tranche di codice) comprensivi dei pacchetti di estensione e degli altri linguaggi supportati.

L’ordine di esecuzione dei code chunks è vincolante: ad esempio, un calcolo programmato in linea di testo, se presuppone oggetti creati in un chunk precedente, muta secondo le modifiche a quel chunk.

È possibile configurare anche le opzioni per gestire o affrancarsi da questo vincolo: impostare il blocco dell’esecuzione (argomento eval = FALSE), eseguire il codice nascostamente (include = FALSE), oppure mostrare il solo risultato (echo = FALSE) o il solo contenuto del codice eseguito (results = "hide").

Lo stile del lavoro di ricerca diventa così non solo rieseguibile ma tracciabile, verificabile nei diversi passaggi dal momento di estrazione del dato a quello di divulgazione.

Trasporto intermodale in Lombardia

Contesto

La nota statistica di esempio riporta le fonti di dati e le operazioni di controllo e trattamento curate in PoliS-Lombardia per restituire una base di conoscenza del trasporto merci intermodale nella sua dimensione geografica.

Dote Merci Ferroviaria

La Dote Merci Ferroviaria è un incentivo di Regione Lombardia al trasporto merci intermodale su ferro per treni completi. L’erogazione della misura di incentivo è programmata per tre annualità. La nota è riferita al secondo anno, dal 31/08/2018 al 30/08/2019.
Il contributo consiste nel monitoraggio dei trasporti effettuati dai beneficiari nell’anno di riferimento 2018/19.

Mappatura del trasporto intermodale

Si programma in R (R Core Team, 2020) la lettura e la predisposizione dei dati originali, disaggregati per singolo treno, allo scopo di rappresentare in cartografia i pattern di trasporto intermodale emergenti.

Trattamento dati

Le procedure di controllo preliminare e le elaborazioni tabellari e cartografiche riportate di seguito sono consultabili nello script "read.R".
Attraverso l’import di "read.R", viene tenuta traccia percorribile delle trasformazioni che dal dato originale hanno condotto al dato elaborato e rappresentato.

knitr::read_chunk("read.R")

Campo di osservazione

In continuità con il primo studio di rendicontazione 2017/18, si propone di rilevare i pattern di trasporto che nel periodo di finanziamento si sono verificati anche parzialmente in territorio regionale.
Le analisi comprendono, al netto delle ambiguità di input, tutti i trasporti effettuati dagli operatori beneficiari (i FER) nell’anno 2018/19.
I trasporti fuori campo di osservazione hanno le seguenti caratteristiche:

Rientrano dunque nel campo di osservazione i record che identificano le relazioni ferroviarie avvenute in Lombardia per km maggiori di 0, indipendentemente dalla localizzazione intra/extra regionale/nazionale di origini, destinazioni e frontiere.
Questo è un aspetto importante allo scopo di ottenere una visione nazionale e sovranazionale dei pattern di trasporto per i quali la Lombardia è stata bacino di attrazione.

Fonti di informazione analizzate

La sintesi dei pattern ricercata in questo lavoro assume due fonti di dati: la rendicontazione dei trasporti in Lombardia nell’anno 2018/19, e le liste di località intermedie tra coppie di origini-destinazioni uniche. Il primo è il dato di frequenza. Il secondo è il dato di sequenza del trasporto, indispensabile per adeguare la spazializzazione delle quantità alla realtà geografica dell’infrastruttura.

Entrambe le fonti sono riconosciute e selezionate per essere in formato leggibile dal software di programmazione statistica. A differenza della prima annualità, il database aggiornato al 2018/19 è sprovvisto del dato di sequenza.
Le fonti sono, rispettivamente:

Altre variabili osservate nel file .xlsx sono i codici id degli impianti ferroviari, i codici in formato libero delle relazioni di trasporto e i nomi delle imprese ferroviarie o IF.
Non essendo le liste RFI native della seconda annualità, diviene necessaria un’integrazione manuale con le coppie di origini-destinazioni che scaturiscono dai nuovi pattern di trasporto.

Procedure di controllo

Viene caricato il file di rendicontazione "Copia di RELAZIONI_REGIONE LOMBARDIA.xlsx", rinominato nella forma breve "relazioniRL.xlsx", proveniente dalla directory "wetransfer_91d39c.zip".
Con verifica dei contenuti di ciascuna variabile, sono esplicitati i «data type» delle colonne e vengono impostati nomi più semplici.

library(tidyverse)
library(readxl)
library(fs)

theme_set(theme_void())

file_move(path = "Copia di RELAZIONI_REGIONE LOMBARDIA.xlsx", 
          new_path = "relazioniRL.xlsx")

relazioni <- read_xlsx("relazioniRL.xlsx", 
                       col_types = c("text", "text", "text", "text", "text", "date", 
                                     "text", "text", "text", "text", "text", "text", "text", 
                                     "numeric", "numeric", "text"))

names(relazioni) <- c("FER", "IF", "cod_contratto", "cod_relazione", "cod_trasporto", "data", 
                      "cod_O", "O", "cod_D", "D", "cod_Ft", "Ft", "IF_O", 
                      "km_It", "km_RL", "NOTE")

La struttura del dataset è in formato tabellare, di classe data.frame/tibble, e le sue dimensioni sono 39074 trasporti e 16 variabili.

glimpse(relazioni)

L’estrazione dei primi 5 record/trasporti si presenta in questo modo. Per il display vengono selezionate e rinominate le variabili chiave: operatore, data, stringa impianto e km Italia/Lombardia.

relazioni %>% 
  select(Operatore = FER, 
         Data = data, 
         Origine = O, 
         Destinazione = D, 
         Frontiera = Ft, 
         `km Italia` = km_It, 
         `km Lombardia` = km_RL
         ) %>% 
  head(5)

I data type sono interpretati coerentemente con il contenuto. Il numero degli operatori beneficiari della dote, riscontrabili nel file "Copia di RENDICONTAZ_REGIONE LOMBARDIA.xlsx", è 34.
Nelle relazioni sono inclusi due operatori non beneficiari, il FER "31" e il FER "71", con i seguenti trasporti:

relazioni %>% 
  group_by(Operatore = FER) %>% 
  summarize(`n. trasporti` = n()) %>% 
  filter(Operatore == "31" | Operatore == "71")

Prima delle operazioni di controllo, non si dispone l’eliminazione tout court di questi due operatori perché ancora appartenenti al campo di osservazione.

Sono rilevati ed esclusi attraverso i controlli 2 trasporti con coppie di origini-destinazioni indeterminate nonostante l’indicazione dei punti di frontiera. Non è infatti possibile distinguere da quale dei due impianti sia partito il trasporto e verso quale si sia diretto. Sono entrambi record del FER "55".

relazioni %>% 
  filter(O == D) %>% 
  select(Operatore = FER, 
         Origine = O, Destinazione = D, 
         Frontiera = Ft)

relazioni <- relazioni %>% filter(O != D)

Una query simile viene eseguita per isolare i trasporti non ambigui ma ridondanti nei punti di frontiera rispetto alle indicazioni di origine-destinazione. Gli operatori che ricadono in questa categoria sono il FER "16", con complessivamente 91 trasporti, e il FER "71", con 1 trasporto.
Per questi record viene eseguita l’elisione delle frontiere in sede di normalizzazione delle stringhe.

relazioni %>% 
  filter(O == Ft | D == Ft) %>% 
  group_by(Operatore = FER, 
           Origine = O, 
           Destinazione = D, 
           Frontiera = Ft
           ) %>% 
  summarize(`n. trasporti` = n())

Estrapolando le frontiere, oltre alle località già osservate nella prima rendicontazione, risulta un impianto inedito che tuttavia occorre 1 volta e non si considera strutturante per nuovi pattern di trasporto.

relazioni %>% 
  count(Frontiera = Ft, 
        name = "n. trasporti") %>% 
  slice_min(`n. trasporti`, n = 1)

"Vallorbe" è un punto di frontiera tra Francia e Svizzera è non è riconducibile a un unico omologo italiano. Sono tre i possibili valichi italiani che conducono a "Vallorbe": "Domo II", "Luino" e "Chiasso" o "Chiasso Sm.".
Viene perciò escluso questo trasporto, che corrisponde anche al solo del FER "31".
Con l’eliminazione del FER "31" il numero degli operatori scende a 35, di cui 34 beneficiari.

relazioni %>% 
  filter(Ft == "Vallorbe") %>% 
  select(Operatore = FER, 
         Origine = O, 
         Destinazione = D, 
         Frontiera = Ft)

relazioni <- relazioni %>% filter(FER != "31")

L’ultima operazione di controllo, prima di passare alla normalizzazione delle stringhe, viene effettuata sulle date per filtrare i trasporti conformi al periodo di incentivazione.
Il dato originale include una colonna NOTE, in cui sono contrassegnati i trasporti fuori periodo. Viene verificata la correttezza dell’attributo in colonna calcolando il range delle date rimanenti. Sono 15 i trasporti cancellati in questo passaggio, per un totale di 18 record su 39074 pervenuti.

library(lubridate)

relazioni <- relazioni %>% 
  filter(is.na(NOTE) | NOTE != "FUORI DAL PERIODO DI INCENTIVAZIONE")

relazioni %>% 
  summarize(`Periodo incentivazione` = date(min(data)) %--% date(max(data)))

A differenza della scorsa rendicontazione, la variabile data immagazzina correttamente la variabile tempo, dimostrandosi priva di anomalie.
Potendola sfruttare, si prevede una più agevole scomposizione in anno, mese e giorno per elaborazioni cumulative di futuro interesse, mentre in questo lavoro si sceglie di esplorare i trasporti avvenuti per quadrimestre del periodo, assimilando al primo dei tre l’unico giorno di agosto 2018.

relazioni <- relazioni %>% 
  mutate(anno = year(data), 
         mese = month(data), 
         giorno = day(data)
         ) %>% 
  select(- data) %>% 
  relocate(c(anno, mese, giorno), .after = FER)

Da una rapida verifica su tutte le variabili, viene appurata la completezza delle stringhe essenziali di origine-destinazione attraverso i 39056 record e, per contro, l’inadeguata copertura dei codici id, che idealmente consentirebbero di eludere lo sforzo di normalizzazione.

relazioni %>% 
  map_df(~ mean(is.na(.))) %>% 
  pivot_longer(everything(), 
               names_to = "Variabili", 
               values_to = "Not Available"
               ) %>% 
  arrange(desc(`Not Available`)) %>% 
  mutate(`Not Available` = scales::percent(`Not Available`))

Al termine dei prossimi passaggi, quello di normalizzazione e quello di sistematizzazione, sarà esportata la versione essenziale del dataset per le successive elaborazioni statistiche.

Elaborazioni

Le elaborazioni si articolano nelle fasi di normalizzazione, sistematizzazione e match del dataset.

Normalizzazione

Le dizioni RFI disponibili nel file "OD-RFI_località.csv" sono uniformi internamente e pertanto costituiscono una guida per la normalizzazione delle stringhe. Questo nella consapevolezza che nel giro di un anno le linee di forza del trasporto intermodale si mantengono relativamente stabili, salvo sconvolgimenti nei rapporti economici.

Tutte le dizioni RFI sono in maiuscolo ad eccezione di "Ortona", che è una località intermedia. La prima parte della procedura riguarda la conversione in maiuscolo in modo da restringere fin da subito le difformità di stile.

relazioni <- relazioni %>% 
  mutate(O_old = O, 
         D_old = D, 
         Ft_old = Ft
         ) %>% 
  mutate(across(c(O, D, Ft), toupper))

library(magrittr)

tibble(`n. stringhe originali` = relazioni %$% 
         coalesce(O_old, D_old, Ft_old) %>% 
         n_distinct(), 
       `n. stringhe maiuscolo` = relazioni %$% 
         coalesce(O, D, Ft) %>% 
         n_distinct())

Vengono create copie \\w+_old delle variabili di origine, destinazione e frontiera per raffronti col dato originale. Le origini, destinazioni e frontiere uniche sottoposte a prima conversione si riducono da 162 a 154.

Con la prima fase di normalizzazione, le combinazioni uniche tra operatore, origine, destinazione e frontiera ammontano a 355.

relazioni %>% 
  group_by(Operatore = FER, 
           Origine = O, 
           Destinazione = D, 
           Frontiera = Ft
           ) %>% 
  summarize() %>% 
  nrow()

Il calcolo delle combinazioni tra stringhe apre alla seconda e più profonda fase di normalizzazione.
Una volta oggettificate e scompattate per operatore, le combinazioni vengono esaminate in modo da predisporre le stringhe al più alto tasso di associazione con le origini-destinazioni RFI e il dato di sequenza di cui sono portatrici.
La procedura è funzionale anche alla successiva fase di sistematizzazione del dataset, quella in cui le origini/destinazioni estere sono proiettate in colonne apposite e sostituite dai punti di frontiera. L’operazione si rende necessaria perché la copertura delle liste RFI è nazionale e parte (o si interrompe) in corrispondenza delle frontiere.

Si riporta un esempio di combinazione calcolata per l’operatore "14" limitatamente ai primi 5 raggruppamenti di origine, destinazione e frontiera. È un caso rappresentativo di questa fase perché contiene sia nomi facilmente uniformabili che meritevoli di ponderazione ai fini del match.
Le dizioni "MELZO" e "GENOVA VOLTRI" vengono ricondotte agli equivalenti RFI "MELZO SCALO" e "GENOVA VOLTRI MARE". "LA SPEZIA" può invece corrispondere a stringhe RFI diverse benché afferenti allo stesso nodo logistico, "LA SPEZIA MIGLIARINA"/"LA SPEZIA MARITTIMA". Lo stesso vale per le frontiere, che nel novero degli RFI possono essere "CHIASSO"/"CHIASSO SM." e "MODANE"/"MODANE FOURNEAUX".

combinazioni <- relazioni %>% 
  group_by(Operatore = FER, 
           Origine = O, 
           Destinazione = D, 
           Frontiera = Ft, 
           .add = TRUE
           ) %>% 
  summarize(`n. trasporti` = n()) %>% 
  ungroup() %>% 
  group_split(Operatore)

combinazioni %>% pluck(5) %>% 
  head(5)

Viene lanciato il codice di normalizzazione, che cumula sul dataset di monitoraggio le modifiche per ciascuna delle combinazioni calcolate. Il risultato per la combinazione del FER "14", rispetto ai medesimi raggruppamenti, è il seguente. "MODANE" conserva la sua dizione, mentre "LA SPEZIA" e "CHIASSO" diventano "LA SPEZIA MARITTIMA" e "CHIASSO SM.".

combinazioni %>% pluck(5) %>% 
  head(5)

Sistematizzazione

La fase di sistematizzazione del dataset si compie riconoscendo le origini e destinazioni in Italia e all’estero.
L’obiettivo è passare da una struttura origine/destinazione/frontiera, con impianti italiani misti a impianti esteri, a una struttura che dissolve i punti di frontiera e separa le origini e destinazioni italiane ed estere.
In questo modo sarà possibile incrociare il dato di frequenza con il dato di sequenza del trasporto, rilevato entro i confini nazionali.

Per costruire il «flag» in-Italia, viene importato il file "nodi_RFI.shp" dei nodi ferroviari elaborato nell’ambito della rendicontazione 2017/18, completo di buona parte dei nuovi impianti e delle stringhe in formato RFI.
L’oggetto in_It è un vettore di stringhe corrispondenti ai nomi del dataset relazioni contenuti anche nello shape. Si appendono a in_It alcuni impianti non originariamente mappati.

library(sf)

nodi_rfi_old <- st_read("nodi_rfi_old.shp")

nomi_nodi <- as.character(nodi_rfi_old$nomirfi)

nomi_relazioni <- relazioni %$% 
  coalesce(O, D, Ft) %>% 
  unique()

in_It <- nomi_relazioni[nomi_relazioni %in% nomi_nodi] %>% 
  append(c("MONTIRONE", "JESI INTERPORTO", 
           "CHIGNOLO PO")) %>% 
  append(c("RUBIERA", "S.MARTINO CAVA", "ALESSANDRIA", 
           "POGGIO RUSCO", "MONFALCONE", 
           "TAVAZZANO", "UDINE", 
           "DINAZZANO", "RACCONIGI"))

L’uso condizionale del flag in-Italia permette di popolare le nuove colonne per le origini/destinazioni estere e di sostituirne correttamente i nomi con i punti di frontiera.
La sistematizzazione è così compiuta.

relazioni <- relazioni %>% 
  mutate(O_est = 
           case_when(!(O %in% in_It) ~ O, 
                     TRUE ~ NA_character_), 
         D_est = 
           case_when(!(D %in% in_It) ~ D, 
                     TRUE ~ NA_character_), 
         O = 
           case_when(!(O %in% in_It) & (D %in% in_It) ~ Ft, 
                     TRUE ~ O), 
         D = 
           case_when(!(D %in% in_It) & (O %in% in_It) ~ Ft, 
                     TRUE ~ D)
         ) %>% 
  select(FER, anno, 
         mese, giorno, 
         O_It = O, D_It = D, 
         O_est, D_est, 
         km_It, km_RL)

Parallelamente, vengono estratti e uniformati i nomi esteri di origine e destinazione nelle colonne appena create. Il processo di omogeneizzazione dei nomi esteri assume il municipio come unità di accorpamento unitario dei poli logistici. Le relazioni estere sono cioè identificate dai nomi delle città al netto delle specificazioni pur presenti nel dataset originale.

Similmente ai trasporti in Italia, le stringhe estere vengono adattate per garantire l’incrocio con le fonti di informazione geografica disponibili in formato tabellare o spaziale.

Tra le fonti, conviene il dataset tabellare delle città del mondo con almeno 1000 abitanti, provvisto di nomi già codificati, longitudine, latitudine e paese di appartenenza.

Il dataset contiene anche le città corrispondenti ai punti di frontiera italiani, di cui si duplicano i record nei casi di "CHIASSO"/"CHIASSO SM." e "VENTIMIGLIA"/"VENTIMIGLIA PARCO ROJA". Non necessita di questo passaggio "MODANE"/"MODANE FOURNEAUX". Nei casi di "BRENNERO", "DOMO II" e "TARVISIO BOSCOVERDE", si convertono gli originali "Brennero - Brenner", "Domodossola" e "Tarvisio".

# link: https://bit.ly/36LsRlS

file_move(path = "geonames-all-cities-with-a-population-1000.csv", 
          new_path = "geonames_cities.csv")

esteri <- read_csv2("geonames_cities.csv", 
                    col_types = cols(`Geoname ID` = col_character()))

Di tutte le città, si estraggono quelle nei fusi europei per evitare omonimie, specialmente con gli Stati Uniti. Le coordinate, inizialmente delle stringhe, vengono separate e ricondotte a numeri reali per la successiva fase di plotting.

esteri <- esteri %>% 
  janitor::clean_names() %>% 
  filter(str_detect(timezone, "Europe")) %>% 
  select(id = geoname_id, 
         città = ascii_name, 
         coordinate = coordinates, 
         paese = country
         ) %>% 
  separate(coordinate, 
           c("lat", "long"), 
           ","
           ) %>% 
  mutate(across(c(lat, long), as.numeric), 
         città = toupper(città))

Attraverso il «geocode» presente in tabella (colonna id), si rimuovono ulteriori duplicati europei, mentre in relazioni si rettificano i pochi nomi esteri misinterpretati come città. Tipicamente, sono nomi di quartieri interni ai municipi che influenzano la toponomastica delle stazioni.

L’oggetto relazioni, ora convenientemente trattato ai fini delle analisi, viene esportato nel formato di interscambio .csv.

relazioni %>% write_csv2("relazioniRL.csv")

Uno stralcio dimostrativo delle origini-destinazioni proiettate su territorio nazionale è il seguente.
L’estratto riguarda le prime 6 relazioni intermodali per n. trasporti. Sono le uniche relazioni con n. maggiore di 900. Gli operatori sono il FER "22", già prevalente nella prima rendicontazione, e il FER "23".

relazioni %>% 
  group_by(Operatore = FER, 
           `Origine italiana` = O_It, 
           `Destinazione italiana` = D_It
           ) %>% 
  summarize(`n. trasporti` = n()) %>% 
  ungroup() %>% 
  slice_max(`n. trasporti`, n = 6)

Match

L’operazione di match tra dato di frequenza e dato di sequenza del trasporto merci avviene estrapolando le origini-destinazioni italiane in relazioni per associarle alla lista di origini-destinazioni RFI.
Un match preliminare con il file della prima rendicontazione, "OD-RFI_località.csv", viene eseguito per isolare le origini-destinazioni scoperte e valutarne qualitativamente il recupero con le integrazioni manuali. L’esito delle integrazioni manuali è il file "OD-RFI_località2020.csv", che contiene un totale di 388 liste uniche rispetto alle 336 di partenza.
Altre integrazioni manuali sono state effettuate sul file "nodirfi2020.shp", sia per le origini-destinazioni che per i nodi ferroviari di attraversamento.

Si importano entrambi i file, il primo per il match tabellare, il secondo per il match spaziale. Dalle sequenze RFI organizzate per riga viene selezionata la colonna delle origini-destinazioni.

OD_rfi <- read_csv2("OD_rfi_località.csv", 
                    col_names = FALSE)

OD_rfi <- OD_rfi %>% select(OD_rfi = X1)

nodi_rfi <- st_read("nodi_rfi.shp")

Si verifica che le integrazioni coprono la quasi totalità dei trasporti. Le relazioni non coperte, o negative al match, riportano tutte un n. trasporti minore di 3.

relazioni %>% 
  transmute(`OD italiane` = paste(O_It, D_It, sep = "_")) %>% 
  group_by(`OD italiane`) %>% 
  summarize(`n. trasporti` = n()) %>% 
  filter(!`OD italiane` %in% OD_rfi$OD_rfi) %>% 
  arrange(desc(`n. trasporti`))

A partire dal distinguo positivo/negativo, si calcolano i dati di sintesi del match tabellare. Su 170 origini-destinazioni italiane uniche, risultano associate 157 contro 13 non coperte da RFI corrispondenti. All’alto tasso di match tra origini-destinazioni si allinea quello dei trasporti, che è del 99.95 %.
Viene definito anche un grado di «polverizzazione» pari a 1.38 trasporti per origine-destinazione non associata.
Nei flussogrammi a venire saranno mappati i 39038 trasporti positivi al match.

relazioni %>% 
  transmute(`OD italiane` = paste(O_It, D_It, sep = "_")) %>% 
  mutate(match = 
           case_when((`OD italiane` %in% OD_rfi$OD_rfi) ~ "positivo", 
                     TRUE ~ "negativo")
         ) %>% 
  group_by(match) %>% 
  summarize(`n. OD italiane` = 
              n_distinct(`OD italiane`), 
            `n. trasporti` = n(), 
            polverizzazione = 
              n() / n_distinct(`OD italiane`)
            ) %>% 
  mutate(`% OD italiane` = 
           100 * `n. OD italiane` / sum(`n. OD italiane`),
         `% trasporti` = 
           100 * `n. trasporti` / sum(`n. trasporti`)
         ) %>% 
  arrange(desc(match))

Summarize sui nodi ferroviari

Il match tra i dati di frequenza e sequenza dei trasporti in Italia è il fondamento della procedura di «summarize» sui nodi della rete ferroviaria. Questa fase di aggregazione numerica sintetizza il numero di attraversamenti per ciascun nodo delle 157 origini-destinazioni coperte.

Come osservato nell’ambito della prima annualità, le stringhe RFI denotano le località toccate dal trasporto, che in vari punti del percorso non sono riconducibili a nodi stazione, scalo o terminal della rete. Al contrario, le località generiche comprendono gli spazi tecnici di interscambio, deviazione e confluenza, i c.d. «punti di movimento» e i riferimenti geografici, frapponendosi ai nodi di effettivo interesse.

Per questa ragione, l’import dei nastri di località, distribuiti per riga e accostati alle 388 origini-destinazioni RFI, si accompagna al distinguo tra nodi in mappa e località generiche prima che venga effettuata l’associazione col n. trasporti.

Ciascuna riga in nastri_rfi, letta come list di stringhe unitarie, viene scompattata per mezzo del separatore ";" e trasformata in un vettore caratteri della propria lunghezza.

nastri_rfi <- read_lines("OD_rfi_località.csv")

nastri_rfi <- nastri_rfi %>% str_split(";")

Il flag in_mappa, similmente a quanto avvenuto con in_It, permette di scorporare da ogni vettore caratteri le stringhe cartograficamente presenti in nodi_rfi. Si appendono anche le origini-destinazioni RFI, in modo che siano mantenute nel successivo export del dato di sequenza al netto delle località generiche.

Viene creata contestualmente la funzione di setacciamento sift() per identificare e omettere tutti i nomi estranei all’oggetto in_mappa. Il flag viene così applicato simultaneamente al livello di ogni origine-destinazione.

in_mappa <- as.character(nodi_rfi$nomirfi) %>% append(OD_rfi$OD_rfi)

nastri_rfi <- nastri_rfi %>% 
  map_depth(2, ~ sift(., in_mappa)) %>% 
  map_depth(1, unlist) %>% 
  map(na.omit)

Le stringhe dei nodi ferroviari in sequenza tra ciascuna origine e destinazione vengono ricomposte col separatore e restituite nel file "nastri_rfi_net.csv".

nastri_rfi %>% 
  map_depth(1, ~ paste(., collapse = ";")) %>% 
  unlist() %>% 
  write_lines("nastri_rfi_net.csv")

Si effettua il join del n. trasporti rendicontati con l’elenco delle origini-destinazioni RFI, e se ne accosta l’ammontare ai relativi nastri. Sul totale di 388, si mantengono i 157 nastri aventi il n. trasporti, riconducibili alle altrettante origini-destinazioni che il dato di rendicontazione condivide con quello RFI.

OD_rfi_n <- OD_rfi %>% 
  left_join(relazioni %>% 
              group_by(OD_It = paste(O_It, 
                                     D_It, 
                                     sep = "_")
                       ) %>% 
              summarize(n_trasporti = n()), 
            by = c("OD_rfi" = "OD_It"))

nastri_rfi <- nastri_rfi %>% 
  cbind(OD_rfi_n$n_trasporti) %>% 
  as.data.frame()

names(nastri_rfi) <- c("stazioni", "n_trasporti")

nastri_rfi <- nastri_rfi %>% 
  mutate(n_trasporti = unlist(n_trasporti)) %>% 
  filter(!is.na(n_trasporti))

Avvalendosi del functional programming, viene lanciato il loop di associazione tra tutti gli elementi di stringa contenuti in ogni nastro e il n. trasporti corrispondente.

L’oggetto ottenuto, originariamente una list, viene dissolto e ricomposto in formato data.frame/tibble distinguendo in colonna le stazioni uniche e il summarize del n. trasporti. Dalle stazioni si escludono le origini-destinazioni, riconoscibili dal separatore "_". Ne risulta la tabella df_nodi, ancora sprovvista di coordinate.

df_nodi <- nastri_rfi %>% 
  mutate(stazioni = 
           map2(stazioni, 
                n_trasporti, 
                ~ paste(.x, .y, sep = ":")
                )
         ) %$% 
  unlist(stazioni) %>% 
  as_tibble() %>% 
  separate(value, 
           c("stazioni", "n_trasporti"), 
           sep = ":"
           ) %>% 
  filter(!str_detect(stazioni, "_")) %>% 
  mutate(n_trasporti = 
           as.numeric(n_trasporti)) %>% 
  group_by(stazioni) %>% 
  summarize(n_trasporti = 
              sum(n_trasporti)) %>% 
  ungroup() %>% 
  arrange(desc(n_trasporti))

Viene eseguito un join tra lo shapefile dei nodi ferroviari, contenuto in nodi_rfi, e df_nodi per operare una scrematura dei nodi con almeno un trasporto, pari a 933. Si riassegna a df_nodi lo stesso contenuto di nodi_rfi, resosi equivalente, ma con l’aggiunta delle coordinate di longitudine e latitudine in colonne distinte.

Si sceglie di scorporare i valori di coordinate e accantonare la colonna di informazione geometrica, propria di uno shapefile, per due motivi interdipendenti:

df_nodi <- nodi_rfi %>% 
  select(stazioni = nomirfi, 
         n_trasporti, 
         OD_It = filtro_OD, 
         regioni = reg
         ) %>% 
  cbind(st_coordinates(.)) %>% 
  as_tibble() %>% 
  rename(long = X, 
         lat = Y
         ) %>% 
  select(- geometry) %>% 
  arrange(desc(n_trasporti))

Si riporta di seguito uno stralcio del summarize sui nodi con le variabili rinominate per la lettura. Si estraggono, per regione e flag dei nodi di origine-destinazione, le prime stazioni per n. trasporti. Le regioni selezionate sono Lombardia ("3"), Piemonte ("1") e Liguria ("7").

In Lombardia, tra le stazioni che oltre a essere attraversate sono di origine-destinazione, quella più frequentata è "GALLARATE", mentre spicca tra le stazioni intermedie quella di "MILANO LAMBRATE". Esistono tre ex-aequo in Piemonte e due in Liguria tra le stazioni di solo attraversamento.

df_nodi %>% 
  rename(Stazioni = stazioni, 
         `n. trasporti` = n_trasporti, 
         `OD italiane` = OD_It, 
         Regioni = regioni, 
         Longitudine = long, 
         Latitudine = lat
         ) %>% 
  group_by(Regioni, 
           `OD italiane`) %>% 
  filter(Regioni %in% c("1", "3", "7")) %>% 
  slice_max(`n. trasporti`, n = 1)

Si esporta il .csv del summarize operato sui nodi, completo delle coordinate di longitudine e latitudine.

df_nodi %>% write_csv2("df_nodi.csv")

Si esporta anche la conversione in classe sf («simple feature»/shapefile), proiettando l’oggetto sf_nodi nel file "sf_nodi.shp". A differenza dell’originale "nodi_rfi.shp", questo file contiene i 933 nodi unici attraversati nell’anno, il n. trasporti aggiornato e le coordinate distinte in longitudine e latitudine. Prima di esportare, si imposta per l’output lo stesso sistema di riferimento (crs) dell’input.

sf_nodi <- st_as_sf(df_nodi, 
                    coords = c("long", "lat"), 
                    remove = FALSE
                    ) %>% 
  st_set_crs(st_crs(nodi_rfi))

sf_nodi %>% st_write("sf_nodi.shp")

Generazione del flussogramma

Per realizzare il flussogramma, si ripercorre l’associazione tra origini-destinazioni italiane, n. trasporti, e nastri delle sole località presenti in mappa, cioè le stringhe dei nodi stazione/scalo/terminal già prelevate con sift() degli RFI originali.

Una differenza è il contenuto dell’oggetto in_mappa, da cui si rimuovono le stringhe di origine-destinazione, le prime di ogni nastro. Serve per non confonderle con quelle che si otterrebbero concatenando con lo stesso carattere "_" i nomi consecutivi tra origine e destinazione. Oltre al sift(), rieseguito sui nastri per escludere le origini-destinazioni, viene così creata la funzione concatenate(), che incolla a due a due le stringhe di ciascuna sequenza.
Il risultato viene incolonnato assieme al n. trasporti corrispondente. Una volta mantenuti i record con almeno 1 trasporto, viene lanciato un secondo loop, questa volta tra n. trasporti e ogni coppia di stringhe. Le coppie unificano i nomi dei nodi consecutivi in formato "partenza_arrivo". Sono necessarie a segmentare i percorsi univocamente per restituire in cartografia l’armatura dei flussi.

Si ricava da questo ciclo df_archi, con separazione di coppie e n. trasporti in variabili distinte alla maniera del df_nodi. Dopo il summarize dei trasporti per arco, o coppia unica di nodi, si provvede a scomporre le stringhe in nomi dei nodi di partenza e arrivo.

L’operazione consente il doppio join con df_nodi, che per ogni partenza e arrivo riporta i dati di longitudine, latitudine, flag origine-destinazione, e regione di appartenenza. Qui si esprime l’importanza del df_nodi senza la «sticky column» delle coordinate geografiche, che sarebbe altrimenti da cancellare in ogni passaggio.

I due join, uno sulle partenze e l’altro sugli arrivi, si svolgono elidendo il n. trasporti dei nodi e rinominando coerentemente gli attributi.

df_archi <- df_archi %>% 
  left_join(df_nodi %>% 
              select(- n_trasporti), 
            by = c("partenza" = "stazioni")
            ) %>% 
  rename_with(~ str_c(., "_partenza"), OD_It:lat) %>% 
  left_join(df_nodi %>% 
              select(- n_trasporti), 
            by = c("arrivo" = "stazioni")
            ) %>% 
  rename_with(~ str_c(., "_arrivo"), OD_It:lat)

Si esegue a scopo dimostrativo una query simile a quella costruita sui nodi per visualizzare uno stralcio di df_archi. Sul totale di 1905, vengono estratti gli archi col maggior n. trasporti tra i gruppi interregionali e intraregionali di Lombardia ("3"), Piemonte ("1") e Liguria ("7").

Il nodo di "VOGHERA" appartiene a un arco che registra il più alto numero di attraversamenti tra quelli sul confine Lombardia-Piemonte. "ARQUATA SCRIVIA" è un nodo dell’arco che prevale tra quelli sul confine Piemonte-Liguria. Tre gli ex aequo intraregionali in Lombardia e Piemonte, che comprendono le frontiere "DOMO II" e "LUINO".

df_archi %>% 
  rename(`Stazione partenza` = partenza, 
         `Stazione arrivo` = arrivo, 
         `n. trasporti` = n_trasporti, 
         `OD italiane partenza` = OD_It_partenza, 
         `Regioni partenza` = regioni_partenza, 
         `Longitudine partenza` = long_partenza, 
         `Latitudine partenza` = lat_partenza, 
         `OD italiane arrivo` = OD_It_arrivo, 
         `Regioni arrivo` = regioni_arrivo, 
         `Longitudine arrivo` = long_arrivo, 
         `Latitudine arrivo` = lat_arrivo
         ) %>% 
  group_by(`Regioni partenza`, 
           `Regioni arrivo`
           ) %>% 
  filter(`Regioni partenza` %in% 
           c("1", "3", "7") & 
           `Regioni arrivo` %in% 
           c("1", "3", "7")
         ) %>% 
  slice_max(`n. trasporti`, n = 1)

df_archi %>% write_csv2("df_archi.csv")

A questo punto è possibile processare le variabili di longitudine e latitudine in partenza e arrivo per trarne il dato geografico non più come POINT ma come LINESTRING.

Popolata la colonna geometry, si converte df_archi nell’oggetto spaziale sf_archi, che viene di seguito esportato in .shp come fatto per df_nodi/sf_nodi.

sf_archi <- st_as_sf(df_archi) %>% st_set_crs(st_crs(nodi_rfi))

sf_archi %>% st_write("sf_archi.shp")

Il risultato cartografico è il flussogramma dei trasporti complessivi nell’anno 2018/19.

Pattern di relazioni nel tempo

Il portato informativo di relazioni rende possibile esplorare l’impegno ferroviario per operatore, periodo e origine-destinazione estera. Questi dati aprono a una definizione interrogabile o sfaccettabile del flussogramma.

Relazioni in Italia

Il primo passaggio è il recupero delle stringhe di origine-destinazione in nastri_rfi e coppie_rfi, oggetti dotati delle 157 liste di stazioni/scali/terminal attraversati almeno 1 volta nell’anno, espressi in sequenza come nomi singoli e a due a due.

Le stringhe di origine-destinazione servono anche in questo caso per il match con le frequenze di trasporto, che a differenza del flussogramma complessivo sono disaggregate per operatore, anno e mese.

L’oggetto relazioni_info è interprete delle combinazioni uniche tra queste variabili, una procedura anticipata in fase di normalizzazione delle stringhe ma senza considerare l’anno e il mese di trasporto. La struttura dei dati è tale che l’uso di count() (v. chunk sotto) è alternativo al group_by() %>% summarize() del n. trasporti per combinazione.

Le frequenze vengono ordinate cronologicamente dalla più alta alla più bassa, e viene eseguito il doppio join con nastri_rfi e coppie_rfi sul campo comune delle origini-destinazioni.

relazioni_info <- relazioni %>% 
  count(OD_It = paste(O_It, D_It, sep = "_"), 
        FER, 
        anno,
        mese, 
        name = "n_trasporti"
        ) %>% 
  arrange(anno, 
          mese, 
          desc(n_trasporti)
          ) %>% 
  right_join(nastri_rfi %>% 
               select(- n_trasporti), 
             by = c("OD_It" = "OD_rfi")
             ) %>% 
  na.omit() %>% 
  left_join(coppie_rfi %>% 
              select(- n_trasporti), 
            by = c("OD_It" = "OD_rfi"))

La tabella relazioni_info si compone dunque in c.d. «formato lungo», con le due «list-column» di stazioni singole e a coppie in sequenza, e i dati di operatore, anno, mese con rispettivo n. trasporti.

La struttura in formato lungo è compatibile con la Grammar of Graphics di ggplot2 e permette di operare query sulla somma dei trasporti, sfaccettandoli in base a condizione di interesse (es. n. trasporti del FER 10 nel luglio 2019).

Su relazioni_info vengono compiute le medesime procedure di functional programming, con l’accortezza di abbinare l’insieme degli attributi di trasporto a ciascuna stringa. Per farlo, si uniscono gli attributi FER, anno, mese e n_trasporti nella nuova variabile info_trasporti con separatore ":".

Si procede con il loop per associare info_trasporti alle stazioni/coppie avendo cura che, una volta incolonnati gli attributi per ":", siano ripristinati come numerici quelli di anno, mese e n. trasporti. Dopodiché si opera la somma dei trasporti per tutte le altre variabili (stazioni/coppe, operatore e periodo).

Il join con df_nodi e df_archi si esegue per ereditare i valori di longitudine, latitudine e n. trasporti complessivo per nodo/arco, che viene così affiancato al n. trasporti parziale. Da questa serie di operazioni si ottengono df_nodi_info e df_archi_info, con le copie spaziali sf_nodi_info e sf_archi_info per la resa cartografica.

Agli oggetti df/sf_\\w+_info si decide di applicare una trasformazione simultanea per aggiungere e collocare la nuova variabile quadrimestre. Il 31-08-2018 viene assimilato al primo quadrimestre.

library(zeallot)

c(df_nodi_info, 
  df_archi_info, 
  sf_nodi_info, 
  sf_archi_info
  ) %<-% 
  map(list(df_nodi_info, 
           df_archi_info, 
           sf_nodi_info, 
           sf_archi_info
           ), 
      ~ mutate(., 
               quadrimestre = case_when(anno == 2018 ~ 1, 
                                        anno == 2019 & mese %in% 1:4 ~ 2, 
                                        TRUE ~ 3)
               ) %>% 
        relocate(., 1:mese, quadrimestre))

Con gli attributi di trasporto, il flussogramma diventa un modello diacronico (quadrimestre) e sincronico (FER). Un esempio di questa rappresentazione isola i due operatori principali per n. trasporti, FER "22" e "58", accomunando i minori in categoria residuale. I pattern si prestano a osservazioni che non si colgono nel flussogramma generale.

Il FER "22" in Italia non si spinge oltre le regioni del Nord, e si addensa tra la direttrice del Sempione e i valichi alpini lombardo e piemontese. Solo nel terzo quadrimestre l’operatore si rafforza in Emilia e inaugura il prolungamento dei trasporti allacciandosi a Veneto e Friuli.

Al contrario, il FER "58" si sviluppa in tutta la penisola, ma è discontinuo in Lombardia e la sua attività si concentra in Lombardia e sul Tirreno nell’ultimo periodo dell’anno.

Nel complesso degli altri operatori, si evidenzia la stabilità della regione logistica milanese e delle rotte adriatica e tirrenica, con un incremento del volume di traffico appenninico a partire dal secondo quadrimestre.

Relazioni estere

Per mappare questa componente di relazioni si adopera l’oggetto esteri, già predisposto in sede di normalizzazione e comprensivo dei nomi di città europee con oltre 1000 abitanti, longitudine, latitudine e paese. Attraverso la normalizzazione, si verifica l’esatta corrispondenza tra l’insieme di frontiere italiane-recapiti esteri in relazioni, e le città con i dati di georiferimento in esteri.

Con la sistematizzazione degli impianti italiani ed esteri, passaggio chiave del match con RFI, l’oggetto relazioni presuppone ora il recupero dei punti di frontiera originali per estrarre il complemento dei trasporti in Italia. Sono funzionali allo scopo i valori NA impostati nelle origini/destinazioni estere, in quanto tracce lasciate dagli ex-punti di frontiera nei rispettivi record.

In questo modo relazioni si trasfigura in relazioni_est, esito anche di aggregazione dei trasporti per gli attributi operatore e periodo. Nella procedura vengono eliminate le colonne dei trasporti in Italia e si provvede ai join con i dati geografici di esteri sulle nuove origini-destinazioni, avendo cura di esplicitare con nomi adeguati le variabili di longitudine, latitudine e paese.

relazioni_est <- relazioni %>% 
  count(FER, 
        anno, mese, 
        O_It, D_It, 
        O_est, D_est, 
        name = "n_trasporti", 
        sort = TRUE
        ) %>% 
  filter(!is.na(O_est) | !is.na(D_est)) %>% 
  mutate(O_est = case_when(is.na(O_est) ~ D_It, 
                           TRUE ~ O_est), 
         D_est = case_when(is.na(D_est) ~ O_It, 
                           TRUE ~ D_est)
         ) %>% 
  select(!ends_with("_It")) %>% 
  left_join(esteri %>% 
              rename_with(~ str_c(., "_O_est"), 
                          - città), 
            by = c("O_est" = "città")
            ) %>% 
  left_join(esteri %>% 
              rename_with(~ str_c(., "_D_est"), 
                          - città), 
            by = c("D_est" = "città"))

Anche nelle relazioni estere di decide di introdurre il quadrimestre come unità di analisi temporale. Dopodiché viene costruita la variabile geometry con lo stesso procedimento di segmentazione in archi, questa volta intesi come linee di desiderio (geometricamente delle rette) tra frontiere italiane e recapiti esteri.

Parallelamente, viene scaricato e assegnato a paesi_est lo shapefile dei limiti nazionali europei per la base cartografica. Di questo layer spaziale vengono preservati i nomi internazionali dei paesi, in inglese come nel dataset esteri, sincerandosi che il crs delle loro geometrie POLYGON sia lo stesso dei POINT e LINESTRING negli altri oggetti spaziali prodotti.

# link: https://bit.ly/3ltn8VA

paesi_est <- st_read("Europe/Europe.shp")

paesi_est <- paesi_est %>% 
  janitor::clean_names() %>% 
  mutate(name = 
           case_when(name == 
                       "Bosnia Herzegovina" 
                     ~ "Bosnia and Herzegovina", 
                     TRUE ~ name)
         ) %>% 
  select(paese_est = name)

Esplorando i paesi di origine e destinazione con oltre 10 trasporti, i volumi di traffico per quadrimestre nell’anno appaiono costanti, e privilegiano le avanzate piattaforme logistiche di Germania e Benelux. Tra i partner più rilevanti anche la Svizzera, mentre sono più rarefatti i flussi di Spagna ed Est Europa.

relazioni_est %>% 
  group_by(Quadrimestre = quadrimestre, 
           `Paese origine` = paese_O_est, 
           `Paese destinazione` = paese_D_est
           ) %>% 
  summarize(`n. trasporti` = sum(n_trasporti), .groups = "drop") %>% 
  group_split(D_est = `Paese destinazione` != "Italy", 
              O_est = `Paese origine` != "Italy") %>% 
  map(~ select(., 
               where(~ sum(str_detect(., "Italy")) == 0), 
               - O_est, 
               - D_est
               ) %>% 
        filter(., `n. trasporti` >= 10) %>% 
        pivot_wider(., 
                    names_from = Quadrimestre, 
                    values_from = `n. trasporti`, 
                    values_fill = 0
                    ) %>% 
        relocate(., 
                 c(`1`, `2`, `3`), .after = last_col())
      ) %>% 
  map_at(1, 
         ~ rename_with(., ~ str_c("O q", .), - 1)) %>% 
  map_at(2, 
         ~ rename_with(., ~ str_c("D q", .), - 1)) %>% 
  map(~ rename(., Paese = 1)) %>% 
  reduce(inner_join) %>% 
  arrange(across(- Paese, desc))

Non si ravvisano grandi mutamenti attraverso le soglie quadrimestrali. Nel complesso degli operatori minori per n. trasporti, avviene dal secondo quadrimestre un infittimento dei rapporti con l’Est Europa. Il FER "22", tra i due operatori maggiori, intrattiene rapporti stabili con più paesi del corridoio centrale, approdando anche in Danimarca e Spagna. Il FER "55", orientato come gli altri verso il Mare del Nord, alimenta nel tempo il volume di traffico, che cresce nel secondo e terzo quadrimestre.

R modello di lavoro

Microsoft Excel vs. R

A ogni strumento si deve un giudizio senza ideologie, relativamente al solo scopo per cui è concepito. Microsoft Excel è l’applicativo statistico che convenzionalmente si utilizza per maneggiare database alfanumerici, filtrando ed eventualmente ricomponendo le variabili rispetto alla sintesi delle osservazioni.

Tuttavia, la naturale centratura di Excel su una griglia di celle predispone lo strumento ad alterazioni manuali, rendendolo strutturalmente adeguato all’attività compilativa piuttosto che al governo di progetti statistici complessi.

R opera invece in astratto, nell’accezione di Object-Oriented Programming: la tabella è solo una delle modalità di inscrizione dei dati, e ciascuna cella potrebbe persino custodire a sua volta un’ulteriore tabella.

Efficienza

Fondandosi su un linguaggio, R vanta velocità e capacità di calcolo enormemente superiori se paragonato a Excel. A parità di (buona) memoria, il tasso di efficienza di Excel rispetto a R è di un’ora per migliaia di record contro pochi secondi per milioni.

Per quanto lo stesso Excel sappia aprire prospetti con numerosissime osservazioni, la particolare costruzione a reticolo di celle indipendenti, diversamente dalla linea di comando di R, implica che la scrittura di una formula sia costantemente ripetuta assieme ai suoi riferimenti incrociati. Questo fa la differenza quando si generano tabelle pivot o più semplicemente nuove variabili: in corrispondenza di ciascuna delle colonne estraibili dalle esistenti, se il prospetto contiene 1 milione di record il calcolo deve essere incollato e rieseguito 1 milione di volte, mentre in R tutto questo è vettorizzato in un unico comando.

Eleganza

Ne consegue l’abbattimento dei rischi di inesattezza e perdita dei dati dovuti alla meccanicità delle procedure di Excel, che in R sono preservate da script riproducibili integralmente. Tutto il percorso di elaborazione è rieseguibile dall’inizio alla fine, il che nel lungo periodo agevola il necessario lavoro di riepilogo e reportistica, anche nelle fasi intermedie. Parallela alla riproducibilità è la specializzazione disciplinare di R, con risorse modellistiche declinate attraverso i circa 20 mila pacchetti disponibili e in aumento.

La rappresentazione grafica è altrettanto conveniente. Le visualizzazioni in R sono più professionali con la Grammar of Graphics di ggplot2. Soprattutto, possono essere esportate in .html accessibili dal web senza perdere l’eventuale interattività o animazione.

Collaborazione

R si distingue per questi motivi anche sull’ultimo fronte inerente al modello di lavoro, quello della collaborazione. La prassi in Excel tipicamente consiste nel rimbalzo di e-mail inoltrate, un aspetto ricorrente che non facilita l’elemento umano quando si tratta di riconsolidare il percorso delle modifiche nel tempo.

Emergono con R due strade. Una è RStudio Connect, abbonamento per lo scambio e la pubblicazione dei progetti statistici dove è possibile impostare per uno o più utenti i filtri di interazione col codice e i risultati. L’altra strada, fedele al comune spirito open-source, è la condivisione su GitHub, portale online pensato per costruire collegialmente i progetti di data science e ingegneria del software.

Prospettiva

L’evoluzione di R per il data science, l’esempio di report statistico e il raffronto col classico Excel lasciano intuire l’innovazione di processo che R può introdurre e fertilizzare.

Lo scenario auspicabile con Tidyverse e Markdown è l’approdo a un sistema in sharing di progetti statistici programmabili, dove l’utente di una unità di ricerca riesce ad accedere alla «scatola nera» di un report prodotto in altre unità, e a replicarne le operazioni cumulando nuova conoscenza sugli esiti.

Si prevede attraverso l’upgrade delle attività statistiche un migliore coordinamento dei prodotti di ricerca, soprattutto quelli compositi e transdisciplinari come Rapporto Lombardia, che già si presta concettualmente a un grande Markdown, e un apprendimento diffuso del nesso tra dati originali e informazioni che saranno poi al servizio dei policy maker regionali.

Riferimenti

PoliS-Lombardia (2019). Studi di accompagnamento alle misure di Dote Merci Ferroviaria regionale. (Codice 190717TER). Milano, PoliS-Lombardia.

R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. https://www.R-project.org/.

Programmazione statistica con R: un esempio sui dati del trasporto intermodale in Lombardia