Ricerca del threshold (Serie giornaliere dati di precipitazione SCIA)

27 gennaio 2014

Ricerca del threshold per il metodo POT. La linea verde è il possibile threshold individuato a “occhio” sulla abse dei grafici di diagnostica e del fit del modello.

La linea rossa rappresenta il 95° percentile della distribuzione dei valori non nulli. La linea gialla rappresenta il 90° percentile della distribuzione dei valori non nulli.

## ricerca del threshold per metodo POT: vengono rappresentati su mrlplot e
## tclplot i thrshold definiti a 'occhio' e quelli corrispondenti al 90 e 95
## percentile dei valori della distribuzione maggiore di 0.
rm(list = objects())
library("POT")

threshold <- function(colonna) {

    dati <- matriceDati[, c(1, 2, 3, colonna)]
    ID <- names(matriceDati)[colonna]
    names(dati)[4] <- "prcp"  ##per semplificarci la vita, per evitare dopo di usare parse(text(etc))

    # if(ID!='6721') return()

    index <- which(as.character(annoI$idstaz) == ID)
    if (!length(index)) 
        stop(paste0("--annoInizio-- Codice idstaz non trovato: ", ID))

    YYYY <- annoI[index, ]$annoI

    ## dobbiamo eliminare gli anni precedenti a YYYY e calcolare la percentuale
    ## di dati mancanti
    if (YYYY > 1951) {
        cat(sprintf("Stazione %s -- subset dati\n", ID))
        index <- which(dati$yy < YYYY)
        dati <- dati[-index, ]
    } else {
        cat(sprintf("Stazione %s -- no subset dati, serie dal 1951\n", ID))
    }

    ## creo oggetto dati.ts per
    obs.ts <- ts(dati$prcp, start = c(YYYY), freq = 365.25)
    dati.ts <- data.frame(time = time(obs.ts), obs = obs.ts)


    ## quantile per valori >0
    q98 <- quantile(obs.ts[obs.ts > 0], probs = 0.98, na.rm = T)

    ## eliminiamo i dati NA per le funzioni mrlplot e tclplot
    ## which(is.na(dati.ts$obs))->index if(length(index))
    ## dati.ts[-index,]->dati.ts

    ## leggiamo il threshold precedentemente individuato e i quantili della
    ## distribuzione
    index.thr <- which(as.character(thr$idstaz) == ID)
    if (!length(index.thr)) 
        stop(paste0("--threshold-- Codice idstaz non trovato: ", ID))

    ## se q98 è troppo vicino al valore di myThr il grafico mrplot e tcplot sono
    ## troppo 'brevi'
    if (abs(q98 - thr[index.thr, ]$myThr) <= 10) 
        q98 <- (q98 + 20)


    ## superamenti in corrispondenza di threshold pari a 10 20 30....100
    eccessi <- lapply(seq(10, q98, length = 10), FUN = function(soglia) {
        length(which(dati.ts$obs > soglia))
    })
    ## fine lapply

    eccessi <- unlist(eccessi)

    tryCatch({
        mrlplot(data = dati.ts$obs, u.range = c(10, q98), main = paste0("STAZIONE ", 
            ID), nt = 50)
        abline(v = thr[index.thr, ]$myThr, col = "green", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p90, col = "orange", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p95, col = "red", lty = 2, lwd = 2)
        mtext(text = eccessi, side = 3, at = seq(10, q98, length = 10), line = 0.45, 
            cex = 0.7)
        axis(side = 3, at = seq(10, q98, length = 10), labels = F, tick = T)
        tcplot(data = dati.ts$obs, u.range = c(10, q98), nt = 50, ask = F, which = 1)
        abline(v = thr[index.thr, ]$myThr, col = "green", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p90, col = "orange", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p95, col = "red", lty = 2, lwd = 2)
        mtext(text = eccessi, side = 3, at = seq(10, q98, length = 10), line = 0.45, 
            cex = 0.7)
        axis(side = 3, at = seq(10, q98, length = 10), labels = F, tick = T)

        tcplot(data = dati.ts$obs, u.range = c(10, q98), nt = 50, ask = F, which = 2)
        abline(v = thr[index.thr, ]$myThr, col = "green", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p90, col = "orange", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p95, col = "red", lty = 2, lwd = 2)
        mtext(text = eccessi, side = 3, at = seq(10, q98, length = 10), line = 0.45, 
            cex = 0.7)
        axis(side = 3, at = seq(10, q98, length = 10), labels = F, tick = T)

        diplot(data = dati.ts, u.range = c(10, q98), nt = 50, main = paste0("STAZIONE ", 
            ID))
        abline(v = thr[index.thr, ]$myThr, col = "green", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p90, col = "orange", lty = 2, lwd = 2)
        abline(v = thr[index.thr, ]$p95, col = "red", lty = 2, lwd = 2)
        mtext(text = eccessi, side = 3, at = seq(10, q98, length = 10), line = 0.45, 
            cex = 0.7)
        axis(side = 3, at = seq(10, q98, length = 10), labels = F, tick = T)
        abline(h = 1, col = "red", lwd = 1, lty = 1)

    }, warning = function(e) {
        paste0("Convergenza non riuscita ", ID)
    })
    ## fine tryCatch

}
## fine funzione threshold



