How to Diagnose Change in Temperature and Precipitation with the R
program ?
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
& CHANGE PIONT
#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 !
Aucun commentaire:
Enregistrer un commentaire