0% found this document useful (0 votes)

61 views50 pages

Analyzing Climate Change with R

Q: What methods are employed in the document for detecting homogeneity in climatic datasets and what are their implications?

The document employs methods such as the Craddock test and statistical visualization (like the QQ plot) to detect homogeneity in climatic datasets. The Craddock test examines the homogeneity by comparing climate data from different sources or periods, identifying significant deviations. Graphical methods such as QQ plots assess normality and outliers visually. Homogeneity testing is crucial for ensuring the reliability of comparative climatic analyses and time series modeling .

Q: How is the Mann-Kendall Trend Test used in detecting trends within hydroclimatic data, as illustrated in the document?

The Mann-Kendall Trend Test is used for detecting trends in hydroclimatic data by evaluating whether there is a statistically significant trend (either upward or downward) over time. In the document, it is applied to analyze precipitation data using `mk.test(prectot)`, which assesses the existence of a monotonic trend without assuming normally-distributed data. This test is especially useful in climate studies due to its robustness against non-normality and unequal variances in time series data .

Q: What is the significance of identifying breakpoints in time series analysis, and how is it applied in the context of the provided document?

Identifying breakpoints in time series analysis is crucial for understanding structural changes within the data, such as a shift in mean, variance, or trend. In the context of the provided document, breakpoints are detected using methods like OLS-CUSUM and MOSUM tests. For example, the `prectot` series is analyzed to find optimal number of breaks using the BIC criteria, where a significant breakpoint was identified around 1980, indicating a change in the pattern of annual precipitation .

Q: What role does the Seasonal Decomposition of Time Series (STL) play in understanding precipitation trends as per the document?

Seasonal Decomposition of Time Series (STL) plays a significant role in understanding precipitation trends by separating the time series into seasonal, trend, and irregular components. In the document, STL is applied using the function `stl(s.window='periodic')`, allowing for the adjustment of seasonal patterns. This decomposition aids in isolating and analyzing trends and seasonal effects within the precipitation data, providing clearer insights for forecasting and understanding underlying patterns .

Q: Explain the concept of model selection criteria like BIC in determining breakpoints in time series according to the document.

The concept of model selection criteria like the Bayesian Information Criterion (BIC) involves choosing models with optimal complexity—balancing fit and simplicity. In time series analysis, especially for identifying breakpoints as discussed in the document, BIC helps determine the number of breakpoints by evaluating the model's likelihood and penalizing for complexity. Lower BIC values indicate a preferable model, guiding the determination of structural changes in the time series data .

Q: Evaluate the effectiveness of using ARIMA models for forecasting precipitation trends, based on the techniques described in the document.

ARIMA models are highly effective for forecasting precipitation trends as they account for autoregressive and moving average components, and can adjust for trends and seasonality by differencing. The document details the steps, such as exploratory analysis, model fitting, and diagnostics like the Ljung-Box test, to ensure model adequacy. ARIMA modeling in the document includes fitting models with different orders to optimize prediction, indicating its adaptability and precision in forecasting climate variables .

Q: Analyze the application of regression models, such as linear regression, in studying seasonal climate variations mentioned in the document.

Linear regression models are applied to study seasonal climate variations by examining relationships between climate variables over time. In the document, regression models are fitted to months or seasonal periods (e.g., JFM) against time to estimate trends and assess statistical significance. For example, the equation y = -0.25x + 519.54 captures the trend in a season's data. These regression models allow for the quantification of trends, aiding in understanding seasonal climate behavior and its changes .

Q: How is the Pettitt's test applied to detect change points in a precipitation dataset, and what were the findings according to the document?

Pettitt's test is a non-parametric test used to detect change points in a time series by identifying the point where there is a significant change in the median. In the provided document, Pettitt's test is applied to the `Tave[['prectot']]` dataset to detect change points. The result of the test indicated a change point, which was further validated by plotting with a red vertical line on the graph to highlight the estimated change point .

Q: What are the implications of using the Cox and Stuart Trend test in historical precipitation analysis as mentioned in the document?

The Cox and Stuart Trend test is used to detect monotonic trends in historical data without assuming normality. In historical precipitation analysis, its application as described in the document provides insights into long-term trends by identifying whether there has been a significant and consistent increase or decrease over time. This non-parametric test is especially valuable for climate data, where assumptions of normality can be violated. Such analyses guide policy and strategy for future climate events by revealing systematic trends .

Q: Discuss the importance and application of Quantile Regression in analyzing climate data, with reference to the document's context.

Quantile Regression is important for analyzing climate data because it provides a more comprehensive analysis by estimating the relationship between variables at different quantiles. In the document, quantile regression is used on precipitation data with `rq(prectot ~ Year, tau = 0.5)` to capture the median trend over the years, while also allowing for a comparison with the null model using an ANOVA test. This approach captures variations that might not be evident through the mean, especially relevant in climate studies where extreme weather events are of interest .

