Elaborazione dei dati meteo con R

Elaborazione dei dati meteo con R

Elaborazione dei dati meteo con R

<!DOCTYPE html>



La trattazione dei dati meteorologici con R

Introduzione

Il questo post vedremo come è possibile trattare ed elaborare dati meteorologici provenienti da capannine di meteo. Il software che verrà utilizzato sarà R, disponibile al sito: https://www.r-project.org/; un software statistico molto potente e open source.

In R esistono diverse librerie che sono in grado di gestire dati territoriali e dati spaziali molto potenti. In questo post utilizzeremo la libreria spacetime che utilizza la libreria sp per la gestione dei dati spaziali e la libreria zoo e xts per il dato temporale.

Importazione dei dati

Importazioni dati meteo

I dati delle capannine meteo sono registrati in un database PostgreSQL. L’accesso da R al database è garantito dalla libreria RPostgreSQL

library("RPostgreSQL", lib.loc="~/R/i686-pc-linux-gnu-library/3.2")
# Si stabilisce la connessione a PoststgreSQL con RPostgreSQL
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, dbname="codipg")# Versione semplificata (database su localhost)
# Si esegue la query per importare i dati necessari
dati_pioggie <- dbGetQuery(con, "SELECT nome, data, mm FROM pioggia2 order by 2, 1")
#converto in fattore il nome delle stazioni
dati_pioggie$nome = factor(dati_pioggie$nome)

Vediamo come è strutturato il dato meteo

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
head(dati_pioggie)
##                nome       data  mm
## 1          Allerona 1992-01-01 0.0
## 2          Bastardo 1992-01-01 0.1
## 3        Casigliano 1992-01-01 0.0
## 4           Cerbara 1992-01-01 0.1
## 5 Citta_di_Castello 1992-01-01 0.0
## 6           Corbara 1992-01-01 0.2

Importazione dati spaziali

Adesso risulta necessario importare lo shape puntuale che individua la localizza.

stazioni =  readOGR(dsn=".","stazioni_point")
row.names(stazioni) = as.character(stazioni$nome)

Abbiamo così creato uno SpatialPointDataFrame ovvero l’analogo di uno shape puntuale. Ne controlliamo la tipologia e ne facciamo un sommario:

class(stazioni)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
summary(stazioni)
## Object of class SpatialPointsDataFrame
## Coordinates:
##               min     max
## coords.x1 2272679 2320496
## coords.x2 4698989 4821088
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0
## +ellps=intl +units=m +no_defs]
## Number of points: 30
## Data attributes:
##                 nome   
##  Allerona         : 1  
##  Amelia           : 1  
##  Attigliano       : 1  
##  Avigliano_Umbro  : 1  
##  Bastardo         : 1  
##  Calvi_dell_Umbria: 1  
##  (Other)          :24

Noteremo che esiste una sola colonna attributi chiamata nome

Collegamento dati spaziali e dati meteo

Utilizzeremo la libreria reshape per avere una tabella dati in forma lunga Per apprfondimenti.

library("reshape", lib.loc="~/R/i686-pc-linux-gnu-library/3.2")

capannine=2:31
# faccio il cast per assicurare la struttura Long della tabella
dati_pioggieLong=cast(dati_pioggieMelt,data~nome,value = "value", fill = NA)

Possiamo adesso accoppiare il dato spaziale ad dato temporale in forma long. Utilizziamo la funzione STFDF che crea uno SpaceTimeFullDataFrame ovvero un oggetto spazio tempo di tipo full. Ovvero nei giorni in cui non è presente il dato meteorologico inserisce un NA ovvero segnala il dato mancate.

#Costruisco un oggetto spazio tempo
pioggieST1=STFDF(stazioni, dati_pioggieLong$data,data.frame(values = as.vector(t(dati_pioggieLong[capannine]))))

Diamo uno sguardo a questo nuovo oggetto:

summary(pioggieST1)
## Object of class STFDF
##  with Dimensions (s, t, attr): (30, 8036, 1)
## [[Spatial:]]
## Object of class SpatialPointsDataFrame
## Coordinates:
##               min     max
## coords.x1 2272679 2320496
## coords.x2 4698989 4821088
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0
## +ellps=intl +units=m +no_defs]
## Number of points: 30
## Data attributes:
##                 nome   
##  Allerona         : 1  
##  Amelia           : 1  
##  Attigliano       : 1  
##  Avigliano_Umbro  : 1  
##  Bastardo         : 1  
##  Calvi_dell_Umbria: 1  
##  (Other)          :24  
## [[Temporal:]]
##      Index              timeIndex   
##  Min.   :1992-01-01   Min.   :   1  
##  1st Qu.:1997-07-01   1st Qu.:2010  
##  Median :2002-12-31   Median :4018  
##  Mean   :2002-12-31   Mean   :4018  
##  3rd Qu.:2008-07-01   3rd Qu.:6027  
##  Max.   :2013-12-31   Max.   :8036  
## [[Data attributes:]]
##      values      
##  Min.   :  0.00  
##  1st Qu.:  0.00  
##  Median :  0.00  
##  Mean   :  2.29  
##  3rd Qu.:  0.60  
##  Max.   :196.80  
##  NA's   :82980

Noteremo che esistono 30 punti che sono le 30 capannine di monitoraggio e dati temporali che iniziano il 1-1-1992 e finiscono il 31-12-2013.

I dati presenti in questo oggetto possono essere acquisiti utilizzando la sintassi con la parentesi [ ] nella forma pioggieST1[i,j,k] dove i si riferisce all’i-esima stazione, j al tempo e k alla variabile. Nel nostro caso abbiamo una sola variabile denominata mm che riferisce alla cumulata giornaliera espressa in millimetri di pioggia. La sintassi seguente fornisce ad esempio il dato per la stazione di Amelia nel periodo dal 10-10-2012 al 15-10-2012

pioggieST1["Amelia","2012-10-10/2012-10-15"]
##            values timeIndex
## 2012-10-10    0.2      7589
## 2012-10-11   14.2      7590
## 2012-10-12   65.0      7591
## 2012-10-13   28.2      7592
## 2012-10-14    0.2      7593
## 2012-10-15   26.8      7594

Trattazione dei dati

La funzione aggregate permette appunto l’aggregazione dei dati.

#aggrego i dati annualmente richiedendo il dato massimo annuo
pioggie_annuali = aggregate(pioggieST1, "year", max)
#mensilmente (dato massimo)
pioggie_mensili = aggregate(pioggieST1, "month", max)
#mensilmente dato medio
pioggie_mensili_max = aggregate(pioggieST1, "month", mean)

Grafici degli andamenti

stplot(pioggie_annuali,mode="tp",main="Andamento pluviometrico del massimo annuo",xlab='anni',ylab='mm') # plot degli andamenti delle singole stazioni

stplot(pioggie_annuali,mode="ts",main="Andamenti massimi annui",xlab='anni',ylab='mm') #plot di tutte le stazioni in un unico grafico, si nota come le stazioni hanno un andamento concorde annualmente