########################################################## Inizio programma anno di inizio per ciascuna serie: idstaz, annoI
annoI <- read.csv(file = "annoInizio_serie1951.csv", sep = ";", head = T, check.names = F, 
    stringsAsFactors = F)

## lettura del file che contiene i possibili thresholds
thr <- read.csv(file = "possibleThr.csv", sep = ";", head = T, check.names = F, 
    stringsAsFactors = F)

## lista file di input e codici ID d ciascuna stazione
matriceDati <- read.csv(file = "matriceDati.csv", sep = ";", head = T, check.names = F, 
    stringsAsFactors = F)
nomiMatrice <- names(matriceDati)
numeroColonne <- ncol(matriceDati)

ldati <- lapply(4:numeroColonne, FUN = threshold)

## Stazione 6220 -- no subset dati, serie dal 1951

## Stazione 6225 -- no subset dati, serie dal 1951

## Stazione 6259 -- no subset dati, serie dal 1951

## Stazione 6278 -- no subset dati, serie dal 1951

## Stazione 6429 -- no subset dati, serie dal 1951

## Stazione 6709 -- subset dati

## Stazione 6720 -- no subset dati, serie dal 1951

## Stazione 6721 -- subset dati

## Stazione 6723 -- no subset dati, serie dal 1951

## Stazione 6729 -- no subset dati, serie dal 1951

## Stazione 6736 -- no subset dati, serie dal 1951

## Stazione 6769 -- no subset dati, serie dal 1951

## Stazione 6773 -- no subset dati, serie dal 1951

## Stazione 6779 -- subset dati

## Stazione 6782 -- subset dati

## Stazione 6794 -- subset dati

## Stazione 6812 -- subset dati

## Stazione 6817 -- no subset dati, serie dal 1951

## Stazione 6820 -- subset dati

## Stazione 6822 -- subset dati

## Stazione 6826 -- subset dati

## Stazione 6827 -- subset dati

## Stazione 6841 -- no subset dati, serie dal 1951

## Stazione 6850 -- no subset dati, serie dal 1951

## Stazione 6851 -- subset dati

## Stazione 6856 -- no subset dati, serie dal 1951

## Stazione 6869 -- no subset dati, serie dal 1951

## Stazione 6874 -- no subset dati, serie dal 1951

## Stazione 6881 -- subset dati

## Stazione 6883 -- subset dati

## Stazione 6891 -- subset dati

## Stazione 6897 -- no subset dati, serie dal 1951

## Stazione 6901 -- no subset dati, serie dal 1951

## Stazione 6923 -- no subset dati, serie dal 1951