-
Notifications
You must be signed in to change notification settings - Fork 0
/
creaClimatolData.R
191 lines (129 loc) · 7.84 KB
/
creaClimatolData.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#3 ottobre 2017
#Prepara i dati nel formato richiesto da climatol per l'omogeneizzazione delle serie
rm(list=objects())
library("purrr")
library("readr")
library("broom")
library("stringr")
library("imputeTS")
source("ClimateObjects.R")
source("ClimateData.R")
source("climatol_help.R")
# Funzioni utili per il programma -----------------------------------------
# lettura anagrafica
leggiAna<-function(nomeFile){
#lettura del file con le stazioni da processare
tryCatch({
read_delim(file=nomeFile,delim=";",col_names=TRUE)
},error=function(e){
stop(sprintf("Errore lettura file anagrafica %s",nomeFile))
})
}#fine leggiAna
# Funzione riduci ---------------------------------------------------------
#applica reduce distinguendo tra dati giornalieri e mensili
riduci<-function(listaOut,monthlyJoin=TRUE){
if(monthlyJoin){
c("yy"="yy","mm"="mm")->colJoin
}else{
c("yy"="yy","mm"="mm","dd")->colJoin
}
listaOut %>% reduce(full_join,by=colJoin)
}#fine riduci
# Funzione per riempire gli NA mensili: da evitare ------------------------
# funzione per riempireNA -------------------------------------------------
riempiNA<-function(x){
# Identificazione anni in cui tutte le stazioni sono NA -------------------
#somma va applicata solo alle colonne dei dati (siano essere serie mensili o giornaliere)
apply(x[,3:ncol(x)],1,FUN =function(val){
length(val)->numeroDati
which(is.na(val))->qualiNA
length(qualiNA)->numeroNA
#cerchiamo solo le righe (anni) in cui almeno metà delle serie ha NA,sono gli anni che vanno
#interpolati perchè altrimenti climatol si blocca
if(numeroNA>floor(numeroDati/2)) return(NA)
sum(val,na.rm=TRUE)
})->vettoreSomma
#indexNA identifica i mesi all'interno di un anno che sono NA in tutte le serie in esame.
which(is.na(vettoreSomma))->indexNA
#per interpolare usiamo il pacchetto imputeTS: interpoliamo solo le righe (gli anni di ciascuna serie) "problematiche"
#(ovvero quelle identificate dalla funzione somma)
map_df(daScrivere,.f=function(serie){
for(ii in indexNA){
if(is.na(serie[ii])){
na.interpolation(serie,option="linear")->serieInterpolata
serie[ii]<-serieInterpolata[ii]
}
}#fine ciclo for
serie
})
} #fine riempiNA
# INIZIO PROGRAMMA --------------------------------------------------------
REGIONE<-"Cluster4" #potrebbe anche essere regione Lazio...
FILE_ANAGRAFICA<-"anagrafica_cluster_4_tmax.csv"
#directory in cui si trovano le varie subdirectory con le serie giornaliere da elaborare.
#Si tratta (ad esempio) delle directory in cui sono stati eseguiti i controlli spaziali. Ogni directory ha
#il nome di una regione italiana (compresa Aeronautica) e contiene i file.txt (ad esempio i file prodotti dai controlli spaziali)
DIRSERIE<-"spatial_controls"
#Quindi dentro DIRSERIE si trovano le varie subdirectory con le serie da elaborare.Perchè una struttura così complessa e non mettere i file tutti in una directory?
#Perchè le serie sono identificate da un numero id che non è unico "123456.txt". Se mettessimo tutto in una directory i file verrebbero sovrascritti e perderemmo
#l'informazione sulla fonte del dato (da che regione proviene?).
PARAM<-c("tmax","tmin")[1]
ANNOI<-1961
ANNOF<-2015
#Se vogliamo in output serie mensili:MONTHLY è TRUE, se vogliamo serie giornaliere allora MONTHLY è FALSE
MONTHLY<-FALSE
#Se vogliamo riempire gli NA mensili quando nessuna stazione ha valori disponibili dobbiamo porre: MONTHLY<-TRUE && REFILL_NA_MONTHLY<-TRUE
#In generale: evitare riempir perche si potrebbero riempire anni conuno stesso identico valore
REFILL_NA_MONTHLY<-FALSE
#Crea grafico serie mensili? Sia che si producano in output le serie giornaliere che quelle mensili, un grafico delle serie
#mensili (o annuali) può aiutare a capire i risultati dell'omogeneizzazione
GRAFICO_SERIE_ANNUALI<-TRUE
#overloading della funzione
purrr::partial(...f=creaClimatolData,annoI=ANNOI,annoF=ANNOF)->creaClimatolData_ii_ff
#la colonna percorsoSerie contiene la posizione (directy) della serie in anagrafica
#Facciamo affidamento che le subdirectory in DIRSERIE abbiano un nome analogo a quello che compare in anagrafica (Regione)
leggiAna(nomeFile = FILE_ANAGRAFICA)->ana
#crea il percorso di directory dove trovare il dato giornaliero
purrr::map_chr(ana$regione,.f=function(rr){
list.files(pattern=paste0("^spatial_controls_",rr,"_.+$"),include.dirs = TRUE)->dirDati
stopifnot(length(dirDati)==1)
dirDati
})->listaDirDati
ana$percorsoSerie<-paste0("./",listaDirDati,"/",ana$SiteId,".txt")
# %>%
# mutate(percorsoSerie=paste0(DIRSERIE,"/",regione,"/",SiteId,".txt"))->ana
#listaOut contiene serie mensili o serie giornaliere. Sono serie abbastanza continue e complete (vedere codice). In pratica abbiamo filtrato
#le serie, eliminando serie con troppi buchi
#Marzo 2018: siccome già a priori le serie sono state controllarte in termini di completezza, dentro creaClimatolData_ii_ff mttiamo dei controlli molto laschi in modo
#di non filtrare ulteriormente le serie che passiamo
purrr::map(ana$percorsoSerie,.f=creaClimatolData_ii_ff,parametro=PARAM,returnMonthly=MONTHLY) %>% compact->listaOut
if(!length(listaOut)) stop("Nessuna serie idonea")
riduci(listaOut,monthlyJoin=MONTHLY)->daScrivere
# Solo se si vogliono serie mensili in output -----------------------------
#Climatol non accetta valori mensili che sono NA in tutte le serie in esame. Questi valori li dobbiamo in qualche modo riempire.
#Prima di tutto identifichiamo quei mesi all'interno di un anno in cui tutte le serie sono NA e poi li riempiamo. Questo problema non si pone
#con le serie giornaliere. Per le serie giornaliere Climatol dispone di strumenti per il refill dei dati. Lo stesso dicasi quando utilizzando
#le serie giornaliere e gli strumenti di Climatol si passa alle serie mensili. In sintesi: i mesi problematici vanno riempiti solo
#quando a Climatol si passano serie mensili calcolate senza utilizzare gli strumenti del pacchetto.
if(MONTHLY && REFILL_NA_MONTHLY) riempiNA(daScrivere)->daScrivere
#I nomi delle stazioni debbono essere: regione.codice, la parte che segue va aggiustata di volta in volta, in base a come sono stati letti i file e come
#compaiono i loro nomi nella variabile percorsoSerie. In questo caso le colonne hanno forma: "spatial_controls_piemonte_febbraio2017_serieFino_30"con il terzo elemento e il sesto
#che rappresentano rispettivamente la regione e il codice della stazione (SiteID)
purrr::map(names(daScrivere)[4:ncol(daScrivere)],.f=~(unlist(str_split(.,"_")))) %>% purrr::map_int(.,length)->lunghezza
stopifnot(all(lunghezza==6)) #sei elementi nel nome
purrr::map(names(daScrivere)[4:ncol(daScrivere)],.f=~(unlist(str_split(.,"_")))) %>% purrr::map_chr(.,3)->NOMI_REGIONI
purrr::map(names(daScrivere)[4:ncol(daScrivere)],.f=~(unlist(str_split(.,"_")))) %>% purrr::map_chr(.,6)->CODICI_REGIONI
names(daScrivere)[4:ncol(daScrivere)]<-paste0(NOMI_REGIONI,"_",CODICI_REGIONI)
#fine parte per fissare i nomi
# Scrittura file dei dati -------------------------------------------------
nome_dat<-paste0(Hmisc::capitalize(PARAM),"_",ANNOI,"-",ANNOF)
scriviFile_dat(x=daScrivere,nomeOut=nome_dat,file_monthly=MONTHLY,acmant=FALSE)
#scrivi gli stessi dati in formato RDS, come un normale dataframe invece che nel formato richiesto da climatol
#Questi dati risulteranno utili nel momento in cui vorremo confrontare i dati omogeneizzati con quelli di partenza
saveRDS(daScrivere,file=paste0(nome_dat,"_raw.RDS"))
# Scrittura del file delle coordinate: ------------------------------------
nome_est<-nome_dat
scriviFile_est(x=daScrivere,ana=ana,nomeOut=nome_est,file_monthly=MONTHLY)
# Nel caso in cui si voglia il grafico delle serie ------------------------
nome_pdf<-nome_dat
if(GRAFICO_SERIE_ANNUALI) grafico(x=daScrivere,parametro=PARAM,nomeOut = paste0(nome_pdf,"_serieAnnuali"))