This document discusses using R programs to analyze climate change by diagnosing changes in temperature and precipitation. It summarizes four main R programs: 1. RClimDex - Calculates 27 climate change indices for temperature and precipitation and provides graphs and trend analysis. 2. RClimTool - Performs quality control, identifies missing data, and tests for data homogeneity. 3. RClimTrend - Applies statistical tests like Mann-Kendall to detect trends and changes in climate data. 4. ClimPact - Evaluates climate extremes and performs additional statistical analyses like boxplots to compare results to RClimDex.

Uploaded by

Abdi-Basid ADAN

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as DOCX, PDF, TXT or read online on Scribd

0% found this document useful (0 votes)

61 views50 pages

Analyzing Climate Change with R

Uploaded by

Abdi-Basid ADAN

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as DOCX, PDF, TXT or read online on Scribd

How to Diagnose Change in Temperature and

Precipitation with the R program ?

IntroductionThe 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, [Link].). 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

Max 5-day precipitation amount Monthly maximum consecutive 5-day

Rx5day (mm)
precipitation

Simple daily intensity index Annual total precipitation divided by the

SDII (mm/day) number of wet days (defined as
PRCP>=1.0mm) in the year

Number of heavy precipitation days Annual count of days when

R10 (day)
PRCP>=10mm

Number of very heavy precipitation days Annual count of days when

R20 (day)
PRCP>=20mm

Number of days above nn mm Annual count of days when PRCP>=nn

Rnn (day)
mm, nn is user defined threshold

Consecutive dry days Maximum number of consecutive days

CDD (day)
with RR<1mm

Consecutive wet days Maximum number of consecutive days

CWD (day)
with RR>=1mm

Very wet days Annual total PRCP when RR>95th

R95p (mm)
percentile

Extremely wet days Annual total PRCP when RR>99th

R99p (mm)
percentile

Annual total wet-day precipitation Annual total PRCP in wet days

PRCPTOT (mm) (RR>=1mm)
R Program Analysis Climat Change
TREND & CHANGE PIONT

#TREND ANANLYSIS
data=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)
data
str(data)

#[Link]("trend")
require(trend)

#data(maxau)
prectot <- data[,"prectot"]
[Link](prectot)

dat=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)

data(nottem)
str(dat)
dat=round(dat,1)
dat1 <- ts(dat, frequency=12, start=c(1961,1))
[Link](dat1)

logsouvenirtimeseries <- log(dat1)

[Link](logsouvenirtimeseries)

require(trend)
[Link](dat1)

#correlated witn precedent month

[Link](dat)

#Multivariate Mann-Kendall Test

require(trend)
data(hcb)
> plot(hcb)

s <- maxau[,"s"]
Q <- maxau[,"Q"]
[Link](s,Q, meth="spearman")

[Link](s,Q)
[Link](s,Q, "spearman")

#Cox and Stuart Trend test

[Link](prectot)

s <- maxau[,"s"]
[Link](prectot)

#Pettitt's test for single change-point detection

[Link](dat)

#Buishand range test

require(trend)
(res <- [Link](prectot))

require(trend)
(res <- [Link](Nile))

#Wallis and Moore Phase-Frequency test

frost <- ts(prectot, start=1961)
[Link](frost)

#Bartels's test for randomness

[Link](prectot)

#Wald-Wolfowitz test for independence and stationarity

[Link](prectot)

#Wallis and Moore Phase-Frequency test

[Link](prectot)

