-
Notifications
You must be signed in to change notification settings - Fork 0
/
Script01
257 lines (215 loc) · 10.2 KB
/
Script01
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
########################################################################################################################################################################
# Script 01 - Time-series-specific analysis.
# Haubrock et al. The invasion curve and dynamics of a European-wide introduced species.
# This script includes the code for the following steps:
# (1) compute the monotonic trends for the abiotic variables over the surveyed period,
# (2) extract climatic data for the site,
# (3) compute the monotonic trends of mean annual temperature and total annual precipitation over the surveyed period,
# (4) compute the monotonic trend of abundance of Potamopyrgus antipodarum,
# (5) create meta regression analysis master file
# (6) export the results
#
########################################################################################################################################################################
rm(list=ls()) #clear the working environment
library(vegan)
library(R.utils)
library(raster)
library(ncdf)
library(rgdal)
library(rgeos)
library(lme4)
library(car)
library(ggplot2)
library(tidyverse)
library(RColorBrewer)
library(dplyr)
library(codyn)
library(reshape2)
library(metafor)
library(sp)
library(RColorBrewer)
library(glmulti)
(1) compute the monotonic trends for the abiotic variables over the surveyed period
#read data
NO3 <-read.csv("NITRO_NO3_mean.csv", sep=",", stringsAsFactors=FALSE)
#data formatting
NO3_data <- NO3[,-1]
rownames(NO3_data) <- NO3[,1]
CTLong <- data.frame(rows = rownames(NO3_data), stack(NO3_data))
CTLong <- CTLong[order(CTLong$rows), ]
CTLong$values<-as.numeric(CTLong$values)
NO3 <- split(CTLong$values, CTLong$rows)
NO3 <- NO3[lengths(NO3) >= 4] #select ony time series with length of4 or more occurrences
NO3_MK_acf<-as.data.frame(do.call(rbind,lapply(NO3[1:297],function(x)unlist(acf(x, na.action = na.pass))))) # check for temporal autocorrelation
NO3_MK_acf<-as.data.frame(do.call(rbind,lapply(NO3[1:297],function(x)unlist(pacf(x, na.action = na.pass))))) # check for temporal autocorrelation
NO3_MK<-as.data.frame(do.call(rbind,lapply(NO3[1:297],function(x)unlist(My.mmkh(x)))))# Modified Mann-Kendall Test For Serially Correlated Data Using Hamed and Rao (1998) Variance Correction Approach
NO3_MK <- cbind(rownames(NO3_MK), data.frame(NO3_MK, row.names=NULL)) # transform rownames into first column
colnames(NO3_MK)[1] = "site_id" # rename first column name
write.csv(NO3_MK, "NO3_MK.csv") #export results
# repeat this step for the other abiotic variables
# (2) Gather climatic data
# To run the following code, it is necessary to download the daily mean temperature and daily total precipitation data
# from the Gridded observational dataset for precipitation, temperature and sea level pressure in Europe,
# available at https://www.ecad.eu/download/ensembles/download.php
# 1) Download file "tg_ens_mean_0.1deg_reg_v22.0e.nc" from: https://www.ecad.eu/download/ensembles/download.php#datafiles
# 2) Extract .gz compressed files using:
# library(R.utils)
# R.utils::gunzip("tg_ens_mean_0.1deg_reg_v22.0e.nc")
# 3) Extract data, see:
# http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html
# https://stackoverflow.com/questions/20621200/extract-time-series-of-a-point-lon-lat-from-netcdf-in-r
install.packages("sp")
install.packages("raster")
install.packages("ncdf4")
library(sp)
library(raster)
library(ncdf4)
df<-read.csv("WGS84.csv",sep=",",stringsAsFactors=FALSE, fileEncoding="UTF-8-BOM") # file containing all the required coordinated
#Temperature
my_function <- function(lon, lat, site) {
ncin <- nc_open("tg_ens_mean_0.1deg_reg_v22.0e.nc")
print(ncin)
t <- ncvar_get(ncin,"time")
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(t)
obsoutput <- ncvar_get(ncin, # look for closest lon and lat
start = c(which.min(abs(ncin$dim$longitude$vals - lon)),
which.min(abs(ncin$dim$latitude$vals - lat)),
1),
count = c(1,1,-1))
DataMeanT <- data.frame(DateN = t, MeanDailyT = obsoutput)
nc_close(ncin)
Data=DataMeanT
Data$Date=as.Date(Data$DateN)
Data$Year=format(Data$Date,"%Y")
Data$Month=format(Data$Date,"%m")
Data$YearMonth=format(Data$Date, format="%Y-%b")
# Data <- within(DataMeanT, {
# Date <- as.Date(DateN, origin="2000-01-01")
# Year <- format(Data$Date,"%Y")
# Month <- format(Data$Date,"%m")
# YearMonth <- format(Data$Date, format="%Y-%b")
# })
Data_annual <- with(Data, aggregate(list("AirT" = MeanDailyT), list(Year=YearMonth),#Year or Month
FUN = mean, na.action = na.pass))
Data_annual$Site=site #This line adds a column with the Site name
write.table(Data_annual, file = "ClimaticData.csv", row.names=FALSE,
append = TRUE, col.names = F, sep = ", ", quote = TRUE) # Export table, with append=TRUE it adds coulmns to an existing table
return(Data_annual) # SAVE AGGREGATED DATA FRAME
}
# ITERATE THROUGH EACH LON/LAT PAIR ELEMENTWISE
write.table(data.frame("Year", "AirT", "Site"), file = "ClimaticData.csv", row.names=FALSE,
append = F, col.names = F, sep = ", ", quote = TRUE) # create a blank table template, values will be added to this template during the iteration
# ITERATE THROUGH EACH LON/LAT PAIR ELEMENTWISE
df_list <- Map(my_function, df$E, df$N, df$Site)
#formatting of data
Dv_temp <-read.csv("ClimaticData.csv", sep=",", stringsAsFactors=FALSE)
head(Dv_temp)
Dv_temp2<-do.call(rbind, str_split(Dv_temp$Year, '-'))
head(Dv_temp2)
Dv_temp2<-as.data.frame(Dv_temp2)
colnames(Dv_temp2)[1] <- "Year"
colnames(Dv_temp2)[2] <- "Month"
head(Dv_temp2)
Dv_temp$Year<-NULL
new_temp <- cbind(Dv_temp2, Dv_temp)
head(new_temp)
colnames(new_temp)[3] <- "AirT"
colnames(new_temp)[4] <- "Site"
new_temp$Year <- as.numeric(new_temp$Year)
new_temp$AirT <- as.numeric(new_temp$AirT)
head(new_temp)
new_temp2<-new_temp %>% group_by(Year,Site) %>% summarise(Mean_temp=mean(AirT))
head(new_temp2)
new_temp3 <- transform(new_temp2, LINK=paste(Site, Year, sep="_"))
head(new_temp3)
write.table(new_temp3,"Temp_data.csv",sep=";") # export result
####Precipitation
my_function <- function(lon, lat, site) {
# Mean daily precipitation
ncin <- nc_open("rr_ens_mean_0.1deg_reg_v22.0e.nc")
print(ncin)
r <- ncvar_get(ncin,"time")
tunits <- ncatt_get(ncin,"time","units")
nt <- dim(r)
obsoutput <- ncvar_get(ncin, # look for closest lon and lat
start = c(which.min(abs(ncin$dim$longitude$vals - lon)),
which.min(abs(ncin$dim$latitude$vals - lat)),
1),
count = c(1,1,-1))
DataMeanRR <- data.frame(DateN = r, MeanDailyRR = obsoutput)
nc_close(ncin)
Data=DataMeanRR
Data$Date=as.Date(Data$DateN,origin="1989-01-01")
Data$Year=format(Data$Date,"%Y")
Data$Month=format(Data$Date,"%m")
Data$YearMonth=format(Data$Date, format="%Y-%b")
# Data <- within(DataMeanT, {
# Date <- as.Date(DateN, origin="2000-01-01")
# Year <- format(Data$Date,"%Y")
# Month <- format(Data$Date,"%m")
# YearMonth <- format(Data$Date, format="%Y-%b")
# })
Data_annual <- with(Data, aggregate(list("Prec" = MeanDailyRR), list(Year=YearMonth),#Year or Month
FUN = mean, na.action = na.pass))
Data_annual$Site=site #This line adds a column with the Site name
write.table(Data_annual, file = "PrecData.csv", row.names=FALSE,
append = TRUE, col.names = F, sep = ", ", quote = TRUE) # Export table, with append=TRUE it adds coulmns to an existing table
return(Data_annual) # SAVE AGGREGATED DATA FRAME
}
# ITERATE THROUGH EACH LON/LAT PAIR ELEMENTWISE
write.table(data.frame("Year", "Prec", "Site"), file = "PrecData.csv", row.names=FALSE,
append = F, col.names = F, sep = ", ", quote = TRUE) # create a blank table template, values will be added to this template during the iteration
# ITERATE THROUGH EACH LON/LAT PAIR ELEMENTWISE
df_list <- Map(my_function, df$E, df$N, df$Site)
#formatting of data
Dv_prec <-read.csv("PrecData.csv", sep=",", stringsAsFactors=FALSE)
head(Dv_prec)
Dv2<-do.call(rbind, str_split(Dv_prec$Year, '-'))
head(Dv2)
Dv_2<-as.data.frame(Dv2)
colnames(Dv_2)[1] <- "Year"
colnames(Dv_2)[2] <- "Month"
head(Dv_2)
Dv_prec$Year<-NULL
new <- cbind(Dv_2, Dv_prec)
colnames(new)[3] <- "Prec"
colnames(new)[4] <- "Site"
new$Year <- as.numeric(new$Year)
new$Prec <- as.numeric(new$Prec)
new$Site <- as.character(new$Site)
str(new)
new2<-new %>% group_by(Year,Site) %>% summarise(Mean_prec=mean(Prec))
head(new2)
new3 <- transform(new2, LINK=paste(Site, Year, sep="_"))
head(new3)
write.table(new3,"Prec_data.csv",sep=";") # export result
# (3) compute the monotonic trends of mean annual temperature and total annual precipitation over the surveyed period,
# see step (1)
# (4) compute the monotonic trend of abundance of Potamopyrgus antipodarum,
#read data
Dv <-read.csv("Raw_data.csv", sep=",", stringsAsFactors=FALSE)
#data formatting
Dv_2 <- split(Dv$abundance, Dv$site_id)
Dv_2 <- Dv_2[lengths(Dv_2) >= 4]
# check for temporal autocorrelation
Dv_MK_acf<-as.data.frame(do.call(rbind,lapply(Dv[1:297],function(x)unlist(acf(x, na.action = na.pass)))))
Dv_MK_acf<-as.data.frame(do.call(rbind,lapply(Dv[1:297],function(x)unlist(pacf(x, na.action = na.pass)))))
#apply the modified MK test
MK<-as.data.frame(do.call(rbind,lapply(Dv_2[1:306],function(x)unlist(My.mmkh(x)))))
#data formatting
setDT(MK, keep.rownames = TRUE)[]
colnames(MK)[1] = "site_id"
colnames(MK)[6] = "P_value"
colnames(MK)[11] = "S_statistic"
MK$site_id<-as.numeric(MK$site_id)
write.csv(MK, "MK.csv") # export results
(5) create meta regression analysis master file
# All exported files were checked in excel for possible formatting issues
#read data
NO3_MK <-read.csv("NO3_MK.csv", sep=",", stringsAsFactors=FALSE)
colnames(O3_MK)[11] = "NO3"
NO3_MK_S<-NO3_MK[c("site_id","NO3")] #extract site_id and S-statistic, now renamed NO3
summary(NO3_MK_S) #inspect data
master_file <- merge(MK,NO3_MK_S,by="site_id") #merge the Dv trend analysis result with S-statistic of the MK trend test resuls
#This step is to be repeated with all other abiotic variables