Introduction
The primary purpose of
this study is to demonstrate the scientific feasibility of a climate change
project. To do this, it is imminent to call on a climatic database of local or
satellite stations. We are going to show how to use these programs of RClimDex,
RClimTool, RHtest and ClimPact under R. It is noteworthy to recall that
knowledge of latitude, longitude, types of variables and their faults (missing,
extreme, erroneous) are necessary for to be able to ensure the required
analyzes on climate change.
I.
RCLIMDEX
The import is done in R with an import program
allowing to display a window of RCLIMDEX. With the outputs of the RCLIMDEX
program under R, four documents including the climatic indices; the qualities
of controls and graphics and trends will be provided to us. The indices are 27
indicators divided into precipitation, maximum temperature and minimum
temperature, among which:
Temperature indices including Min and Max temperature
(SU25, TXx, TNx, TMaxMean, TNn, TMinMean, DTR,,… etc.).
Precipitation indices including humidity and rain
(PRCPTOT SDII, CDD, CWD, RX1day, RX5day, R10, R20, R95p, R99p,… etc.).
The log document provides graphs of temperature and
precipitation variations, as well as other graphic illustrations in Excel
format for tables and pdf for images.
The plot file groups together the figures of the
climate change indices associated with their respective slopes and subjected to
the Mann Kendall trend test in order to evaluate the significant trends (sign
of the slope; Student significance; regression line; correlation; coefficient
of determination). The last Trend document sufficiently summarizes the
information in graphic format from the previous file.
I.
RClimTool
L’importation The import is done under R with an
import program allowing to display a window of RClimtool. With the outputs of
the RClimTool program, three files, including homogeneity; Missing_data and
Quality_Control will have to be extracted to use. The first file that emerges
is the Quality_Control, for each of these variables, we study the atypical,
outlier or outliers values presented by our daily, monthly and annual
observations. The second is that of missing value, the processing will be done
separately with values and without missing value. The last homogeneity
document presents the various tests of normality, Mann Kendall significance,
correlation between sites, with the associated graphs. It is possible to deal
with missing values, because this prevents the calculation loop from
functioning normally. Various deterministic methods such as statistics such as
the average; the median or even probabilistic techniques (Amelia, MICE, PPCA,
etc.) are used to supplement missing values if they do not exceed certain
thresholds of significance.
II.
RCLIMTREND
Plusieurs
Méthodes statistiques pour les sciences du climat comme le Tests d'homogénéité
absolue SNHT absolu ; SD différent,
breaks, Buishand, Pettitt, rapport de von Neumann et ratio-rank, Worsley et
Craddock, tests d'homogénéité relative SNHT absolu ; break SD ; la
correlation ; la recherche des valeurs aberrantes sont les diverses atout
que l’on peut bénéficier avec ce programme.
III.
ClimPact
Le
programme de ClimPact s’exécute de la même façon que RClimDex, les issues sont
similaires cette fois-ci, une étude plus approfondie que le précèdent avec une
évaluation de variation extrême. D’autant plus que nous avons comme suppléments
sur le treshold (seuil). Les différents de statistiques descriptives (boxplot, histogramme.etc.).
En comparant les valeurs des tests de Mann Kendall respectifs sur les indices
de calculs retenus, on constate que c’est presque similaire.
IV.
RHtest
With this program, it is possible to detect point
changes of the studied series and to inform us about what needs adjustment. Geoclim ; Anclim ; EdGCM ; MASH ;
ProClimDB…etc
Spatial analysis tools are designed for analysis in climatology.
Provides scientists and non-scientists with more decision support and
facilitates filling observation gaps (precipitation and temp) using a mix of
field observations, readily available satellite data and networks.
Quality Control
The
preliminary study on data quality detects the extreme values (tools) by
realizing the statistics of maximum and minimum. The existence of missing data
disrupts the analysis of the observation series for calculating indices under
the RclimDex program (Zhang and Yang, 2004). A month containing more than 3
missing days is not considered in the year. An incomplete year is also private
in index calculations. Statistical indicators such as the median make it
possible to complete the missing values when the proportion of these values
is negligible.
Homogenization
Several
tests make it possible to detect the homogenization of the series observed. The
bivariate case like the Craddock test (1979) or the univariate case like the
PETTITT test (1979). The programs of RclimDex, Rclimtool, Rclimpact (Zhang and
Yang, 2004). allow the analysis of the homogenization of temperature and
Precipitation series.
1. Climate Change
indices
The
variability and trend of precipitation and temperature is determined with the
calculation of indices under the application of RclimDex. There are 11 for the
Precipitation and 16 for the temperature. The definition of each of these
indices and their units are presented in the following table.
Table1. List of 11 ETCCDMI core Climate Indices for
precipitation
|
Indicator |
Name |
Définitions |
|
RX1day
(mm) |
Max
1-day precipitation amount |
Monthly maximum 1-day precipitation |
|
Rx5day (mm) |
Max 5-day precipitation amount |
Monthly maximum
consecutive 5-day precipitation |
|
SDII
(mm/day) |
Simple
daily intensity index |
Annual total precipitation divided by the number of
wet days (defined as PRCP>=1.0mm) in the year |
|
R10 (day) |
Number of heavy
precipitation days |
Annual count of
days when PRCP>=10mm |
|
R20
(day) |
Number of very heavy precipitation days |
Annual count of days when PRCP>=20mm |
|
Rnn (day) |
Number of days
above nn mm |
Annual count of
days when PRCP>=nn mm, nn is user defined threshold |
|
CDD
(day) |
Consecutive
dry days |
Maximum number of consecutive days with RR<1mm |
|
CWD (day) |
Consecutive wet days |
Maximum number
of consecutive days with RR>=1mm |
|
R95p
(mm) |
Very
wet days |
Annual total PRCP when RR>95th
percentile |
|
R99p (mm) |
Extremely wet days |
Annual total
PRCP when RR>99th percentile |
|
PRCPTOT
(mm) |
Annual total wet-day precipitation |
Annual total PRCP in wet days (RR>=1mm) |
R
Program Analysis Climat Change
#TREND ANANLYSIS
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
data
str(data)
#install.packages("trend")
require(trend)
#data(maxau)
prectot <- data[,"prectot"]
mk.test(prectot)
dat=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
data(nottem)
str(dat)
dat=round(dat,1)
dat1 <- ts(dat, frequency=12,
start=c(1961,1))
plot.ts(dat1)
logsouvenirtimeseries <- log(dat1)
plot.ts(logsouvenirtimeseries)
require(trend)
smk.test(dat1)
#correlated witn precedent month
csmk.test(dat)
#Multivariate Mann-Kendall
Test
require(trend)
data(hcb)
> plot(hcb)
s <- maxau[,"s"]
Q <- maxau[,"Q"]
cor.test(s,Q, meth="spearman")
partial.mk.test(s,Q)
partial.cor.trend.test(s,Q,
"spearman")
#Cox and Stuart Trend test
cs.test(prectot)
s <- maxau[,"s"]
sens.slope(prectot)
dat=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
require(trend)
das=ts(dat, frequency=12,start=1961)
sea.sens.slope(das)
#
Changement Dectection
#data(PagesData)
require(trend)
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
attach(data)
#Pettitt's test for single change-point
detection
pettitt.test(dat)
#Buishand range test
require(trend)
(res <- br.test(prectot))
require(trend)
(res <- snh.test(Nile))
#Wallis and Moore Phase-Frequency test
frost <- ts(prectot, start=1961)
wm.test(frost)
#Bartels's test for randomness
bartels.test(prectot)
#Wald-Wolfowitz test for independence
and stationarity
ww.test(prectot)
#Wallis and Moore Phase-Frequency test
wm.test(prectot)
#Buishand U test
#Nile
(res <- bu.test(prectot))
precto=ts(prectot, frequency=12,
start=12)
par(mfrow=c(1,2))
plot(precto); plot(res)
#Standard Normal Homogeinity Test
require(trend)
#Nile
(res <- snh.test(prectot))
cor.test( x,y,alternative =
c("two.sided", "less", "greater"),
use="complete.obs",
method="kendall",conf.level = 0.95)
require(graphics)
x <- rnorm(50)
y <- runif(30)
# Do x and y come from the same
distribution?
ks.test(x, y)
# Does x come from a shifted gamma
distribution with shape 3 and rate 2?
ks.test(x+2, "pgamma", 3, 2) #
two-sided, exact
ks.test(x+2, "pgamma", 3, 2,
exact = FALSE)
ks.test(x+2, "pgamma", 3, 2,
alternative = "gr")
#
The approach after Pettitt (1979)
#
is commonly applied to detect a single change-point
#
in hydrological series or climate series with continuous data.
#
It tests the H0: The T variables follow one or more
#
distributions that have the same location parameter (no change),
#
against the alternative: a change point exists.
The Pettitt-test is conducted in such a
way
# Change-point Detection: Pettitt’s Test
library(package = dplyr, warn.conflicts
= FALSE )
library(package = ggplot2,
warn.conflicts = FALSE)
library(package = trend, warn.conflicts
= FALSE)
# Load Data
data <- read.csv(file =
'c:/Users/pooya/Downloads/Torbat Heydariyeh - Daily.csv',
header = TRUE)
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
# Compactly Display the Structure of
data
str(object = data)
# Return the First 10 Row of data
head(x = data, n = 10)
# Select Variable
Tave <- data %>%
mutate(date = as.Date(x = paste0(Year, '-', Month, '-', Day),
'%Y-%m-%d')) %>%
select(date, t)
# Compactly Display the Structure of
Tave
str(object = Tave)
# Return the First 10 Row of Tave
head(x = Tave, n = 10)
# Tave Summaries
summary(object = Tave)
# NA Remove
Tave <- na.omit(Tave)
# Plot
require(ggplot2)
ggplot(data = Tave, mapping = aes(x =
date, y = prectot)) +
geom_line()
# Pettitt’s test
pettittTest <- trend::pettitt.test(x
= Tave[['prectot']])
print(pettittTest)
print(Tave[['date']][pettittTest$estimate])
# Plot
ggplot(data = Tave, mapping = aes(x =
date, y = t)) +
geom_line() +
geom_vline(mapping = aes(xintercept =
as.numeric(Tave[['date']][pettittTest$estimate])),
linetype = 2,
colour = "red",
size = 2)
#
Break Point
#first generate some data that has
monthly observations on 6 years of a seasonal
#process followed by 2 years of monthly
that is not seasonal
#-------------------------------------------------------------------------------------
ry=read.table(file("clipboard"),header=T,
sep="\t", row.names=1)
attach(ry)
str(ry)
prectot=ts(prectot, start=1961,end=2017)
plot(prectot, ylab="Precipitation
(mm)", type="o",xlab="Year")
abline(h= mean(prectot),col='blue')
require(strucchange)
prectot=ts(prectot, start=1961,end=2017)
ocus.nile <- efp(prectot ~ 1, type =
"OLS-CUSUM")
omus.nile <- efp(prectot ~ 1, type =
"OLS-MOSUM")
rocus.nile <- efp(prectot ~ 1, type =
"Rec-CUSUM")
opar <- par(mfrow=c(2,2))
plot(prectot, ylab="Precipitation
(mm)", type="o",xlab="Year", main ="Annual
Precipitation")
abline(h= mean(prectot),col='blue')
plot(ocus.nile,
type="o",xlab="Year")
plot(omus.nile,
type="o",xlab="Year")
plot(rocus.nile,
type="o",xlab="Year")
require(strucchange)
attach(rm)
#first we are going to use the supF test
to determine whether the parameters of #our monthly dummy variable regression
are stable through the time series:
fs.days =
Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)
plot(fs.days)
#use the breakpoints functions to find
BIC corresponding to optimal number of #breaks
bp.days =
breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)
summary(bp.days)
x11()
plot(bp.days)
prectot=ts(prectot, start=1961,
end=2017)
bp.nile <- breakpoints(prectot ~ 1)
nile.fac <- breakfactor(bp.nile,
breaks = 1 )
fm1.nile <- lm(prectot~ nile.fac - 1)
plot(bp.nile)
prectot=ts(prectot, start=1961,
end=2017)
opar <- par(mfrow=c(2,1),
mar=c(2,2,0,2))
plot(ocus.nile,
alt.boundary=F,main="")
abline(v= 1980, lty=2, col='red')
plot(prectot, ylab="")
abline(h= mean(prectot),col='blue')
abline(v= 1980, lty=2, col='red')
lines(ts(predict(fm1.nile),start=1980,freq=1),
col='darkgreen',lwd=2)
par(opar)
plot(prectot)
## F statistics indicate one breakpoint
fs.nile <- Fstats(prectot ~ 1)
plot(fs.nile)
breakpoints(fs.nile)
lines(breakpoints(fs.nile))
## or
bp.nile <- breakpoints(prectot ~ 1)
summary(bp.nile)
## the BIC also chooses one breakpoint
plot(bp.nile)
breakpoints(bp.nile)
## confidence interval
ci.nile <- confint(bp.nile)
ci.nile
lines(ci.nile)
## fit null hypothesis model and model
with 1 breakpoint
fm0 <- lm(prectot ~ 1)
fm1 <- lm(prectot ~
breakfactor(bp.nile, breaks = 1))
plot(prectot, xlab="Year",
ylab="Preciptation (mm)")
lines(ts(fitted(fm0), start = 2006), col
= 3)
lines(ts(fitted(fm1), start = 2006), col
= 4)
lines(bp.nile)
abline(v=2006, lty=2, col=2)
## confidence intervals
ci.seat2 <- confint(bp.nile, breaks =
2)
ci.seat2
lines(ci.seat2)
#rm(list=ls())
data=read.csv(file("clipboard"),header=T,
sep="\t")
str(data)
if(!require(mblm)){install.packages("mblm")}
if(!require(ggplot2)){install.packages("ggplot2")}
data$Year = round(data$Year)
data$Year = data$Year - 1961
library(mblm)
model= mblm(prectot ~ Year, data=data)
summary(model)
Sum = summary(model)$coefficients
library(ggplot2)
ggplot(data, aes(x=Year, y=prectot)) +
geom_point() +
geom_abline(intercept = Sum[1], slope = Sum[2], color="blue",
size=1.2) +
labs(x = "Years after 1961")
### Optional quantile regression follows
if(!require(quantreg)){install.packages("quantreg")}
library(quantreg)
model.q = rq(prectot ~ Year, data =
data, tau = 0.5)
summary(model.q)
library(ggplot2)
model.null = rq(prectot ~ 1, data =
data, tau = 0.5)
anova(model.q, model.null)
Sumq = summary(model.q)$coefficients
ggplot(data, aes(x=Year, y=prectot)) +
geom_point() +
geom_abline(intercept = Sumq[1], slope = Sumq[2], color="red",
size=1.2) +
labs(x = "Years after 1961")
1)Seasonality
###############
r=read.table(file("clipboard"),header=T,
sep="\t", dec=".")
attach(r)
reg=lm(JFM~year)
summary(reg)
JFM=ts(JFM, start=1953, end=2017,
freq=1)
plot(JFM, type="b", lty=3)
y = -0.25x + 519.54
p-value = .09054
rm=read.table(file("clipboard"),header=T,
sep="\t", row.names=1)
getwd()
setwd("C:/Users/Lenovo/Documents")
write.table(rm,
file="anclim.txt",row.names=F)
rm=read.table(file("clipboard"),header=T,
sep="\t")
str(rm)
attach(rm)
#plot(0,0,xlim=c(1961,2017),
ylim=c(0,300))
DJF=ts(DJF, start=1961,end=2017)
Y=ts(Y, start=1961, end=2017)
#par(new=TRUE)
#plot(DJF, type="o",
xlab="", ylab="",axes=FALSE)
#par(new=TRUE)
#plot(Y, type="l",
xlab="", ylab="", add=TRUE, axes=FALSE)
names(r)
reg=lm(SNO.tav~X,data=r)
summary(reg)
plot(X,MAM, type="o",lty=9,
xlab="Year", ylab="Precipitation(mm)")
d=abline(reg)
legend("topright",c("precipitation(mm)","adjusted"),col=c("black","black"),cex=c(0.7,0.7),lty=c(9,1))
legend(locator(1),c("y= -0.323
x+677
P=.526
R²=.0073"),cex=0.65,bg="NA",box.col="white")
str(r)
library(climatol) # load the functions
of the package
homogen("r", 1980, 2017,
expl=TRUE)
#Missing Value
library(naniar)
vis_miss(r[,-1])
#dd2m("r", 1981, 2000,
homog=TRUE)
#dahstat(’Ttest’, 1981, 2000,
stat=’series’)
library(naniar)
require("UpSetR")
gg_miss_upset(r)
library(mice)
md.pattern(r)
library(reshape2)
library(ggplot2)
#Homogeneous
############
#Homogeneous
boxplot(data)
hist(data)
#median ajustement
pairs(datatable, col="blue",
main="Scatterplots")
Y=cbind(LRY)
X=cbind(LRV, LRC, INT)
#
hist(Y, prob=TRUE, col =
"blue", border = "black")
lines(density(Y))
#
OLSreg=lm(Y~X)
summary(OLSreg)
#
Qreg25=rq(Y~X, tau=0.25)
summary(Qreg25)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
anova(Qreg25, Qreg75)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
plot(sumQR)
# Check
for Break point
library("bfast")
library("zoo")
library("strucchange")
plot(Nile, ylab="Annual Flow of the river
Nile")
abline(h= mean(Nile),col='blue')
plot(merge(
Nile = as.zoo(Nile),
zoo(mean(Nile), time(Nile)),
CUSUM = cumsum(Nile - mean(Nile)),
zoo(0, time(Nile)),
MOSUM = rollapply(Nile - mean(Nile), 15, sum),
zoo(0, time(Nile))
), screen = c(1, 1, 2, 2, 3, 3), main =
"", xlab = "Time",
col = c(1, 4, 1, 4, 1, 4)
)
plot(merge(
Nile = as.zoo(Nile),
zoo(c(NA, cumsum(head(Nile, -1))/1:99),
time(Nile)),
CUSUM = cumsum(c(0, recresid(lm(Nile ~
1)))),
zoo(0, time(Nile))
), screen = c(1, 1, 2, 2), main =
"", xlab = "Time",
col = c(1, 4, 1, 4)
)
ocus.nile <- efp(Nile ~ 1, type =
"OLS-CUSUM")
omus.nile <- efp(Nile ~ 1, type =
"OLS-MOSUM")
rocus.nile <- efp(Nile ~ 1, type =
"Rec-CUSUM")
opar <- par(mfrow=c(2,2))
plot(ocus.nile)
plot(omus.nile)
plot(rocus.nile)
par(opar)
plot(1870 + 5:95, sapply(5:95,
function(i) {
before <- 1:i
after <- (i+1):100
res <- c(Nile[before] -
mean(Nile[before]), Nile[after] - mean(Nile[after]))
sum(res^2)
}), type = "b", xlab =
"Time", ylab = "RSS")
bp.nile <- breakpoints(Nile ~ 1)
nile.fac <- breakfactor(bp.nile,
breaks = 1 )
fm1.nile <- lm(Nile ~ nile.fac - 1)
plot(bp.nile)
opar <- par(mfrow=c(2,1),
mar=c(2,2,0,2))
plot(ocus.nile,
alt.boundary=F,main="")
abline(v= 1898, lty=2, col='red')
plot(Nile, ylab="Annual Flow of the
river Nile")
abline(h= mean(Nile),col='blue')
abline(v= 1898, lty=2, col='red')
lines(ts(predict(fm1.nile),start=1871,freq=1),
col='darkgreen',lwd=2)
par(opar)
#TREND ANANLYSIS
#################
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
data
str(data)
#install.packages("trend")
require(trend)
#data(maxau)
prectot <- data[,"prectot"]
mk.test(prectot)
dat=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
data(nottem)
str(dat)
dat=round(dat,1)
dat1 <- ts(dat, frequency=12,
start=c(1961,1))
plot.ts(dat1)
logsouvenirtimeseries <- log(dat1)
plot.ts(logsouvenirtimeseries)
require(trend)
smk.test(dat1)
#correlated witn precedent month
csmk.test(dat)
#Multivariate Mann-Kendall Test
require(trend)
data(hcb)
> plot(hcb)
s <- maxau[,"s"]
Q <- maxau[,"Q"]
cor.test(s,Q, meth="spearman")
partial.mk.test(s,Q)
partial.cor.trend.test(s,Q,
"spearman")
#Cox and Stuart Trend test
cs.test(prectot)
s <- maxau[,"s"]
sens.slope(prectot)
dat=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
require(trend)
das=ts(dat, frequency=12,start=1961)
sea.sens.slope(das)
#
Changement Dectection
#######################
#data(PagesData)
require(trend)
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
attach(data)
#Pettitt's test for single change-point
detection
pettitt.test(dat)
#Buishand range test
require(trend)
(res <- br.test(prectot))
require(trend)
(res <- snh.test(Nile))
#Wallis and Moore Phase-Frequency test
frost <- ts(prectot, start=1961)
wm.test(frost)
#Bartels's test for randomness
bartels.test(prectot)
#Wald-Wolfowitz test for independence
and stationarity
ww.test(prectot)
#Wallis and Moore Phase-Frequency test
wm.test(prectot)
#Buishand U test
#Nile
(res <- bu.test(prectot))
precto=ts(prectot, frequency=12,
start=12)
par(mfrow=c(1,2))
plot(precto); plot(res)
#Standard Normal Homogeinity Test
require(trend)
#Nile
(res <- snh.test(prectot))
cor.test( x,y,alternative =
c("two.sided", "less", "greater"),
use="complete.obs",
method="kendall",conf.level = 0.95)
require(graphics)
x <- rnorm(50)
y <- runif(30)
# Do x and y come from the same
distribution?
ks.test(x, y)
# Does x come from a shifted gamma
distribution with shape 3 and rate 2?
ks.test(x+2, "pgamma", 3, 2) #
two-sided, exact
ks.test(x+2, "pgamma", 3, 2,
exact = FALSE)
ks.test(x+2, "pgamma", 3, 2,
alternative = "gr")
#
The approach after Pettitt (1979)
#
is commonly applied to detect a single change-point
#
in hydrological series or climate series with continuous data.
#
It tests the H0: The T variables follow one or more
#
distributions that have the same location parameter (no change),
#
against the alternative: a change point exists.
The Pettitt-test is conducted in such a
way
# Change-point Detection: Pettitt’s Test
library(package = dplyr, warn.conflicts
= FALSE )
library(package = ggplot2,
warn.conflicts = FALSE)
library(package = trend, warn.conflicts
= FALSE)
# Load Data
data <- read.csv(file =
'c:/Users/pooya/Downloads/Torbat Heydariyeh - Daily.csv',
header = TRUE)
data=read.csv(file("clipboard"),header=T,
sep="\t", row.names=1)
# Compactly Display the Structure of
data
str(object = data)
# Return the First 10 Row of data
head(x = data, n = 10)
# Select Variable
Tave <- data %>%
mutate(date = as.Date(x = paste0(Year, '-', Month, '-', Day),
'%Y-%m-%d')) %>%
select(date, t)
# Compactly Display the Structure of
Tave
str(object = Tave)
# Return the First 10 Row of Tave
head(x = Tave, n = 10)
# Tave Summaries
summary(object = Tave)
# NA Remove
Tave <- na.omit(Tave)
# Plot
require(ggplot2)
ggplot(data = Tave, mapping = aes(x =
date, y = prectot)) +
geom_line()
# Pettitt’s test
pettittTest <- trend::pettitt.test(x
= Tave[['prectot']])
print(pettittTest)
print(Tave[['date']][pettittTest$estimate])
# Plot
ggplot(data = Tave, mapping = aes(x =
date, y = t)) +
geom_line() +
geom_vline(mapping = aes(xintercept = as.numeric(Tave[['date']][pettittTest$estimate])),
linetype = 2,
colour = "red",
size = 2)
#
Break Point
#first generate some data that has
monthly observations on 6 years of a seasonal
#process followed by 2 years of monthly
that is not seasonal
#-------------------------------------------------------------------------------------
ry=read.table(file("clipboard"),header=T,
sep="\t", row.names=1)
attach(ry)
str(ry)
prectot=ts(prectot, start=1961,end=2017)
plot(prectot, ylab="Precipitation
(mm)", type="o",xlab="Year")
abline(h= mean(prectot),col='blue')
require(strucchange)
prectot=ts(prectot, start=1961,end=2017)
ocus.nile <- efp(prectot ~ 1, type =
"OLS-CUSUM")
omus.nile <- efp(prectot ~ 1, type =
"OLS-MOSUM")
rocus.nile <- efp(prectot ~ 1, type =
"Rec-CUSUM")
opar <- par(mfrow=c(2,2))
plot(prectot, ylab="Precipitation
(mm)", type="o",xlab="Year", main ="Annual
Precipitation")
abline(h= mean(prectot),col='blue')
plot(ocus.nile,
type="o",xlab="Year")
plot(omus.nile,
type="o",xlab="Year")
plot(rocus.nile,
type="o",xlab="Year")
require(strucchange)
attach(rm)
#first we are going to use the supF test
to determine whether the parameters of #our monthly dummy variable regression
are stable through the time series:
fs.days =
Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)
plot(fs.days)
#use the breakpoints functions to find
BIC corresponding to optimal number of #breaks
bp.days =
breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)
summary(bp.days)
x11()
plot(bp.days)
prectot=ts(prectot, start=1961,
end=2017)
bp.nile <- breakpoints(prectot ~ 1)
nile.fac <- breakfactor(bp.nile,
breaks = 1 )
fm1.nile <- lm(prectot~ nile.fac - 1)
plot(bp.nile)
prectot=ts(prectot, start=1961,
end=2017)
opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))
plot(ocus.nile,
alt.boundary=F,main="")
abline(v= 1980, lty=2, col='red')
plot(prectot, ylab="")
abline(h= mean(prectot),col='blue')
abline(v= 1980, lty=2, col='red')
lines(ts(predict(fm1.nile),start=1980,freq=1),
col='darkgreen',lwd=2)
par(opar)
plot(prectot)
## F statistics indicate one breakpoint
fs.nile <- Fstats(prectot ~ 1)
plot(fs.nile)
breakpoints(fs.nile)
lines(breakpoints(fs.nile))
## or
bp.nile <- breakpoints(prectot ~ 1)
summary(bp.nile)
## the BIC also chooses one breakpoint
plot(bp.nile)
breakpoints(bp.nile)
## confidence interval
ci.nile <- confint(bp.nile)
ci.nile
lines(ci.nile)
## fit null hypothesis model and model
with 1 breakpoint
fm0 <- lm(prectot ~ 1)
fm1 <- lm(prectot ~
breakfactor(bp.nile, breaks = 1))
plot(prectot, xlab="Year",
ylab="Preciptation (mm)")
lines(ts(fitted(fm0), start = 2006), col
= 3)
lines(ts(fitted(fm1), start = 2006), col
= 4)
lines(bp.nile)
abline(v=2006, lty=2, col=2)
## confidence intervals
ci.seat2 <- confint(bp.nile, breaks =
2)
ci.seat2
lines(ci.seat2)
#rm(list=ls())
data=read.csv(file("clipboard"),header=T,
sep="\t")
str(data)
if(!require(mblm)){install.packages("mblm")}
if(!require(ggplot2)){install.packages("ggplot2")}
data$Year = round(data$Year)
data$Year = data$Year - 1961
library(mblm)
model= mblm(prectot ~ Year, data=data)
summary(model)
Sum = summary(model)$coefficients
library(ggplot2)
ggplot(data, aes(x=Year, y=prectot)) +
geom_point() +
geom_abline(intercept = Sum[1], slope = Sum[2], color="blue",
size=1.2) +
labs(x = "Years after 1961")
### Optional quantile regression follows
if(!require(quantreg)){install.packages("quantreg")}
library(quantreg)
model.q = rq(prectot ~ Year, data =
data, tau = 0.5)
summary(model.q)
library(ggplot2)
model.null = rq(prectot ~ 1, data =
data, tau = 0.5)
anova(model.q, model.null)
Sumq = summary(model.q)$coefficients
ggplot(data, aes(x=Year, y=prectot)) +
geom_point() +
geom_abline(intercept = Sumq[1], slope = Sumq[2], color="red",
size=1.2) +
labs(x = "Years after 1961")
#Calculus for Index CLIMATE CHANGE WITH
R
#
0> >> <<
Importation Data
rm(list=ls())
djibclim=read.table(file("clipboard"),header=T,
dec=".",sep="\t")
#File Change dir
#getwd()
#setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/A.ClimDjibouti/0.Data")
#save(djiclim, file =
"djiclim.rda")
#load(file = "djibclim.rda")
#View(djibclim)
#str(djibclim);summary(djiclim)
#sapply(djibclim,sd)
#sapply(djibclim, mean, na.rm=TRUE)
#library(psych)
#describe(djiclim[,4:6])
#Different Plot Types:
#plot(r,type="l",col=2)
#plot(r,type="b")
#plot(r,type="p")
#plot(r,type="o")
#plot(r,type="h")
#plot(r,type="s", col=4)
#plot(r,type="c")
plot(r,type="o",
main="", xlab="Year", ylab="Precipitation (mm)")
attach(r)
summary(lm(prectot~Year, data=r))
#
#
1> >> <<
source import Rclimtool
#File change dir
#getwd()
setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project
2")
source("RClimTool\\RClimTool.R")
#source("RClimTool\\RClimTool_v2.R")
#folder na gen, na, missing before
homogeneous
#day month year without label
#
#tmax = cargar(gfile(type =
"open"))
#tmin = cargar(gfile(type =
"open"))
#precip = cargar(gfile(type =
"open"))
#
#plot(x = day, main = "Title",
xlab = "Label for the x axis", ylab = "Label for the y
axis")
#hist(x = day, main = "Title",
xlab = "Label for the x axis", ylab = "Label for the y
axis")
#boxplot(x = day, main =
"Title", ylab = "Label for the y axis")
#
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\Rclimtrends\\R")
#library(climtrends)
source("climtrends.R")
names(djibclim)
attach(djibclim)
CraddockTest(Days[,Temp_Tn])
#
# Craddock test for testing the
homogeneity of yearly average temperatures of Milan vs Turin
data(yearly.average.temperature.Turin.Milan)
craddock.result <-
CraddockTest(yearly.average.temperature.Turin.Milan,2,3)
plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')
#
# Example of inhomogeneity
tmp <-
yearly.average.temperature.Turin.Milan
tmp[20:43,3] <- tmp[20:43,3]+1 #
adding an error of +1 degree to Milan's yearly average temperatures from 1980
craddock.result <-
CraddockTest(tmp,2,3)
dev.new()
plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')
#
#Homogénéity Analysis
names(djiclim)
library(ggpubr)
ggqqplot(djiclim$Rain_Rrr)
#
library("car")
qqPlot(djiclim$Rain_Rrr)
#
#shapiro.test(rnorm(100, mean = 5, sd =
3))
shapiro.test(runif(100, min = 2, max =
4))
#require(stats)
shapiro.test(djibclim$Rain_Rrr)
#
library(tseries)
jarque.bera.test(djibclim$Rain_Rrr)
#
#Kolmogorov-Smirnov Tests
t.test(djibclim$Rain_Rrr)
wilcox.test(djibclim$Rain_Rrr)
ks.test(djibclim$Rain_Rrr)
#
#
#
#
2> >> << Program Rclim Dex
#library("RClimDex")
#rclimdex.start()
getwd()
setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project
2/Climat/RclimDex")
source("1.0.RclimDex.R.txt")
require("climdex.pcic")
require("Kendall")
require(tcltk)
#year month day without label
#I+1
i-1
#
#
3> >> << Program ClimPact2
#Download the above zip file to your
computer
#(https://github.com/ARCCSS-extremes/climpact2/archive/master.zip)
# and extract it. This will create a
directory named climpact2-master.
#In Windows: open R and select
"File -> Change dir..."
#and select the climpact2-master
directory created in step 1.
#Then type source('climpact2.GUI.r').
#(NOTE: If nothing happens, try run the
additional command startss()).
#In Linux/macOS: change to the
climpact2-master directory created in step 1,
#then open R in a terminal window and
type source('climpact2.GUI.r').
#The first time Climpact2 is run it will
install required R packages.
#This will likely require you to select
a mirror from which to download.
getwd()
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\2.Project
2\\Climat\\RClimPact\\climpact2-master")
source("climpact2.GUI.r")
#
#
#
4> >> << Program RHtest
#setwd("D:\\Docs.Abdi-Basid.ADAN\\Doc.LMME\\RHtest")
#source ("RHtest_V4.r")
#source ("RHtest.txt")
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\RHtest\\RHtests-master")
source("RHtestsV4_20180301.r")
StartGUI()
#
4'> >> << Program RHtest
library(PCICt)
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
CDD index.
cdd <- climdex.cdd(ci)
#
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
cold spell duration index.
csdi <- climdex.csdi(ci)
#
#
#
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
CWD index.
cdd <- climdex.cdd(ci)
#
#
5> >> <<
ANCLIM
at=read.table(file("clipboard"),header=T,
sep="\t")
write(at, file = "data.txt",
sep = " ")
#
5> >> <<
Geoclim
#
6> >> <<
GeoDa
# ANNEXE
########
require(ggplot2)
ggplot(data = quakes) +
geom_density(mapping = aes(x = mag, fill = region), alpha = 0.5)
library("RColorBrewer")
display.brewer.all()
# Missing Data Imputation
#########################
library(mice)
Mousey = mice(Greenwich)
Greenwich = complete(Mousey)
misp=read.csv(file("clipboard"),
header=T, sep="\t")
str(misp)
misp=misp[,-c(1:3)]
library(naniar)
vis_miss(misp)
#Prediction
###########
#Import Data
data=read.table(file("clipboard"),
header=T, sep="\t", dec=".")
str(data)
summary(data)
library(forecast)
fit.series3 <- auto.arima(prectot)
fcast.series3 <-
forecast(fit.series3)
plot(fcast.series3)
ts_passengers =
ts(prectot,start=1961,end=2017,frequency=1)
plot(ts_passengers,xlab="",
ylab="")
m_ets = ets(ts_passengers)
f_ets = forecast(m_ets, h=24) # forecast
24 months into the future
plot(f_ets, type="o",
xlab="Year", ylab="Precipitaion(mm)")
m_tbats = tbats(ts_passengers)
f_tbats = forecast(m_tbats, h=24)
plot(f_tbats)
m_aa = auto.arima(ts_passengers)
f_aa = forecast(m_aa, h=24)
plot(f_aa)
barplot(c(ETS=m_ets$aic, ARIMA=m_aa$aic,
TBATS=m_tbats$AIC),
col="light blue",
ylab="AIC")
#https://a-little-book-of-r-for-time-series.readthedocs.io/en/latest/src/timeseries.html
#https://otexts.com/fpp2/arima-r.html
prectot %>% stl(s.window='periodic')
%>% seasadj() -> prectot
autoplot(prectot)
#Steps to be followed for ARIMA
modeling:
#1. Exploratory analysis
#2. Fit the model
#3. Diagnostic measures
tsData = ts(prectot, start = c(1961,1),
end=c(2017,1),frequency = 12)
#1. Exploratory analysis
components.ts = decompose(tsData)
plot(components.ts)
library("fUnitRoots")
urkpssTest(tsData, type =
c("tau"), lags = c("short"),use.lag = NULL, doplot = TRUE)
tsstationary = diff(tsData,
differences=1)
plot(tsstationary)
acf(tsstationary, lag.max=34)
pacf(tsstationary, lag.max=34)
fitARIMA <- arima(tsData,
order=c(1,1,1),seasonal = list(order = c(1,0,0), period =
12),method="ML")
library(lmtest)
coeftest(fitARIMA)
confint(fitARIMA)
acf(fitARIMA$residuals)
library(FitAR)
LjungBoxTest (fitARIMA$residuals,k=2,StartLag=1)
qqnorm(fitARIMA$residuals)
qqline(fitARIMA$residuals)
auto.arima(tsData, trace=TRUE)
predict(fitARIMA,n.ahead = 20)
require(forecast)
futurVal <-
forecast.Arima(fitARIMA,h=10, level=c(99.5))
plot.forecast(futurVal)
plot(fitARIMA,h=10, level=c(99.5))
A=auto.arima(prectot)
autoplot(forecast(A))
autoplot(forecast(prectot), n=100)
library(forecast)
AutoArimaModel=auto.arima(prectot)
AutoArimaModel
plot(AutoArimaModel)
require(tseries); require(astsa)
acf(prectot)
pacf(prectot)
precto=ts(prectot, frequency=12,
start=1961, end=2017)
fit = arima(precto, order = c(0, 0, 0))
pred = predict(fit, n.ahead = 50)
fit = arima(log(precto), c(1, 1, 0),
seasonal = list(order = c(0, 1, 1), period = 12))
pred <- predict(fit, n.ahead = 10*12)
ts.plot(precto,exp(pred$pred), log =
"y", lty = c(1,3))
od <- options(digits = 5) # avoid too
much spurious accuracy
q=predict(arima(lh, order = c(3,0,0)),
n.ahead = 12)
(fit <- arima(USAccDeaths, order =
c(0,1,1),
seasonal = list(order =
c(0,1,1))))
predict(fit, n.ahead = 6)
ts(vector, start=, end=, frequency=)
# subset the time series (June 2014 to
December 2014)
myts2 <- window(myts, start=c(2014,
6), end=c(2014, 12))
# plot series
plot(myts)
# Seasonal decomposition
fit <- stl(myts,
s.window="period")
plot(fit)
# additional plots
monthplot(myts)
library(forecast)
seasonplot(myts)
# simple exponential - models level
fit <- HoltWinters(myts, beta=FALSE,
gamma=FALSE)
# double exponential - models level and
trend
fit <- HoltWinters(myts, gamma=FALSE)
# triple exponential - models level,
trend, and seasonal components
fit <- HoltWinters(myts)
# predictive accuracy
library(forecast)
accuracy(fit)
# predict next three future values
library(forecast)
forecast(fit, 3)
plot(forecast(fit, 3))
# fit an ARIMA model of order P, D, Q
fit <- arima(myts, order=c(p, d, q)
# predictive accuracy
library(forecast)
accuracy(fit)
# predict next 5 observations
library(forecast)
forecast(fit, 5)
plot(forecast(fit, 5))
library(forecast)
# Automated forecasting using an
exponential model
fit <- ets(myts)
# Automated forecasting using an ARIMA
model
fit <- auto.arima(myts)
#The forecast package provides functions
for the automatic selection
#of exponential and ARIMA models. The
ets() function
#supports both additive and
multiplicative models.
#The auto.arima() function can handle
both seasonal
#and nonseasonal ARIMA models. Models
are chosen to maximize one of
#several fit criteria.
#lag(ts, k) lagged version of time series, shifted back k observations
#diff(ts, differences=d) difference the time series d times
#ndiffs(ts) Number of differences required to achieve stationarity (from
the forecast package)
#acf(ts) autocorrelation
function
#pacf(ts) partial autocorrelation function
#adf.test(ts) Augemented Dickey-Fuller test. Rejecting the null hypothesis
suggests that a time series is stationary (from the tseries package)
#Box.test(x, type="Ljung-Box") Pormanteau test that observations in
vector or time series x are independent
x <-
c(32,64,96,118,126,144,152.5,158)
y <-
c(99.5,104.8,108.5,100,86,64,35.3,15)
x <- 1:10
y <- x + c(-0.5,0.5)
plot(x,y, xlim=c(0,11), ylim=c(-1,12))
fit1 <- lm( y~offset(x) -1 )
fit2 <- lm( y~x )
fit3 <- lm( y~poly(x,3) )
fit4 <- lm( y~poly(x,9) )
library(splines)
fit5 <- lm( y~ns(x, 3) )
fit6 <- lm( y~ns(x, 9) )
fit7 <- lm( y ~ x + cos(x*pi) )
xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1,
data.frame(x=xx)), col='blue')
lines(xx, predict(fit2,
data.frame(x=xx)), col='green')
lines(xx, predict(fit3,
data.frame(x=xx)), col='red')
lines(xx, predict(fit4,
data.frame(x=xx)), col='purple')
lines(xx, predict(fit5,
data.frame(x=xx)), col='orange')
lines(xx, predict(fit6,
data.frame(x=xx)), col='grey')
lines(xx, predict(fit7,
data.frame(x=xx)), col='black')
polyfit <- function(i) x <-
AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval =
c(1,length(x)-1))$minimum)
for (i in 2:length(x)-1)
print(polyfit(i))
lm(y ~ x + I(x^2) + I(x^3))
lm(y ~ poly(x, 3, raw=TRUE))
#Reseau Neurone Prediction
# reseaux neurone application sous R
library(MASS)
library(nnet)
# apprentissage
nnet.reg=nnet(O3obs~.,data=datappr,size=5,decay=1,linout=TRUE,maxit=500)
summary(nnet.reg)
library(e1071)
plot(tune.nnet(O3obs~.,data=datappr,size=c(2,3,4),
decay=c(1,2,3),maxit=200,linout=TRUE))
plot(tune.nnet(O3obs~.,data=datappr,size=4:5,decay=1:10)
nnet.reg=nnet(O3obs~.,data=datappr,size=3,decay=2,
linout=TRUE,maxit=200)
# calcul et graphe des résidus
fit.nnetr=predict(nnet.reg,data=datappr)
res.nnetr=fit.nnetr-datappr[,"O3obs"]
plot.res(fit.nnetr,res.nnetr)
# apprentissage
nnet.dis=nnet(DepSeuil~.,data=datappq,size=5,decay=1)
summary(nnet.reg)
#matrice de confusion
table(nnet.dis$fitted.values>0.5,datappq$DepSeuil)
function(formula, data, size, niter = 1,
nplis = 10, decay = 0, maxit = 100)
{
n = nrow(data)
tmc=0
un = rep(1, n)
ri = sample(nplis, n, replace = T)
cat(" k= ")
for(i in sort(unique(ri))) {
cat(" ", i, sep =
"")
for(rep in 1:niter) {
learn = nnet(formula, data[ri != i, ],
size = size,trace = F, decay = decay,
maxit = maxit)
tmc = tmc + sum(un[(data$DepSeuil[ri ==
i]
== "TRUE") != (predict(learn,
data[ri == i,
]) > 0.5)])
}
}
cat("\n", "Taux de mal
classes")
tmc/(niter * length(unique(ri)) * n)
}
CVnn(DepSeuil~.,data=datappq,size=7,
decay=0)
...
# exécuter pour différentes valeur du decay
pred.nnetr=predict(nnet.reg,newdata=datestr)
pred.nnetq=predict(nnet.dis,newdata=datestq)
# Erreur quadratique moyenne de prévision
sum((pred.nnetr-datestr[,"O3obs"])^2)/nrow(datestr)
# Matrice de confusion pour la prévision du
# dépassement de seuil (régression)
table(pred.nnetr>150,datestr[,"O3obs"]>150)
# Même chose pour la discrimination
table(pred.nnetq>0.5,datestq[,"DepSeuil"])
library(ROCR)
rocnnetr=pred.nnetr/300
prednnetr=prediction(rocnnetr,datestq$DepSeuil)
perfnnetr=performance(prednnetr,"tpr","fpr")
rocnnetq=pred.nnetq
prednnetq=prediction(rocnnetq,datestq$DepSeuil)
perfnnetq=performance(prednnetq,"tpr","fpr")
# tracer les courbes ROC en les superposant
# pour mieux comparer
plot(perflogit,col=1)
plot(perfnnetr,col=2,add=TRUE)
plot(perfnnetq,col=3,add=TRUE)
#Extreme Value Analysis for Climate
Change
r=read.table(file("clipboard"),header=TRUE,sep="\t",
row.names=1)
str(r)
attach(r)
prcp=precipitation
moy=mean(precipitation)
sigma=sqrt(1/(length(precipitation)-1))*sum((precipitation-moy)^2)
k=(sigma/moy)^-1.086
x=precipitation
GEV=(1/sigma)*exp(-
(1+(k*(x-moy)/sigma))^-(1/k)
)*(1+(k*(x-moy)/sigma ))^(-1-(1/k))
GEV
barplot(GEV,type="o")
GP=(1/sigma)*(1+(k*((x-moy)/sigma))
^(-1-(1/k)))
GP
plot(GP,type="o")
GEV/GP
# load packages
library(extRemes)
library(xts)
# get data from eHYD
ehyd_url <-
"http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=105700&file=2"
precipitation_xts <-
read_ehyd(ehyd_url)
# mean residual life plot:
mrlplot(precipitation_xts,
main="Mean Residual Life Plot")
# The mean residual life plot depicts
the Thresholds (u) vs Mean Excess flow.
# The idea is to ?nd the lowest
threshold where the plot is nearly linear;
# taking into account the 95% con?dence
bounds.
# fitting the GPD model over a range of
thresholds
threshrange.plot(precipitation_xts, r =
c(30, 45), nint = 16)
# ismev implementation is faster:
# ismev::gpd.fitrange(precipitation_xts,
umin=30, umax=45, nint = 16)
# set threshold
th <- 40
# maximum likelihood estimation
pot_mle <-
fevd(as.vector(precipitation_xts), method = "MLE",
type="GP", threshold=th)
# diagnostic plots
plot(pot_mle)
rl_mle <- return.level(pot_mle, conf
= 0.05, return.period= c(2,5,10,20,50,100))
# L-moments estimation
pot_lmom <-
fevd(as.vector(precipitation_xts), method = "Lmoments",
type="GP", threshold=th)
# diagnostic plots
plot(pot_lmom)
rl_lmom <- return.level(pot_lmom,
conf = 0.05, return.period= c(2,5,10,20,50,100))
# return level plots
par(mfcol=c(1,2))
# return level plot w/ MLE
plot(pot_mle, type="rl",
main="Return Level Plot for Oberwang w/ MLE",
ylim=c(0,200), pch=16)
loc <-
as.numeric(return.level(pot_mle, conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col=
'midnightblue',lty=6)
segments(0.01,loc,100, loc,
col='midnightblue', lty=6)
# return level plot w/ LMOM
plot(pot_lmom, type="rl",
main="Return Level Plot for Oberwang w/ L-Moments",
ylim=c(0,200))
loc <- as.numeric(return.level(pot_lmom,
conf = 0.05,return.period=100))
segments(100, 0, 100, loc, col=
'midnightblue',lty=6)
segments(0.01,loc,100, loc,
col='midnightblue', lty=6)
# comparison of return levels
results <-
t(data.frame(mle=as.numeric(rl_mle),
lmom=as.numeric(rl_lmom)))
colnames(results) <-
c(2,5,10,20,50,100)
round(results,1)
library(extRemes)
library(ismev)
# Southwest-England rainfall data from
Coles
data(rain)
# simple threshold (usually one should
not use fixed quantiles)
u < - quantile(rain, 0.9)
# fit GP using Bayesian parameter
estimation
x <- fevd(rain, threshold = u ,
type='GP', method='Bayesian')
B1.fit=fevd(B1[,2], data=B1,
threshold=B1.th, type="GEV",method ="MLE",
units="KVA",time.units="seconds",period =
"Seconds")
B1.fit1=fevd(B1[,2], data=B1,
threshold=B1.th,type="GP",method ="MLE",
units="KVA",time.units="seconds",period =
"Seconds")
B1.fit2=fevd(B1[,2], data=B1,
threshold=B1.th,type="Gumbel", method
="MLE",units="KVA",time.units="seconds",period =
"Seconds")
B1.fit3=fevd(B1[,2], data=B1,
threshold=B1.th,type="Exponential",method ="MLE",
units="KVA",time.units="seconds",period =
"Seconds")
fit.AIC=summary(B1.fit, silent=TRUE)$AIC
fit1.AIC=summary(B1.fit1,
silent=TRUE)$AIC
fit2.AIC=summary(B1.fit2,
silent=TRUE)$AIC
fit3.AIC=summary(B1.fit3,
silent=TRUE)$AIC
fit.AIC
# [1] 39976258
fit1.AIC
# [1] 466351.5
fit2.AIC
# [1] 13934878
fit3.AIC
# [1] 466330.8
plot(B1.fit)
plot(B1.fit1)
plot(B1.fit2)
plot(B1.fit3)
rm(list=ls())
rm(list=ls())
rm(list=ls())
rm=read.table(file("clipboard"),header=T,
sep="\t", row.names=1)
getwd()
setwd("C:/Users/Lenovo/Documents")
write.table(rm,
file="anclim.txt",row.names=F)
rm=read.table(file("clipboard"),header=T,
sep="\t")
str(rm)
attach(rm)
#plot(0,0,xlim=c(1961,2017),
ylim=c(0,300))
DJF=ts(DJF, start=1961,end=2017)
Y=ts(Y, start=1961, end=2017)
#par(new=TRUE)
#plot(DJF, type="o",
xlab="", ylab="",axes=FALSE)
#par(new=TRUE)
#plot(Y, type="l",
xlab="", ylab="", add=TRUE, axes=FALSE)
names(r)
reg=lm(SNO.tav~X,data=r)
summary(reg)
plot(X,MAM, type="o",lty=9,
xlab="Year", ylab="Precipitation(mm)")
d=abline(reg)
legend("topright",c("precipitation(mm)","adjusted"),col=c("black","black"),cex=c(0.7,0.7),lty=c(9,1))
legend(locator(1),c("y= -0.323
x+677
P=.526
R²=.0073"),cex=0.65,bg="NA",box.col="white")
str(r)
library(climatol) # load the functions
of the package
homogen("r", 1980, 2017,
expl=TRUE)
#Missing Value
library(naniar)
vis_miss(r[,-1])
#dd2m("r", 1981, 2000,
homog=TRUE)
#dahstat(’Ttest’, 1981, 2000,
stat=’series’)
library(naniar)
require("UpSetR")
gg_miss_upset(r)
library(mice)
md.pattern(r)
library(reshape2)
library(ggplot2
Abdi-Basid Analysis
abdi-basid@outlook.com
#Import Data
data=read.table(file("clipboard"),
header=T, sep="\t", dec=".")
str(data)
summary(data)
#Homogeneous
boxplot(data)
hist(data)
#median ajustement
pairs(datatable, col="blue",
main="Scatterplots")
Y=cbind(LRY)
X=cbind(LRV, LRC, INT)
#
hist(Y, prob=TRUE, col =
"blue", border = "black")
lines(density(Y))
#
OLSreg=lm(Y~X)
summary(OLSreg)
#
Qreg25=rq(Y~X, tau=0.25)
summary(Qreg25)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
anova(Qreg25, Qreg75)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
plot(sumQR)
# Check for Break point
library("bfast")
library("zoo")
library("strucchange")
ts(vector, start=, end=, frequency=)
# subset the time series (June 2014 to
December 2014)
myts2 <- window(myts, start=c(2014,
6), end=c(2014, 12))
# plot series
plot(myts)
# Seasonal decomposition
fit <- stl(myts,
s.window="period")
plot(fit)
# additional plots
monthplot(myts)
library(forecast)
seasonplot(myts)
# simple exponential - models level
fit <- HoltWinters(myts, beta=FALSE,
gamma=FALSE)
# double exponential - models level and
trend
fit <- HoltWinters(myts, gamma=FALSE)
# triple exponential - models level,
trend, and seasonal components
fit <- HoltWinters(myts)
# predictive accuracy
library(forecast)
accuracy(fit)
# predict next three future values
library(forecast)
forecast(fit, 3)
plot(forecast(fit, 3))
# fit an ARIMA model of order P, D, Q
fit <- arima(myts, order=c(p, d, q)
# predictive accuracy
library(forecast)
accuracy(fit)
# predict next 5 observations
library(forecast)
forecast(fit, 5)
plot(forecast(fit, 5))
library(forecast)
# Automated forecasting using an
exponential model
fit <- ets(myts)
# Automated forecasting using an ARIMA
model
fit <- auto.arima(myts)
#The forecast package provides functions
for the automatic selection
#of exponential and ARIMA models. The
ets() function
#supports both additive and
multiplicative models.
#The auto.arima() function can handle
both seasonal
#and nonseasonal ARIMA models. Models
are chosen to maximize one of
#several fit criteria.
#lag(ts, k) lagged version of time series, shifted back k observations
#diff(ts, differences=d) difference the time series d times
#ndiffs(ts) Number of differences required to achieve stationarity (from
the forecast package)
#acf(ts) autocorrelation
function
#pacf(ts) partial autocorrelation function
#adf.test(ts) Augemented Dickey-Fuller test. Rejecting the null hypothesis
suggests that a time series is stationary (from the tseries package)
#Box.test(x, type="Ljung-Box") Pormanteau test that observations in
vector or time series x are independent
plot(Nile, ylab="Annual Flow of the river
Nile")
abline(h= mean(Nile),col='blue')
plot(merge(
Nile = as.zoo(Nile),
zoo(mean(Nile), time(Nile)),
CUSUM = cumsum(Nile - mean(Nile)),
zoo(0, time(Nile)),
MOSUM = rollapply(Nile - mean(Nile), 15, sum),
zoo(0, time(Nile))
), screen = c(1, 1, 2, 2, 3, 3), main =
"", xlab = "Time",
col = c(1, 4, 1, 4, 1, 4)
)
plot(merge(
Nile = as.zoo(Nile),
zoo(c(NA, cumsum(head(Nile, -1))/1:99),
time(Nile)),
CUSUM = cumsum(c(0, recresid(lm(Nile ~
1)))),
zoo(0, time(Nile))
), screen = c(1, 1, 2, 2), main =
"", xlab = "Time",
col = c(1, 4, 1, 4)
)
ocus.nile <- efp(Nile ~ 1, type =
"OLS-CUSUM")
omus.nile <- efp(Nile ~ 1, type =
"OLS-MOSUM")
rocus.nile <- efp(Nile ~ 1, type =
"Rec-CUSUM")
opar <- par(mfrow=c(2,2))
plot(ocus.nile)
plot(omus.nile)
plot(rocus.nile)
par(opar)
plot(1870 + 5:95, sapply(5:95,
function(i) {
before <- 1:i
after <- (i+1):100
res <- c(Nile[before] -
mean(Nile[before]), Nile[after] - mean(Nile[after]))
sum(res^2)
}), type = "b", xlab =
"Time", ylab = "RSS")
bp.nile <- breakpoints(Nile ~ 1)
nile.fac <- breakfactor(bp.nile,
breaks = 1 )
fm1.nile <- lm(Nile ~ nile.fac - 1)
plot(bp.nile)
opar <- par(mfrow=c(2,1),
mar=c(2,2,0,2))
plot(ocus.nile,
alt.boundary=F,main="")
abline(v= 1898, lty=2, col='red')
plot(Nile, ylab="Annual Flow of the
river Nile")
abline(h= mean(Nile),col='blue')
abline(v= 1898, lty=2, col='red')
lines(ts(predict(fm1.nile),start=1871,freq=1),
col='darkgreen',lwd=2)
par(opar)
x <-
c(32,64,96,118,126,144,152.5,158)
y <-
c(99.5,104.8,108.5,100,86,64,35.3,15)
x <- 1:10
y <- x + c(-0.5,0.5)
plot(x,y, xlim=c(0,11), ylim=c(-1,12))
fit1 <- lm( y~offset(x) -1 )
fit2 <- lm( y~x )
fit3 <- lm( y~poly(x,3) )
fit4 <- lm( y~poly(x,9) )
library(splines)
fit5 <- lm( y~ns(x, 3) )
fit6 <- lm( y~ns(x, 9) )
fit7 <- lm( y ~ x + cos(x*pi) )
xx <- seq(0,11, length.out=250)
lines(xx, predict(fit1,
data.frame(x=xx)), col='blue')
lines(xx, predict(fit2,
data.frame(x=xx)), col='green')
lines(xx, predict(fit3, data.frame(x=xx)),
col='red')
lines(xx, predict(fit4,
data.frame(x=xx)), col='purple')
lines(xx, predict(fit5,
data.frame(x=xx)), col='orange')
lines(xx, predict(fit6,
data.frame(x=xx)), col='grey')
lines(xx, predict(fit7,
data.frame(x=xx)), col='black')
polyfit <- function(i) x <-
AIC(lm(y~poly(x,i)))
as.integer(optimize(polyfit,interval =
c(1,length(x)-1))$minimum)
for (i in 2:length(x)-1)
print(polyfit(i))
lm(y ~ x + I(x^2) + I(x^3))
lm(y ~ poly(x, 3, raw=TRUE))
library(mice)
Mousey = mice(Greenwich)
Greenwich = complete(Mousey)
misp=read.csv(file("clipboard"),
header=T, sep="\t")
str(misp)
misp=misp[,-c(1:3)]
library(naniar)
vis_miss(misp)
# PROGRAM CLIMATE CHANGE
WITH R
#
abdi-basid@outlook.com
#
#
#
#
0> >> <<
Importation Data
rm(list=ls())
djibclim=read.table(file("clipboard"),header=T,
dec=".",sep="\t")
#File Change dir
#getwd()
#setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/A.ClimDjibouti/0.Data")
#save(djiclim, file =
"djiclim.rda")
#load(file = "djibclim.rda")
#View(djibclim)
#str(djibclim);summary(djiclim)
#sapply(djibclim,sd)
#sapply(djibclim, mean, na.rm=TRUE)
#library(psych)
#describe(djiclim[,4:6])
#Different Plot Types:
#plot(r,type="l",col=2)
#plot(r,type="b")
#plot(r,type="p")
#plot(r,type="o")
#plot(r,type="h")
#plot(r,type="s", col=4)
#plot(r,type="c")
plot(r,type="o",
main="", xlab="Year", ylab="Precipitation (mm)")
attach(r)
summary(lm(prectot~Year, data=r))
#
#
#
1> >> <<
source import Rclimtool
#File change dir
#getwd()
setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project
2")
source("RClimTool\\RClimTool.R")
#source("RClimTool\\RClimTool_v2.R")
#folder na gen, na, missing before
homogeneous
#day month year without label
#
#tmax = cargar(gfile(type =
"open"))
#tmin = cargar(gfile(type =
"open"))
#precip = cargar(gfile(type =
"open"))
#
#plot(x = day, main = "Title",
xlab = "Label for the x axis", ylab = "Label for the y
axis")
#hist(x = day, main = "Title",
xlab = "Label for the x axis", ylab = "Label for the y
axis")
#boxplot(x = day, main =
"Title", ylab = "Label for the y axis")
#
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\Rclimtrends\\R")
#library(climtrends)
source("climtrends.R")
names(djibclim)
attach(djibclim)
CraddockTest(Days[,Temp_Tn])
#
# Craddock test for testing the
homogeneity of yearly average temperatures of Milan vs Turin
data(yearly.average.temperature.Turin.Milan)
craddock.result <-
CraddockTest(yearly.average.temperature.Turin.Milan,2,3)
plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')
#
# Example of inhomogeneity
tmp <-
yearly.average.temperature.Turin.Milan
tmp[20:43,3] <- tmp[20:43,3]+1 #
adding an error of +1 degree to Milan's yearly average temperatures from 1980
craddock.result <-
CraddockTest(tmp,2,3)
dev.new()
plot(yearly.average.temperature.Turin.Milan[,1],craddock.result,type='l',xlab='',ylab='Craddock')
#
#Homogénéity Analysis
names(djiclim)
library(ggpubr)
ggqqplot(djiclim$Rain_Rrr)
#
library("car")
qqPlot(djiclim$Rain_Rrr)
#
#shapiro.test(rnorm(100, mean = 5, sd =
3))
shapiro.test(runif(100, min = 2, max =
4))
#require(stats)
shapiro.test(djibclim$Rain_Rrr)
#
library(tseries)
jarque.bera.test(djibclim$Rain_Rrr)
#
#Kolmogorov-Smirnov Tests
t.test(djibclim$Rain_Rrr)
wilcox.test(djibclim$Rain_Rrr)
ks.test(djibclim$Rain_Rrr)
#
#
#
2> >> << Program Rclim Dex
#library("RClimDex")
#rclimdex.start()
getwd()
setwd("D:/Docs.Abdi-Basid.ADAN/14.Doc.LMME/2.Project
2/Climat/RclimDex")
source("1.0.RclimDex.R.txt")
require("climdex.pcic")
require("Kendall")
require(tcltk)
#year month day without label
#I+1
i-1
#
#
#
3> >> << Program ClimPact2
#Download the above zip file to your
computer
#(https://github.com/ARCCSS-extremes/climpact2/archive/master.zip)
# and extract it. This will create a
directory named climpact2-master.
#In Windows: open R and select
"File -> Change dir..."
#and select the climpact2-master
directory created in step 1.
#Then type source('climpact2.GUI.r').
#(NOTE: If nothing happens, try run the
additional command startss()).
#In Linux/macOS: change to the
climpact2-master directory created in step 1,
#then open R in a terminal window and
type source('climpact2.GUI.r').
#The first time Climpact2 is run it will
install required R packages.
#This will likely require you to select
a mirror from which to download.
getwd()
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\2.Project
2\\Climat\\RClimPact\\climpact2-master")
source("climpact2.GUI.r")
#
#
4> >> << Program RHtest
#setwd("D:\\Docs.Abdi-Basid.ADAN\\Doc.LMME\\RHtest")
#source ("RHtest_V4.r")
#source ("RHtest.txt")
setwd("D:\\Docs.Abdi-Basid.ADAN\\14.Doc.LMME\\RHtest\\RHtests-master")
source("RHtestsV4_20180301.r")
StartGUI()
#
4'> >> << Program RHtest
library(PCICt)
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
CDD index.
cdd <- climdex.cdd(ci)
#
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
cold spell duration index.
csdi <- climdex.csdi(ci)
#
## Create a climdexInput object from
some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
tmax.dates <- as.PCICt(do.call(paste,
ec.1018935.tmax[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
tmin.dates <- as.PCICt(do.call(paste,
ec.1018935.tmin[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
prec.dates <- as.PCICt(do.call(paste,
ec.1018935.prec[,c("year",
"jday")]), format="%Y
%j", cal="gregorian")
## Load the data in.
ci <-
climdexInput.raw(ec.1018935.tmax$MAX_TEMP,
ec.1018935.tmin$MIN_TEMP,
ec.1018935.prec$ONE_DAY_PRECIPITATION,
tmax.dates, tmin.dates, prec.dates,
base.range=c(1971, 2000))
## Create an annual timeseries of the
CWD index.
cdd <- climdex.cdd(ci)
#
#
5> >> <<
ANCLIM
at=read.table(file("clipboard"),header=T,
sep="\t")
write(at, file = "data.txt",
sep = " ")
#
6> >> <<
ANNEXE
pairs(datatable, col="blue",
main="Scatterplots")
Y=cbind(LRY)
X=cbind(LRV, LRC, INT)
#
hist(Y, prob=TRUE, col =
"blue", border = "black")
lines(density(Y))
#
OLSreg=lm(Y~X)
summary(OLSreg)
#
Qreg25=rq(Y~X, tau=0.25)
summary(Qreg25)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
anova(Qreg25, Qreg75)
#
QR=rq(Y~X, tau=seq(0.2, 0.8, by=0.1))
sumQR=summary(QR)
#
plot(sumQR)
r=read.table(file("clipboard"),header=T,
sep="\t",dec=",",row.names=1)
str(r)
attach(r)
names(r)
dim(r)
summary(r)
library(Hmisc)
describe(r,num.desc=c("mean","median","var","sd","valid.n"),horizontal=TRUE)
sapply(r, sd)
sapply(r,range)
t=na.omit(r)
round(cor(t),3)
#e=scale(r, center=TRUE, scale=T)
library(caret)
preproc1 <- preProcess(r,
method=c("center", "scale"))
n <- predict(preproc1, r);summary(n)
library(corrplot)
mydata.cor = cor(t, method =
c("pearson"))
corrplot(mydata.cor, method =
"number", type = "lower")
#methode = number, ellipse, square,
shade, color, pie,
#type = full" "upper"
"lower"
library(RColorBrewer)
corrplot(mydata.cor, type =
"upper", order = "hclust",col = brewer.pal(n = 8, name =
"PuOr"))
palette =
colorRampPalette(c("yellow", "orange", "red"))
(20)
heatmap(x = mydata.cor, col = palette,
symm = TRUE)
#https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html
names(t)
#http://www.sthda.com/french/wiki/parametres-graphiques
s=read.table(file("clipboard"),header=T,
sep="\t",dec=",",row.names=1)
t=na.omit(s)
attach(t)
reg1=lm(OB~CH,data=t)
reg2=lm(OB~PE,data=t)
reg3=lm(OB~TA,data=t)
reg4=lm(OB~AR,data=t)
reg5=lm(OB~ER,data=t)
limx=c(0,500);limy=c(0,900) #windows(height = 70, width = 70)
plot(-20,-20,xaxt="n",yaxt="n",xlim=limx,ylim=limy,cex.lab=1.2,xlab="Ground
Based (mm/month)",ylab="Remotely sensed (mm/month)")
axis(1,cex.axis=1.3);axis(2,cex.axis=1.3)
par(new=TRUE)
plot(OB~CH,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="red",xlab="",ylab="");abline(reg1,col="red",lwd=2)
par(new=TRUE)
plot(OB~PE,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="blue",xlab="",ylab="");abline(reg2,col="blue",lwd=2)
par(new=TRUE)
plot(OB~TA,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="brown",xlab="",ylab="");abline(reg3,col="brown",lwd=2)
par(new=TRUE)
plot(OB~AR,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="orange",xlab="",ylab="");abline(reg4,col="yellow",lwd=2)
par(new=TRUE)
plot(OB~ER,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="darkgreen",xlab="",ylab="");abline(reg5,col="turquoise1",lwd=2)
par(new=TRUE)
abline(0,1.33,col="black",lwd=2,lty=2)
legend(30,400, legend=c("CH",
"PE", "TA", "AR", "ER"),horiz=T,bty
="n",
col=c("red","blue","brown","orange","darkgreen"),
lty=c(1,1,1,1,1), lwd=c(3,3,3,3,3), cex=0.8)
summary(reg1)
summary(reg2)
summary(reg3)
summary(reg4)
summary(reg5)
require(hydroGOF)
names(r);attach(r)
rmsec=sqrt(sum((OB-CH)^2)/length(OB));round(rmsec,3)
rmsec=sqrt(sum((OB-PE)^2)/length(OB));round(rmsec,3)
rmsec=sqrt(sum((OB-TA)^2)/length(OB));round(rmsec,3)
rmsec=sqrt(sum((OB-AR)^2)/length(OB));round(rmsec,3)
rmsec=sqrt(sum((OB-ER)^2)/length(OB));round(rmsec,3)
#rmse(sim=Data_Mon_11.325,
obs=Datat_Obser)
require(hydroGOF)
PBIAS = 100 * (abs(sum(OB-CH)) /
sum(OB));round(PBIAS,3)
PBIAS = 100 * (abs(sum(OB-CH)) /
sum(OB));round(PBIAS,3)
PBIAS = 100 * (abs(sum(OB-CH)) /
sum(OB));round(PBIAS,3)
PBIAS = 100 * (abs(sum(OB-CH)) /
sum(OB));round(PBIAS,3)
PBIAS = 100 * (abs(sum(OB-CH)) /
sum(OB));round(PBIAS,3)
#pbias(sim=CH, obs=OB, na.rm=TRUE)
require(hydroGOF)
# mse(CH, OB)
# mae(CH, OB)
gof(CH, OB)
ggof(sim=CH, obs=OB,
ftype="dm", FUN=mean)
#
STDE / R2
/ SLOPE /
r / CD / CV RSD
################################################
r=read.table(file("clipboard"),header=T,
sep="\t",dec=",",row.names=1)
attach(r)
names(r)
# "OB" "CH"
"PE" "TA" "AR" "ER"
round(sd(PE),3)
reg=lm(OB~ER,data=r);summary(reg)
round(cor(r),3)
BIAS =
sum(OB-TA)/length(OB);round(BIAS,3)
RMSE=sqrt(sum((OB-CH)^2)/length(OB));round(RMSE,3)
RMSE=sqrt(sum((OB-PE)^2)/length(OB));round(RMSE,3)
RMSE=sqrt(sum((OB-TA)^2)/length(OB));round(RMSE,3)
RMSE=sqrt(sum((OB-AR)^2)/length(OB));round(RMSE,3)
RMSE=sqrt(sum((OB-ER)^2)/length(OB));round(RMSE,3)
RMSD =
sqrt(sum((R$OB-R$PE)^2)/length(R$OB));round(RMSD,3)
names(n)
#########
[1] "JF.OB" "JF.CH" "JF.PE" "JF.TA" "JF.AR" "JF.ER" "MAM.OB"
[8] "MAM.CH" "MAM.PE" "MAM.TA" "MAM.AR" "MAM.ER" "JJAS.OB" "JJAS.CH"
[15] "JJAS.PE"
"JJAS.TA" "JJAS.AR" "JJAS.ER"
"OND.OB"
"OND.CH"
"OND.PE"
[22] "OND.TA" "OND.AR" "OND.ER"
##########
RSD = abs(sd(CH)/mean(CH));round(RSD,
digits = 3)
###########################################
# display the diagram with the better
model
#oldpar<-taylor.diagram(OB,CH,pch=19)
#taylor.diagram.modified(OB,CH,
text="Model 1")
s=read.table(file("clipboard"),header=T,
sep="\t",dec=",",row.names=1)
attach(s)
names(s)
str(s)
cbind(sapply(s, sd, na.rm=TRUE))
CHrmse=sqrt(sum((OB-CH)^2)/length(OB));round(CHrmse,3)
PErmse=sqrt(sum((OB-PE)^2)/length(OB));round(PErmse,3)
TArmse=sqrt(sum((OB-TA)^2)/length(OB));round(TArmse,3)
ARrmse=sqrt(sum((OB-AR)^2)/length(OB));round(ARrmse,3)
ERrmse=sqrt(sum((OB-ER)^2)/length(OB));round(ERrmse,3)
rbind(CHrmse,PErmse,TArmse,ARrmse,ERrmse)
OB=AIRPORT
CH=CHIRPS
PE=PERSIANN
TA=TAMSATV3
AR=ARCV2
ER=ERA5
#
par(new=TRUE)
require(plotrix)
oldpar<-taylor.diagram(OB,CH,add=F,
pch=19,pos.cor=TRUE,
xlab="Standard
deviation",ylab="Standard Deviation",main="",
show.gamma=TRUE,ngamma=10,gamma.col="green",sd.arcs=1,
ref.sd=TRUE,sd.method="sample",grad.corr.lines=c(0.1,
0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.99),col="blue",
pcex=1.5,cex.axis=1.1,normalize=F)
# now add the worse model
taylor.diagram(OB,PE,add=TRUE,col="black",
pcex=1.5,cex.axis=1.1,normalize=F)
taylor.diagram(OB,TA,add=TRUE,col="pink",
pcex=1.5,cex.axis=1.1,normalize=F)
taylor.diagram(OB,AR,add=TRUE,col="brown",
pcex=1.5,cex.axis=1.1,normalize=F)
taylor.diagram(OB,ER,add=TRUE,col="green",
pcex=1.5,cex.axis=1.1,normalize=F)
# get approximate legend position # add
a legend
legend(1.3,1.6,legend=c("JF.OB","JF.CH","JF.PE","JF.TA","JF.AR","JF.ER"),horiz=FALSE,
pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","orange","darkgreen"),
cex=0.7)
# SAISONNALITY
#par(mfrow=c(1,2))
#"red","blue","brown","orange","darkgreen"
s=read.table(file("clipboard"),header=T,
sep="\t",dec=",",row.names=1)
names(s)
attach(s)
cor(s)
str(s)
# show the "all correlation"
display
require(plotrix)
taylor.diagram(JF.OB, JF.CH, col="blue", pos.cor=F,pcex=1.5,normalize=F)
taylor.diagram(JF.OB, JF.PE,add=TRUE, col="black", pcex=1.5,normalize=F)
taylor.diagram(JF.OB, JF.TA,add=TRUE, col="pink", pcex=1.5,normalize=F)
taylor.diagram(JF.OB, JF.AR,add=TRUE, col="brown", pcex=1.5,normalize=F)
taylor.diagram(JF.OB, JF.ER,add=TRUE, col="green", pcex=1.5,normalize=F)
require(plotrix)
taylor.diagram(MAM.OB, MAM.CH, col="blue", pos.cor=F,pcex=1.5,normalize=F)
taylor.diagram(MAM.OB, MAM.PE,add=TRUE, col="black", pcex=1.5,normalize=F)
taylor.diagram(MAM.OB, MAM.TA,add=TRUE, col="pink", pcex=1.5,normalize=F)
taylor.diagram(MAM.OB, MAM.AR,add=TRUE, col="brown", pcex=1.5,normalize=F)
taylor.diagram(MAM.OB, MAM.ER,add=TRUE, col="green", pcex=1.5,normalize=F)
require(plotrix)
taylor.diagram(JJAS.OB, JJAS.CH, col="blue", pos.cor=F,pcex=1.5,normalize=F)
taylor.diagram(JJAS.OB, JJAS.PE,add=TRUE, col="black", pcex=1.5,normalize=F)
taylor.diagram(JJAS.OB, JJAS.TA,add=TRUE, col="pink", pcex=1.5,normalize=F)
taylor.diagram(JJAS.OB, JJAS.AR,add=TRUE, col="brown", pcex=1.5,normalize=F)
taylor.diagram(JJAS.OB, JJAS.ER,add=TRUE, col="green", pcex=1.5,normalize=F)
require(plotrix)
taylor.diagram(OND.OB, OND.CH, col="blue", pos.cor=F,pcex=1.5,normalize=F)
taylor.diagram(OND.OB, OND.PE,add=TRUE, col="black", pcex=1.5,normalize=F)
taylor.diagram(OND.OB, OND.TA,add=TRUE, col="pink", pcex=1.5,normalize=F)
taylor.diagram(OND.OB, OND.AR,add=TRUE, col="brown", pcex=1.5,normalize=F)
taylor.diagram(OND.OB, OND.ER,add=TRUE, col="green", pcex=1.5,normalize=F)
legend(30,53
,legend=c("OND.OBSERVATON","OND.CHIRPS","OND.PESERIANNCDR","OND.TAMSAT","OND.ARC","OND.ERA"),horiz=FALSE,
pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","orange","darkgreen"),
cex=0.7)
CHrmse=sqrt(sum((OND.OB-OND.CH)^2)/length(OND.OB))
PErmse=sqrt(sum((OND.OB-OND.PE)^2)/length(OND.OB))
TArmse=sqrt(sum((OND.OB-OND.TA)^2)/length(OND.OB))
ARrmse=sqrt(sum((OND.OB-OND.AR)^2)/length(OND.OB))
ERrmse=sqrt(sum((OND.OB-OND.ER)^2)/length(OND.OB))
round(rbind(CHrmse,PErmse,TArmse,ARrmse,ERrmse),3)
#RMS_Diff =
sum(((OND.CH-mean(OND.CH))-(OND.OB-mean(OND.OB)))^2)/length(OND.OB)
OBsd=sd(OND.OB)
CHsd=sd(OND.CH)
PEsd=sd(OND.PE)
TAsd=sd(OND.TA)
ARsd=sd(OND.AR)
ERsd=sd(OND.ER)
round(rbind(CHsd,PEsd,TAsd,ARsd,ERsd),3)
CHBIAS =
sum(OND.OB-OND.CH)/length(OND.OB)
PEBIAS =
sum(OND.OB-OND.PE)/length(OND.OB)
TABIAS = sum(OND.OB-OND.TA)/length(OND.OB)
ARBIAS =
sum(OND.OB-OND.AR)/length(OND.OB)
ERBIAS =
sum(OND.OB-OND.ER)/length(OND.OB)
round(rbind(CHBIAS,PEBIAS,TABIAS,ARBIAS,ERBIAS),3)
R=read.table(file("clipboard"),header=T,
sep="\t",dec=".",row.names=1)
library(caret)
preproc1 <- preProcess(R,
method=c("center"))
#"center", "scale"
norm1 <- predict(preproc1,R)
print(norm1)
require(openair)
## in the examples below, most effort
goes into making some artificial data
## the function itself can be run very
simply
## Not run:
## dummy model data for 2003
dat <- selectByDate(mydata, year =
2003)
dat <- data.frame(date = mydata$date,
obs = mydata$nox, mod = mydata$nox)
## now make mod worse by adding bias and
noise according to the month
## do this for 3 different models
dat <- transform(dat, month =
as.numeric(format(date, "%m")))
mod1 <- transform(dat, mod = mod + 10
* month + 10 * month * rnorm(nrow(dat)),
model = "model 1")
## lag the results for mod1 to make the
correlation coefficient worse
## without affecting the sd
mod1 <- transform(mod1, mod =
c(mod[5:length(mod)], mod[(length(mod) - 3) :
length(mod)]))
## model 2
mod2 <- transform(dat, mod = mod + 7
* month + 7 * month * rnorm(nrow(dat)),
model = "model 2")
## model 3
mod3 <- transform(dat, mod = mod + 3
* month + 3 * month * rnorm(nrow(dat)),
model = "model 3")
mod.dat <- rbind(mod1, mod2, mod3)
## basic Taylor plot
TaylorDiagram(mod.dat, obs =
"obs", mod = "mod", group = "model")
## Taylor plot by season
TaylorDiagram(mod.dat, obs =
"obs", mod = "mod", group = "model", type =
"season")
## now show how to evaluate model
improvement (or otherwise)
mod1a <- transform(dat, mod = mod + 2
* month + 2 * month * rnorm(nrow(dat)),
model = "model 1")
mod2a <- transform(mod2, mod = mod *
1.3)
mod3a <- transform(dat, mod = mod +
10 * month + 10 * month * rnorm(nrow(dat)),
model = "model 3")
mod.dat2 <- rbind(mod1a, mod2a,
mod3a)
mod.dat$mod2 <- mod.dat2$mod
## now we have a data frame with 3
models, 1 set of observations
## and TWO sets of model predictions
(mod and mod2)
## do for all models
TaylorDiagram(mod.dat, obs =
"obs", mod = c("mod", "mod2"), group =
"model")
## End(Not run)
## Not run:
## all models, by season
TaylorDiagram(mod.dat, obs =
"obs", mod = c("mod", "mod2"), group =
"model",
type = "season")
## consider two groups (model/month). In
this case all months are shown by model
## but are only differentiated by model.
TaylorDiagram(mod.dat, obs =
"obs", mod = "mod", group = c("model",
"month"))
## End(Not run)
# ANNEXE
####################################################3
#Taylor, K.E. (2001) Summarizing
multiple aspects of model performance in a single diagram. Journal of
Geophysical Research, 106: 7183-7192.
library(datasets)
library(ncdf4)
library(plotrix)
taylor.diagram(r,r,add=FALSE,col="red",pch=4,pos.cor=TRUE,xlab="MERRA
SD (Normalised)",ylab="RCA4 runs SD
(normalised)",main="Taylor
Diagram",show.gamma=TRUE,ngamma=3,sd.arcs=1,ref.sd=TRUE,grad.corr.lines=c(0.2,0.4,0.6,0.8,0.9),pcex=1,cex.axis=1,normalize=TRUE,mar=c(5,4,6,6),lwd=10,font=5,lty=3)
lpos<-1.5*sd(Data_Mon_11.275)
legend(1.5,1.5,cex=1.2,pt.cex=1.2,legend=c("volcano"),pch=4,col=c("red"))
taylor.diagram(data,data,normalize=TRUE)
# fake some reference data
ref<-rnorm(30,sd=2)
#
add a little noise
model1<-ref+rnorm(30)/2
#
add more noise
model2<-ref+rnorm(30)
#
display the diagram with the better model
oldpar<-taylor.diagram(ref,model1)
#
now add the worse model
taylor.diagram(ref,model2,add=TRUE,col="blue")
#
get approximate legend position
lpos<-1.5*sd(ref)
#
add a legend
legend(lpos,lpos,legend=c("Better","Worse"),pch=19,col=c("red","blue"))
#
now restore par values
par(oldpar)
#
show the "all correlation" display
taylor.diagram(ref,model1,pos.cor=FALSE)
taylor.diagram(ref,model2,add=TRUE,col="blue")
FIN DU PROGRAMME
Enjoy It !
 |
| The Abdi-Basid Courses Institute |