#Buishand U test
#Nile
(res <- [Link](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 <- [Link](prectot))

[Link]( x,y,alternative = c("[Link]", "less", "greater"),

use="[Link]", method="kendall",[Link] = 0.95)

require(graphics)
x <- rnorm(50)
y <- runif(30)
# Do x and y come from the same distribution?
[Link](x, y)
# Does x come from a shifted gamma distribution with shape 3 and
rate 2?
[Link](x+2, "pgamma", 3, 2) # two-sided, exact
[Link](x+2, "pgamma", 3, 2, exact = FALSE)
[Link](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, [Link] = FALSE )
library(package = ggplot2, [Link] = FALSE)
library(package = trend, [Link] = FALSE)

# Load Data
data <- [Link](file = 'c:/Users/pooya/Downloads/Torbat Heydariyeh
- [Link]',
header = TRUE)
data=[Link](file("clipboard"),header=T, sep="\t", [Link]=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 = [Link](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 <- [Link](Tave)
# Plot
require(ggplot2)
ggplot(data = Tave, mapping = aes(x = date, y = prectot)) +
geom_line()

# Pettitt’s test
pettittTest <- trend::[Link](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 = [Link](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=[Link](file("clipboard"),header=T, sep="\t", [Link]=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)
[Link] <- efp(prectot ~ 1, type = "OLS-CUSUM")
[Link] <- efp(prectot ~ 1, type = "OLS-MOSUM")
[Link] <- 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([Link], type="o",xlab="Year")
plot([Link], type="o",xlab="Year")
plot([Link], 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:

[Link] =
Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)
plot([Link])

#use the breakpoints functions to find BIC corresponding to optimal

number of #breaks
[Link] = breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)
summary([Link])
x11()
plot([Link])

prectot=ts(prectot, start=1961, end=2017)

[Link] <- breakpoints(prectot ~ 1)
[Link] <- breakfactor([Link], breaks = 1 )
[Link] <- lm(prectot~ [Link] - 1)
plot([Link])

prectot=ts(prectot, start=1961, end=2017)

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))
plot([Link], [Link]=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([Link]),start=1980,freq=1),
col='darkgreen',lwd=2)
par(opar)

plot(prectot)

## F statistics indicate one breakpoint

[Link] <- Fstats(prectot ~ 1)
plot([Link])
breakpoints([Link])
lines(breakpoints([Link]))

## or
[Link] <- breakpoints(prectot ~ 1)
summary([Link])

## the BIC also chooses one breakpoint

plot([Link])
breakpoints([Link])

## confidence interval
[Link] <- confint([Link])
[Link]
lines([Link])

## fit null hypothesis model and model with 1 breakpoint

fm0 <- lm(prectot ~ 1)
fm1 <- lm(prectot ~ breakfactor([Link], 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([Link])
abline(v=2006, lty=2, col=2)

## confidence intervals
ci.seat2 <- confint([Link], breaks = 2)
ci.seat2
lines(ci.seat2)

#rm(list=ls())
data=[Link](file("clipboard"),header=T, sep="\t")
str(data)

if(!require(mblm)){[Link]("mblm")}
if(!require(ggplot2)){[Link]("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)){[Link]("quantreg")}
library(quantreg)

model.q = rq(prectot ~ Year, data = data, tau = 0.5)

summary(model.q)

library(ggplot2)
[Link] = rq(prectot ~ 1, data = data, tau = 0.5)
anova(model.q, [Link])

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=[Link](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=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)

getwd()
setwd("C:/Users/Lenovo/Documents")
[Link](rm, file="[Link]",[Link]=F)

rm=[Link](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([Link]~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","b
lack"),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",[Link]="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)
[Link](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 = [Link](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 = [Link](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)
)

[Link] <- efp(Nile ~ 1, type = "OLS-CUSUM")

[Link] <- efp(Nile ~ 1, type = "OLS-MOSUM")
[Link] <- efp(Nile ~ 1, type = "Rec-CUSUM")

opar <- par(mfrow=c(2,2))

plot([Link])
plot([Link])
plot([Link])
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")

[Link] <- breakpoints(Nile ~ 1)

[Link] <- breakfactor([Link], breaks = 1 )
[Link] <- lm(Nile ~ [Link] - 1)
plot([Link])

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

plot([Link], [Link]=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([Link]),start=1871,freq=1),
col='darkgreen',lwd=2)
par(opar)

#TREND ANANLYSIS
#################
data=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)
data
str(data)

#[Link]("trend")
require(trend)

#data(maxau)
prectot <- data[,"prectot"]
[Link](prectot)

dat=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)

data(nottem)
str(dat)
dat=round(dat,1)
dat1 <- ts(dat, frequency=12, start=c(1961,1))
[Link](dat1)

logsouvenirtimeseries <- log(dat1)

[Link](logsouvenirtimeseries)

require(trend)
[Link](dat1)

#correlated witn precedent month

[Link](dat)

#Multivariate Mann-Kendall Test

require(trend)
data(hcb)
> plot(hcb)

s <- maxau[,"s"]
Q <- maxau[,"Q"]
[Link](s,Q, meth="spearman")

[Link](s,Q)
[Link](s,Q, "spearman")

#Cox and Stuart Trend test

[Link](prectot)

s <- maxau[,"s"]
[Link](prectot)

dat=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)

require(trend)
das=ts(dat, frequency=12,start=1961)
[Link](das)

# Changement Dectection
#######################
#data(PagesData)
require(trend)
data=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)
attach(data)

#Pettitt's test for single change-point detection

[Link](dat)

#Buishand range test

require(trend)
(res <- [Link](prectot))

require(trend)
(res <- [Link](Nile))

#Wallis and Moore Phase-Frequency test

frost <- ts(prectot, start=1961)
[Link](frost)

#Bartels's test for randomness

[Link](prectot)

#Wald-Wolfowitz test for independence and stationarity

[Link](prectot)

#Wallis and Moore Phase-Frequency test

[Link](prectot)

#Buishand U test
#Nile
(res <- [Link](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 <- [Link](prectot))

[Link]( x,y,alternative = c("[Link]", "less", "greater"),

use="[Link]", method="kendall",[Link] = 0.95)

# The approach after Pettitt (1979)

The Pettitt-test is conducted in such a way

# Change-point Detection: Pettitt’s Test

library(package = dplyr, [Link] = FALSE )
library(package = ggplot2, [Link] = FALSE)
library(package = trend, [Link] = FALSE)

# Load Data
data <- [Link](file = 'c:/Users/pooya/Downloads/Torbat Heydariyeh
- [Link]',
header = TRUE)
data=[Link](file("clipboard"),header=T, sep="\t", [Link]=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 = [Link](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 <- [Link](Tave)

# Plot
require(ggplot2)
ggplot(data = Tave, mapping = aes(x = date, y = prectot)) +
geom_line()

# Pettitt’s test
pettittTest <- trend::[Link](x = Tave[['prectot']])
print(pettittTest)

print(Tave[['date']][pettittTest$estimate])

# 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=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)
attach(ry)
str(ry)

opar <- par(mfrow=c(2,2))

[Link] =
Fstats(year~Dec+Feb+Mar+Apr+May+Jun+Jul+Aug+Sep+Oct+Nov,data=rm)
plot([Link])

#use the breakpoints functions to find BIC corresponding to optimal

number of #breaks
[Link] = breakpoints(year~Dec+Feb+Mar+Apr+May,data=rm)
summary([Link])
x11()
plot([Link])

prectot=ts(prectot, start=1961, end=2017)

[Link] <- breakpoints(prectot ~ 1)
[Link] <- breakfactor([Link], breaks = 1 )
[Link] <- lm(prectot~ [Link] - 1)
plot([Link])

prectot=ts(prectot, start=1961, end=2017)

plot(prectot)

## F statistics indicate one breakpoint

[Link] <- Fstats(prectot ~ 1)
plot([Link])
breakpoints([Link])
lines(breakpoints([Link]))

## or
[Link] <- breakpoints(prectot ~ 1)
summary([Link])

## the BIC also chooses one breakpoint

plot([Link])
breakpoints([Link])

## confidence interval
[Link] <- confint([Link])
[Link]
lines([Link])

## fit null hypothesis model and model with 1 breakpoint

## confidence intervals
ci.seat2 <- confint([Link], breaks = 2)
ci.seat2
lines(ci.seat2)

#rm(list=ls())
data=[Link](file("clipboard"),header=T, sep="\t")
str(data)

if(!require(mblm)){[Link]("mblm")}
if(!require(ggplot2)){[Link]("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)){[Link]("quantreg")}
library(quantreg)

model.q = rq(prectot ~ Year, data = data, tau = 0.5)

summary(model.q)

library(ggplot2)
[Link] = rq(prectot ~ 1, data = data, tau = 0.5)
anova(model.q, [Link])
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=[Link](file("clipboard"),header=T, dec=".",sep="\t")
#File Change dir
#getwd()
#setwd("D:/[Link]/[Link]/[Link]/[Link]")
#save(djiclim, file = "[Link]")
#load(file = "[Link]")
#View(djibclim)
#str(djibclim);summary(djiclim)
#sapply(djibclim,sd)
#sapply(djibclim, mean, [Link]=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:/[Link]/[Link]/[Link] 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:\\[Link]\\[Link]\\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([Link])
[Link] <-
CraddockTest([Link],2,3)
plot([Link][,1],[Link],type
='l',xlab='',ylab='Craddock')
#
# Example of inhomogeneity
tmp <- [Link]
tmp[20:43,3] <- tmp[20:43,3]+1 # adding an error of +1 degree to
Milan's yearly average temperatures from 1980
[Link] <- CraddockTest(tmp,2,3)
[Link]()
plot([Link][,1],[Link],type
='l',xlab='',ylab='Craddock')
#
#Homogénéity Analysis
names(djiclim)
library(ggpubr)
ggqqplot(djiclim$Rain_Rrr)
#
library("car")
qqPlot(djiclim$Rain_Rrr)
#
#[Link](rnorm(100, mean = 5, sd = 3))
[Link](runif(100, min = 2, max = 4))
#require(stats)
[Link](djibclim$Rain_Rrr)
#
library(tseries)
[Link](djibclim$Rain_Rrr)
#
#Kolmogorov-Smirnov Tests
[Link](djibclim$Rain_Rrr)
[Link](djibclim$Rain_Rrr)
[Link](djibclim$Rain_Rrr)
#
#
#
# 2> >> << Program Rclim Dex
#library("RClimDex")
#[Link]()
getwd()
setwd("D:/[Link]/[Link]/[Link]
2/Climat/RclimDex")
source("[Link]")
require("[Link]")
require("Kendall")
require(tcltk)
#year month day without label
#I+1 i-1
#
# 3> >> << Program ClimPact2
#Download the above zip file to your computer
#([Link]
# 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('[Link].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('[Link].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:\\[Link]\\[Link]\\[Link] 2\\Climat\\
RClimPact\\climpact2-master")
source("[Link].r")
#
#
# 4> >> << Program RHtest
#setwd("D:\\[Link]\\[Link]\\RHtest")
#source ("RHtest_V4.r")
#source ("[Link]")
setwd("D:\\[Link]\\[Link]\\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.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the CDD index.
cdd <- [Link](ci)
#
## Create a climdexInput object from some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the cold spell duration index.
csdi <- [Link](ci)
#
#
#
## Create a climdexInput object from some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the CWD index.
cdd <- [Link](ci)

#
# 5> >> << ANCLIM

at=[Link](file("clipboard"),header=T, sep="\t")
write(at, file = "[Link]", sep = " ")

# 5> >> << Geoclim

# 6> >> << GeoDa

# ANNEXE
########

require(ggplot2)
ggplot(data = quakes) +
geom_density(mapping = aes(x = mag, fill = region), alpha =
0.5)

library("RColorBrewer")
[Link]()

# Missing Data Imputation

#########################
library(mice)
Mousey = mice(Greenwich)
Greenwich = complete(Mousey)

misp=[Link](file("clipboard"), header=T, sep="\t")

str(misp)
misp=misp[,-c(1:3)]
library(naniar)
vis_miss(misp)

#Prediction
###########

#Import Data
data=[Link](file("clipboard"), header=T, sep="\t", dec=".")
str(data)
summary(data)
library(forecast)

fit.series3 <- [Link](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 = [Link](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")

#[Link]
latest/src/[Link]
#[Link]

prectot %>% stl([Link]='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

[Link] = decompose(tsData)
plot([Link])

library("fUnitRoots")
urkpssTest(tsData, type = c("tau"), lags = c("short"),[Link] =
NULL, doplot = TRUE)
tsstationary = diff(tsData, differences=1)
plot(tsstationary)

acf(tsstationary, [Link]=34)
pacf(tsstationary, [Link]=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)
[Link](tsData, trace=TRUE)

predict(fitARIMA,[Link] = 20)
require(forecast)
futurVal <- [Link](fitARIMA,h=10, level=c(99.5))
[Link](futurVal)
plot(fitARIMA,h=10, level=c(99.5))

A=[Link](prectot)
autoplot(forecast(A))

autoplot(forecast(prectot), n=100)
library(forecast)
AutoArimaModel=[Link](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, [Link] = 50)

fit = arima(log(precto), c(1, 1, 0), seasonal = list(order = c(0, 1,

1), period = 12))
pred <- predict(fit, [Link] = 10*12)
[Link](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)), [Link] = 12)

(fit <- arima(USAccDeaths, order = c(0,1,1),

seasonal = list(order = c(0,1,1))))

predict(fit, [Link] = 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, [Link]="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 <- [Link](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 [Link]() 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
#[Link](ts) Augemented Dickey-Fuller test. Rejecting the null
hypothesis suggests that a time series is stationary (from the
tseries package)
#[Link](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, [Link]=250)

lines(xx, predict(fit1, [Link](x=xx)), col='blue')
lines(xx, predict(fit2, [Link](x=xx)), col='green')
lines(xx, predict(fit3, [Link](x=xx)), col='red')
lines(xx, predict(fit4, [Link](x=xx)), col='purple')
lines(xx, predict(fit5, [Link](x=xx)), col='orange')
lines(xx, predict(fit6, [Link](x=xx)), col='grey')
lines(xx, predict(fit7, [Link](x=xx)), col='black')

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))

[Link](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
[Link]=nnet(O3obs~.,data=datappr,size=5,decay=1,linout=TRUE,maxit=
500)
summary([Link])

library(e1071)
plot([Link](O3obs~.,data=datappr,size=c(2,3,4),
decay=c(1,2,3),maxit=200,linout=TRUE))
plot([Link](O3obs~.,data=datappr,size=4:5,decay=1:10)

[Link]=nnet(O3obs~.,data=datappr,size=3,decay=2,
linout=TRUE,maxit=200)
# calcul et graphe des résidus
[Link]=predict([Link],data=datappr)
[Link]=[Link]-datappr[,"O3obs"]
[Link]([Link],[Link])

# apprentissage
[Link]=nnet(DepSeuil~.,data=datappq,size=5,decay=1)
summary([Link])
#matrice de confusion
table([Link]$[Link]>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

[Link]=predict([Link],newdata=datestr)
[Link]=predict([Link],newdata=datestq)
# Erreur quadratique moyenne de prévision
sum(([Link]-datestr[,"O3obs"])^2)/nrow(datestr)
# Matrice de confusion pour la prévision du
# dépassement de seuil (régression)
table([Link]>150,datestr[,"O3obs"]>150)
# Même chose pour la discrimination
table([Link]>0.5,datestq[,"DepSeuil"])

library(ROCR)
rocnnetr=[Link]/300
prednnetr=prediction(rocnnetr,datestq$DepSeuil)
perfnnetr=performance(prednnetr,"tpr","fpr")
rocnnetq=[Link]
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=[Link](file("clipboard"),header=TRUE,sep="\t", [Link]=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 <- "[Link]
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

[Link](precipitation_xts, r = c(30, 45), nint = 16)
# ismev implementation is faster:
# ismev::[Link](precipitation_xts, umin=30, umax=45, nint =
16)
# set threshold
th <- 40

# maximum likelihood estimation

pot_mle <- fevd([Link](precipitation_xts), method = "MLE",
type="GP", threshold=th)
# diagnostic plots
plot(pot_mle)
rl_mle <- [Link](pot_mle, conf = 0.05, [Link]=
c(2,5,10,20,50,100))

# L-moments estimation
pot_lmom <- fevd([Link](precipitation_xts), method = "Lmoments",
type="GP", threshold=th)
# diagnostic plots
plot(pot_lmom)
rl_lmom <- [Link](pot_lmom, conf = 0.05, [Link]=
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 <- [Link]([Link](pot_mle, conf =
0.05,[Link]=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 <- [Link]([Link](pot_lmom, conf =
0.05,[Link]=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([Link](mle=[Link](rl_mle),
lmom=[Link](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')

[Link]=fevd(B1[,2], data=B1, threshold=[Link], type="GEV",method

="MLE", units="KVA",[Link]="seconds",period = "Seconds")
B1.fit1=fevd(B1[,2], data=B1, threshold=[Link],type="GP",method
="MLE", units="KVA",[Link]="seconds",period = "Seconds")
B1.fit2=fevd(B1[,2], data=B1, threshold=[Link],type="Gumbel", method
="MLE",units="KVA",[Link]="seconds",period = "Seconds")
B1.fit3=fevd(B1[,2], data=B1,
threshold=[Link],type="Exponential",method ="MLE",
units="KVA",[Link]="seconds",period = "Seconds")

[Link]=summary([Link], silent=TRUE)$AIC
[Link]=summary(B1.fit1, silent=TRUE)$AIC
[Link]=summary(B1.fit2, silent=TRUE)$AIC
[Link]=summary(B1.fit3, silent=TRUE)$AIC

[Link]
# [1] 39976258
[Link]
# [1] 466351.5
[Link]
# [1] 13934878
[Link]
# [1] 466330.8

plot([Link])
plot(B1.fit1)
plot(B1.fit2)
plot(B1.fit3)

rm(list=ls())
rm(list=ls())
rm(list=ls())

rm=[Link](file("clipboard"),header=T, sep="\t", [Link]=1)

getwd()
setwd("C:/Users/Lenovo/Documents")
[Link](rm, file="[Link]",[Link]=F)
rm=[Link](file("clipboard"),header=T, sep="\t")

str(rm)
attach(rm)

names(r)
reg=lm([Link]~X,data=r)
summary(reg)

plot(X,MAM, type="o",lty=9, xlab="Year", ylab="Precipitation(mm)")

library(naniar)
require("UpSetR")
gg_miss_upset(r)

library(mice)
[Link](r)
library(reshape2)
library(ggplot2
Abdi-Basid Analysis
abdi-basid@[Link]

#Import Data
data=[Link](file("clipboard"), header=T, sep="\t", dec=".")
str(data)
summary(data)
#Homogeneous
boxplot(data)
hist(data)

# 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, [Link]="period")
plot(fit)

# additional plots
monthplot(myts)
library(forecast)
seasonplot(myts)

# simple exponential - models level

# 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 <- [Link](myts)

#The forecast package provides functions for the automatic selection

#lag(ts, k) lagged version of time series, shifted back k

plot(Nile, ylab="Annual Flow of the river Nile")

abline(h= mean(Nile),col='blue')

[Link] <- efp(Nile ~ 1, type = "OLS-CUSUM")

[Link] <- efp(Nile ~ 1, type = "OLS-MOSUM")
[Link] <- efp(Nile ~ 1, type = "Rec-CUSUM")

opar <- par(mfrow=c(2,2))

plot([Link])
plot([Link])
plot([Link])
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")

[Link] <- breakpoints(Nile ~ 1)

[Link] <- breakfactor([Link], breaks = 1 )
[Link] <- lm(Nile ~ [Link] - 1)
plot([Link])

opar <- par(mfrow=c(2,1), mar=c(2,2,0,2))

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, [Link]=250)

polyfit <- function(i) x <- AIC(lm(y~poly(x,i)))

[Link](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=[Link](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@[Link]
#
#
#
# 0> >> << Importation Data
rm(list=ls())
djibclim=[Link](file("clipboard"),header=T, dec=".",sep="\t")
#File Change dir
#getwd()
#setwd("D:/[Link]/[Link]/[Link]/[Link]")
#save(djiclim, file = "[Link]")
#load(file = "[Link]")
#View(djibclim)
#str(djibclim);summary(djiclim)
#sapply(djibclim,sd)
#sapply(djibclim, mean, [Link]=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:/[Link]/[Link]/[Link] 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:\\[Link]\\[Link]\\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([Link])
[Link] <-
CraddockTest([Link],2,3)
plot([Link][,1],[Link],type
='l',xlab='',ylab='Craddock')
#
# Example of inhomogeneity
tmp <- [Link]
tmp[20:43,3] <- tmp[20:43,3]+1 # adding an error of +1 degree to
Milan's yearly average temperatures from 1980
[Link] <- CraddockTest(tmp,2,3)
[Link]()
plot([Link][,1],[Link],type
='l',xlab='',ylab='Craddock')
#
#Homogénéity Analysis
names(djiclim)
library(ggpubr)
ggqqplot(djiclim$Rain_Rrr)
#
library("car")
qqPlot(djiclim$Rain_Rrr)
#
#[Link](rnorm(100, mean = 5, sd = 3))
[Link](runif(100, min = 2, max = 4))
#require(stats)
[Link](djibclim$Rain_Rrr)
#
library(tseries)
[Link](djibclim$Rain_Rrr)
#
#Kolmogorov-Smirnov Tests
[Link](djibclim$Rain_Rrr)
[Link](djibclim$Rain_Rrr)
[Link](djibclim$Rain_Rrr)
#
#
# 2> >> << Program Rclim Dex
#library("RClimDex")
#[Link]()
getwd()
setwd("D:/[Link]/[Link]/[Link]
2/Climat/RclimDex")
source("[Link]")
require("[Link]")
require("Kendall")
require(tcltk)
#year month day without label
#I+1 i-1
#
#
# 3> >> << Program ClimPact2
#Download the above zip file to your computer
#([Link]
# 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('[Link].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('[Link].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:\\[Link]\\[Link]\\[Link] 2\\Climat\\
RClimPact\\climpact2-master")
source("[Link].r")
#
# 4> >> << Program RHtest
#setwd("D:\\[Link]\\[Link]\\RHtest")
#source ("RHtest_V4.r")
#source ("[Link]")
setwd("D:\\[Link]\\[Link]\\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.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the CDD index.
cdd <- [Link](ci)
#
## Create a climdexInput object from some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the cold spell duration index.
csdi <- [Link](ci)
#
## Create a climdexInput object from some data already loaded in and
## ready to go.
## Parse the dates into PCICt.
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
[Link] <- [Link]([Link](paste, [Link][,c("year",
"jday")]), format="%Y %j", cal="gregorian")
## Load the data in.
ci <- [Link]([Link]$MAX_TEMP,
[Link]$MIN_TEMP, [Link]$ONE_DAY_PRECIPITATION,
[Link], [Link], [Link], [Link]=c(1971, 2000))
## Create an annual timeseries of the CWD index.
cdd <- [Link](ci)

#
# 5> >> << ANCLIM

at=[Link](file("clipboard"),header=T, sep="\t")
write(at, file = "[Link]", 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=[Link](file("clipboard"),header=T, sep="\
t",dec=",",[Link]=1)
str(r)
attach(r)
names(r)
dim(r)

summary(r)

library(Hmisc)
describe(r,[Link]=c("mean","median","var","sd","valid.n"),horizont
al=TRUE)

sapply(r, sd)
sapply(r,range)

t=[Link](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)
[Link] = cor(t, method = c("pearson"))
corrplot([Link], method = "number", type = "lower")
#methode = number, ellipse, square, shade, color, pie,
#type = full" "upper" "lower"

library(RColorBrewer)
corrplot([Link], type = "upper", order = "hclust",col =
[Link](n = 8, name = "PuOr"))

palette = colorRampPalette(c("yellow", "orange", "red")) (20)

heatmap(x = [Link], col = palette, symm = TRUE)
#[Link]
[Link]
names(t)
#[Link]
s=[Link](file("clipboard"),header=T, sep="\
t",dec=",",[Link]=1)
t=[Link](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,[Link]=1.2,xlab="Ground
Based (mm/month)",ylab="Remotely sensed (mm/month)")
axis(1,[Link]=1.3);axis(2,[Link]=1.3)

par(new=TRUE)
plot(OB~CH,xlim=limx,ylim=limy,pch=20,xaxt="n",yaxt="n",col="red",xl
ab="",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",x
lab="",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="darkgre
en",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, [Link]=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=[Link](file("clipboard"),header=T, sep="\
t",dec=",",[Link]=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] "[Link]" "[Link]" "[Link]" "[Link]" "[Link]" "[Link]"
"[Link]"
[8] "[Link]" "[Link]" "[Link]" "[Link]" "[Link]" "[Link]"
"[Link]"
[15] "[Link]" "[Link]" "[Link]" "[Link]" "[Link]" "[Link]"
"[Link]"
[22] "[Link]" "[Link]" "[Link]"

##########
RSD = abs(sd(CH)/mean(CH));round(RSD, digits = 3)

###########################################
# display the diagram with the better model
#oldpar<-[Link](OB,CH,pch=19)
#[Link](OB,CH, text="Model 1")
s=[Link](file("clipboard"),header=T, sep="\
t",dec=",",[Link]=1)
attach(s)
names(s)
str(s)
cbind(sapply(s, sd, [Link]=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<-[Link](OB,CH,add=F, pch=19,[Link]=TRUE,
xlab="Standard deviation",ylab="Standard Deviation",main="",
[Link]=TRUE,ngamma=10,[Link]="green",[Link]=1,
[Link]=TRUE,[Link]="sample",[Link]=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,[Link]=1.1,normalize=F)

# now add the worse model

[Link](OB,PE,add=TRUE,col="black",
pcex=1.5,[Link]=1.1,normalize=F)
[Link](OB,TA,add=TRUE,col="pink",
pcex=1.5,[Link]=1.1,normalize=F)
[Link](OB,AR,add=TRUE,col="brown",
pcex=1.5,[Link]=1.1,normalize=F)
[Link](OB,ER,add=TRUE,col="green",
pcex=1.5,[Link]=1.1,normalize=F)

# get approximate legend position # add a legend

legend(1.3,1.6,legend=c("[Link]","[Link]","[Link]","[Link]","[Link]","JF.
ER"),horiz=FALSE,

pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","ora
nge","darkgreen"), cex=0.7)

# SAISONNALITY
#par(mfrow=c(1,2))
#"red","blue","brown","orange","darkgreen"
s=[Link](file("clipboard"),header=T, sep="\
t",dec=",",[Link]=1)
names(s)
attach(s)
cor(s)
str(s)

# show the "all correlation" display

require(plotrix)
[Link]([Link], [Link], col="blue",
[Link]=F,pcex=1.5,normalize=F)
[Link]([Link], [Link],add=TRUE, col="black",
pcex=1.5,normalize=F)
[Link]([Link], [Link],add=TRUE, col="pink",
pcex=1.5,normalize=F)
[Link]([Link], [Link],add=TRUE, col="brown",
pcex=1.5,normalize=F)
[Link]([Link], [Link],add=TRUE, col="green",
pcex=1.5,normalize=F)

legend(30,53
,legend=c("[Link]","[Link]","[Link]","[Link]
AT","[Link]","[Link]"),horiz=FALSE,
pch=c(15,19,19,19,19,19),col=c("darkgreen","red","blue","brown","ora
nge","darkgreen"), cex=0.7)

CHrmse=sqrt(sum(([Link])^2)/length([Link]))
PErmse=sqrt(sum(([Link])^2)/length([Link]))
TArmse=sqrt(sum(([Link])^2)/length([Link]))
ARrmse=sqrt(sum(([Link])^2)/length([Link]))
ERrmse=sqrt(sum(([Link])^2)/length([Link]))
round(rbind(CHrmse,PErmse,TArmse,ARrmse,ERrmse),3)

#RMS_Diff =
sum((([Link]-mean([Link]))-([Link]-mean([Link])))^2)/length([Link])
OBsd=sd([Link])
CHsd=sd([Link])
PEsd=sd([Link])
TAsd=sd([Link])
ARsd=sd([Link])
ERsd=sd([Link])
round(rbind(CHsd,PEsd,TAsd,ARsd,ERsd),3)

CHBIAS = sum([Link])/length([Link])
PEBIAS = sum([Link])/length([Link])
TABIAS = sum([Link])/length([Link])
ARBIAS = sum([Link])/length([Link])
ERBIAS = sum([Link])/length([Link])
round(rbind(CHBIAS,PEBIAS,TABIAS,ARBIAS,ERBIAS),3)

R=[Link](file("clipboard"),header=T, sep="\
t",dec=".",[Link]=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 <- [Link](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 = [Link](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")

[Link] <- rbind(mod1, mod2, mod3)

## basic Taylor plot

TaylorDiagram([Link], obs = "obs", mod = "mod", group = "model")

## Taylor plot by season

TaylorDiagram([Link], 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)
[Link]$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([Link], obs = "obs", mod = c("mod", "mod2"), group =
"model")

## End(Not run)
## Not run:
## all models, by season
TaylorDiagram([Link], 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([Link], 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)

[Link](r,r,add=FALSE,col="red",pch=4,[Link]=TRUE,xlab="MERR
A SD (Normalised)",ylab="RCA4 runs SD (normalised)",main="Taylor
Diagram",[Link]=TRUE,ngamma=3,[Link]=1,[Link]=TRUE,[Link]
nes=c(0.2,0.4,0.6,0.8,0.9),pcex=1,[Link]=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,[Link]=1.2,legend=c("volcano"),pch=4,col=c("r
ed"))
[Link](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<-[Link](ref,model1)
# now add the worse model
[Link](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
[Link](ref,model1,[Link]=FALSE)
[Link](ref,model2,add=TRUE,col="blue")

FIN DU PROGRAMME
Enjoy It !

Common questions