Thursday-PAAC RTC Training Course
E. Di Giuseppe, M. Pasqui, M. Bacci
November 23, 2017
- Computation of main ETCCDI precipitation-related indices
- Let us importing again one of the time series of precipitation of Chirps database (say P12):
- Recall the method for creating
climdexInputobject from vectors of data - Get available indices by name
- Index n. 27: PRCPTOT-Annual total precipitation in wet days
- Index n. 26: R99pTOT. Annual total PRCP when RR > 99p.
- Comparison of probability density functions (pdf)
- Visualizing with
{ggplot2}package - How to analyse a specific period
Computation of main ETCCDI precipitation-related indices
Let us importing again one of the time series of precipitation of Chirps database (say P12):
Chirps <- read.csv("Data/CHIRPS/RHO-P12.csv",header = TRUE,sep=",")
#### 1) set column names:
names(Chirps)## [1] "day" "month" "year" "value"names(Chirps) <- c("day","month","year", "prec")
#### 2) set missing data:
Chirps[Chirps==-9999] <- NA # use NA for "not available" data
#### 2) transform month abbreviation from string to number:
## Function::
MonthStr2MonthNum<- function(x){ # x is a vector of string (month abbr in English)
xNum=match(x, month.abb)
return(xNum)
}
Chirps[,"month"]<-MonthStr2MonthNum(Chirps[,"month"])
summary(Chirps)## day month year prec
## Min. : 1.00 Min. : 1.000 Min. :1981 Min. : 0.000
## 1st Qu.: 8.00 1st Qu.: 4.000 1st Qu.:1990 1st Qu.: 0.000
## Median :16.00 Median : 6.000 Median :1999 Median : 0.000
## Mean :15.73 Mean : 6.487 Mean :1999 Mean : 2.295
## 3rd Qu.:23.00 3rd Qu.: 9.000 3rd Qu.:2008 3rd Qu.: 1.148
## Max. :31.00 Max. :12.000 Max. :2017 Max. :63.390
## NA's :19Recall the method for creating climdexInput object from vectors of data
library(climdex.pcic)## Loading required package: PCICtPREC <- Chirps$prec
PREC.DATES <- as.PCICt(do.call(paste, Chirps[,c("year","month", "day")]),
format="%Y %m %d",cal="gregorian")## Warning in strptime(x, format, tz = tz): unknown timezone 'zone/tz/2017c.
## 1.0/zoneinfo/Europe/Rome'# Create a climdexInput object:
PREC.CLIMDEX <- climdexInput.raw(prec = PREC,
prec.dates = PREC.DATES,
base.range = c(1971,2000),
northern.hemisphere = TRUE)Get available indices by name
There is a function that returns a list of (function) names of available indices for the climdexInput object.
climdex.get.available.indices(PREC.CLIMDEX, function.names = TRUE)## [1] "climdex.rx1day" "climdex.rx5day" "climdex.sdii"
## [4] "climdex.r10mm" "climdex.r20mm" "climdex.rnnmm"
## [7] "climdex.cdd" "climdex.cwd" "climdex.r95ptot"
## [10] "climdex.r99ptot" "climdex.prcptot"Index n. 27: PRCPTOT-Annual total precipitation in wet days
Let \(RR_{ij}\) be the daily precipitation amount on day i in period j. If I represents the number of days in j (days where precipitation is at least 1mm), then \[ PRCPTOT = \sum_{i}^{I} RR_{ij} \].
climdex.prcptot(PREC.CLIMDEX)## 1981 1982 1983 1984 1985 1986 1987
## NA 755.9296 704.6531 649.5334 777.4621 769.2207 679.3705
## 1988 1989 1990 1991 1992 1993 1994
## 956.4889 770.7285 719.6728 765.6774 884.7143 777.4212 1114.1505
## 1995 1996 1997 1998 1999 2000 2001
## 712.4442 797.8416 786.0353 899.2332 975.2630 746.8264 837.7817
## 2002 2003 2004 2005 2006 2007 2008
## 641.3847 962.3959 653.7390 717.6724 952.8755 883.0599 803.1679
## 2009 2010 2011 2012 2013 2014 2015
## 845.2735 1110.5289 721.4674 977.2265 839.2103 828.8681 1124.7626
## 2016 2017
## 818.3850 NAIndex n. 26: R99pTOT. Annual total PRCP when RR > 99p.
Let \(RR_{wj}\) be the daily precipitation amount on a wet day w (\(RR \ge 1.0\)mm) in period i and let RRwn99 be the 99th percentile of precipitation on wet days in the selected period. If W represents the number of wet days in the period, then: \[R99p_{j} = \sum_{w=1}^{W} RR_{wj} \; where \; RR_{wj} \ge RR_{wn}99\] .
r99ptot.71.00<-climdex.r99ptot(PREC.CLIMDEX)
r99ptot.71.00## 1981 1982 1983 1984 1985 1986 1987 1988
## NA 0.0000 0.0000 0.0000 78.4106 0.0000 0.0000 126.9771
## 1989 1990 1991 1992 1993 1994 1995 1996
## 0.0000 0.0000 35.4896 211.4968 35.3022 159.2233 0.0000 0.0000
## 1997 1998 1999 2000 2001 2002 2003 2004
## 35.0606 0.0000 37.9793 34.7534 78.6521 0.0000 121.2835 0.0000
## 2005 2006 2007 2008 2009 2010 2011 2012
## 53.8560 73.9371 157.6478 57.0891 0.0000 51.8861 0.0000 270.2457
## 2013 2014 2015 2016 2017
## 75.1054 0.0000 169.4250 71.1378 NALet us changing the base line for the computation of the percentile. Then, we have to re-build the climdexInput object.
PREC.CLIMDEX2 <- climdexInput.raw(prec = PREC,
prec.dates = PREC.DATES,
base.range = c(1961,1990),
northern.hemisphere = TRUE)and compute the annual time series of the index:
r99ptot.61.90<-climdex.r99ptot(PREC.CLIMDEX2)
r99ptot.61.90## 1981 1982 1983 1984 1985 1986 1987 1988
## NA 0.0000 0.0000 34.4260 112.8079 34.6100 0.0000 126.9771
## 1989 1990 1991 1992 1993 1994 1995 1996
## 0.0000 0.0000 35.4896 211.4968 35.3022 159.2233 34.2049 0.0000
## 1997 1998 1999 2000 2001 2002 2003 2004
## 35.0606 34.4703 37.9793 34.7534 78.6521 0.0000 121.2835 0.0000
## 2005 2006 2007 2008 2009 2010 2011 2012
## 53.8560 107.6661 191.2892 57.0891 0.0000 85.8612 33.7278 270.2457
## 2013 2014 2015 2016 2017
## 109.3338 0.0000 169.4250 71.1378 NAComparison of probability density functions (pdf)
We use density() function that uses approx to do linear interpolation:
# Compute pdf:
r99ptot.71.00.d<-density(r99ptot.71.00,na.rm=TRUE,from=0,
to=max(r99ptot.71.00,na.rm=TRUE))
plot(r99ptot.71.00.d)
r99ptot.61.90.d<-density(r99ptot.61.90,na.rm=TRUE,from=0,
to=max(r99ptot.61.90,na.rm=TRUE))
plot(r99ptot.61.90.d)
Visualizing with {ggplot2} package
library(ggplot2)
# build a dataframe:
n=length(r99ptot.71.00)
dat<-data.frame(
dens=c(r99ptot.71.00,r99ptot.61.90), lines=rep(c("r99ptot.71.00","r99ptot.61.90"),each=n)
)
ggplot(dat, aes(x = dens, fill = lines)) + geom_density(alpha = 0.5)## Warning: Removed 4 rows containing non-finite values (stat_density).
How to analyse a specific period
We use the {dplyr} package to subset data and we build again the climdexInput object setting the argument max.missing.days=c(annual = 365,monthly = 31) in order to allow the computation. In fact, the climdexInput object built from subsetted dataset is composed of the complete time series, where the values of days outside the period of interest are set to NA.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Let's say June-July-August-September:
JJAS.daily <- filter(Chirps, month==6 | month==7 | month==8 | month==9)
## summary(JJAS.daily)
# Create a climdexInput object:
PREC <- JJAS.daily$prec
PREC.DATES <- as.PCICt(do.call(paste, JJAS.daily[,c("year","month", "day")]),
format="%Y %m %d",cal="gregorian")
PREC.CLIMDEX <- climdexInput.raw(prec = PREC,
prec.dates = PREC.DATES,
base.range = c(1971,2000),
northern.hemisphere = TRUE,
max.missing.days = c(annual = 365,
monthly = 31))
climdex.r95ptot(PREC.CLIMDEX)## 1981 1982 1983 1984 1985 1986 1987 1988
## 123.9153 85.8389 28.2332 60.8925 193.8586 96.2591 56.4576 183.1738
## 1989 1990 1991 1992 1993 1994 1995 1996
## 93.7203 0.0000 35.4896 211.4968 87.8871 300.7150 34.2049 87.3719
## 1997 1998 1999 2000 2001 2002 2003 2004
## 92.8665 122.4810 124.6032 62.6685 139.2760 27.0517 177.0089 27.2312
## 2005 2006 2007 2008 2009 2010 2011 2012
## 144.9007 166.5927 250.5355 141.9651 55.8779 268.9372 33.7278 299.4028
## 2013 2014 2015 2016 2017
## 199.4926 57.4597 369.5718 215.1158 55.5294