.Rhistory

)
)
)
)
)
")"
")
save(c(starttime, endtime), file="OutputFolder/timerange_AMSdata.RData")
starttime
endtime
save(starttime, endtime, file="OutputFolder/timerange_AMSdata.RData")	# save starttime and endtime of dataset in an R Object
load("PWSQC Rnotebook/OutputFolder/timerange_AMSdata.RData")
rm(endtime)
load("PWSQC Rnotebook/OutputFolder/timerange_AMSdata.RData")
rm(list=(ls()))
library(geosphere) # to calculate distance bewteen lat/lon coordinates
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
neighbourlist <- vector("list", nrow(Meta))
for(i in 1:nrow(Meta)){
dist <- distm(cbind(Meta$lon, Meta$lat), c(Meta$lon[i], Meta$lat[i]), fun = distHaversine)	# make a list of distances to all stations including itself
neighbourlist[[i]] <- Meta$id[which(dist > 0 & dist <= range)] }	# select the ID's of stations where the distance is smaller than 'range' but larger than zero to avoid being matched with itself
save(neighbourlist, file=paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))	# save 'neighbourlist' as list in an R Object
rm(list=(ls()))
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
rm(list=(ls()))
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/FZ_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/HI_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
Ndataset2 <- Ndataset * defaultbiascorrection	# the multiplicationfactor does not impact SO in any way
Ndataset2[which((HI_flags == 1)|(FZ_flags == 1))] <- NA
rm(Ndataset) # delete Nmatrix from working memory can speed code up for large datasets.
for(i in 1:nrow(Meta)){
Nint <- Ndataset2[,i]
if((length(neighbourlist[[i]]) < nstat)|(length(which(is.na(Nint)==F)) < 1)){SOflag <- rep(-1, times=length(Nint))	# if there are not enough stations nearby or no observations in 'Nint', all intervals get flagged as -1
}else{
Nintrain <- rep(0, length=length(Nint))
Nintrain[which(Nint > 0)] <- 1
Nintraincum <- cumsum(Nintrain)	# cumulative intervals with nonzero rainmeasurements
comparestartrowA <- match((Nintraincum-rainyint+1), Nintraincum)-1	# row from which the window should start to have at least 'rainyint' rainy intervals. match the first value that belongs to a number of rainy intervals equal to 'rainyint'-1 (double minus becomes +), and detract 1 after. This makes sure that if multiple rows have same amount, the last one (nearest to matching interval) is chosen. Note that this may result in a rownumber of 0 (1 minus 1), which should be replaced by NA (see next line), as the 0th row does not exist.
comparestartrowA[which(comparestartrowA == 0)] <- NA
comparestartrowB <- c(rep(NA, times=compareint-1), 1:(length(Nint)-compareint+1))	# row from which the window should start to have 'compareint' number of intervals
comparestartrow <- ifelse(is.na(comparestartrowB), NA, ifelse((comparestartrowA < comparestartrowB), comparestartrowA, comparestartrowB))	# choose either 'compareint' steps before, or where 'rainyint' was reached
NeighbourVal <- Ndataset2[,which(Meta$id %in% neighbourlist[[i]])]
NeighbourVal[which(is.na(Nint)),] <- NA		# replace rows where 'Nint' is NA with NA, this is needed for checking 'matchint' later
cortable <- biastable <-  matrix(NA, ncol=ncol(NeighbourVal), nrow=nrow(NeighbourVal)) # table of size 'NeighbourVal' to fill in the correlation values
for(t in 1:length(Nint)){
print(paste("SO filter: station", i, "out of", nrow(Meta), " -  t",t, "out of", length(Nint)))
if(is.na(comparestartrow[t])){next}
NeighbourValselec <- NeighbourVal[comparestartrow[t]:t,]
columnselec <- which(apply(NeighbourValselec, 2, function(x) length(which(is.na(x)==F))) > matchint)	# find columns that have at least 'matchint' overlapping intervals
if(length(columnselec) < nstat){next}
cortable[t,columnselec] <- apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs')) # determine the correlations with all neighbour stations over the past period that was chosen, this can yield a NA when the comparing stations measures zeroes only
biastable[t,columnselec] <- apply(NeighbourValselec[,columnselec], 2, function(x) mean(Nint[comparestartrow[t]:t]/defaultbiascorrection - x, na.rm=T)/mean(x, na.rm=T) ) # calculate the relative bias in the mean of the raw observations with all biascorrected neighbour stations over the past period between 'comparestartrow[t]' and 't'. To do so, 'Nint', which is based on 'Ndataset' times 'defaultbiascorrection', should be divided by 'defaultbiascorrection'.
} # end of r-loop
SOflag <- rep(0, times=length(Nint))
SOflag[which(apply(cortable, 1, function(x) median(x, na.rm=T)) < corthres)] <- 1	# rows where the median of all neighbour correlations was below 'corthres' are flagged as Station Outlier
SOflag[which(apply(cortable, 1, function(x) length(which(is.na(x)==F))) < nstat)] <- -1	# 'cortable' will have more or equal amounts of NA values than 'biastable'. Where 'SOflag' is 0, enough stations were included in the calculation to consider that row in 'biastable'.
} # end of ifelse loop
if(exists("SO_flags")==F){ SO_flags <- SOflag
}else{ SO_flags <- cbind(SO_flags, SOflag) }
biascorrectiontimeline <- rep(defaultbiascorrection, times=length(Nint)) # start out with the 'defaultbiascorrection' for all stations at all time intervals
if(length(which(SOflag == 0)) > 0){
biasmed <- apply(biastable, 1, function(x) median(x, na.rm=T))	# the median of the bias between raw timeserie 'Nint' with neighbouring stations that are multiplied with the 'defaultbiascorrection'
for(brow in which(SOflag == 0)){
biasprev <- biascorrectiontimeline[brow]
biasnew <- 1 / (1+biasmed[brow])
if( abs(log(biasnew / biasprev)) > log(1+biasthres) ){	# means that if [1/(1+biasthres) > BCFnew/BCFprev > 1+biasthres], change it for the remainder of the timeline
biascorrectiontimeline[(brow+1):length(biascorrectiontimeline)] <- biasnew }
} # end of brow-loop
} # end of if-loop
if(exists("biascorrectiontable")==F){ biascorrectiontable <- biascorrectiontimeline
}else{ biascorrectiontable <- cbind(biascorrectiontable, biascorrectiontimeline) }
} # end of i loop
rm(list=(ls()))
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/FZ_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/HI_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
Ndataset2 <- Ndataset * defaultbiascorrection	# the multiplicationfactor does not impact SO in any way
#Ndataset2[which((HI_flags == 1)|(FZ_flags == 1))] <- NA
rm(Ndataset) # delete Nmatrix from working memory can speed code up for large datasets.
for(i in 1:nrow(Meta)){
Nint <- Ndataset2[,i]
if((length(neighbourlist[[i]]) < nstat)|(length(which(is.na(Nint)==F)) < 1)){SOflag <- rep(-1, times=length(Nint))	# if there are not enough stations nearby or no observations in 'Nint', all intervals get flagged as -1
}else{
Nintrain <- rep(0, length=length(Nint))
Nintrain[which(Nint > 0)] <- 1
Nintraincum <- cumsum(Nintrain)	# cumulative intervals with nonzero rainmeasurements
comparestartrowA <- match((Nintraincum-rainyint+1), Nintraincum)-1	# row from which the window should start to have at least 'rainyint' rainy intervals. match the first value that belongs to a number of rainy intervals equal to 'rainyint'-1 (double minus becomes +), and detract 1 after. This makes sure that if multiple rows have same amount, the last one (nearest to matching interval) is chosen. Note that this may result in a rownumber of 0 (1 minus 1), which should be replaced by NA (see next line), as the 0th row does not exist.
comparestartrowA[which(comparestartrowA == 0)] <- NA
comparestartrowB <- c(rep(NA, times=compareint-1), 1:(length(Nint)-compareint+1))	# row from which the window should start to have 'compareint' number of intervals
comparestartrow <- ifelse(is.na(comparestartrowB), NA, ifelse((comparestartrowA < comparestartrowB), comparestartrowA, comparestartrowB))	# choose either 'compareint' steps before, or where 'rainyint' was reached
NeighbourVal <- Ndataset2[,which(Meta$id %in% neighbourlist[[i]])]
NeighbourVal[which(is.na(Nint)),] <- NA		# replace rows where 'Nint' is NA with NA, this is needed for checking 'matchint' later
cortable <- biastable <-  matrix(NA, ncol=ncol(NeighbourVal), nrow=nrow(NeighbourVal)) # table of size 'NeighbourVal' to fill in the correlation values
for(t in 1:length(Nint)){
print(paste("SO filter: station", i, "out of", nrow(Meta), " -  t",t, "out of", length(Nint)))
if(is.na(comparestartrow[t])){next}
NeighbourValselec <- NeighbourVal[comparestartrow[t]:t,]
columnselec <- which(apply(NeighbourValselec, 2, function(x) length(which(is.na(x)==F))) > matchint)	# find columns that have at least 'matchint' overlapping intervals
if(length(columnselec) < nstat){next}
cortable[t,columnselec] <- apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs')) # determine the correlations with all neighbour stations over the past period that was chosen, this can yield a NA when the comparing stations measures zeroes only
biastable[t,columnselec] <- apply(NeighbourValselec[,columnselec], 2, function(x) mean(Nint[comparestartrow[t]:t]/defaultbiascorrection - x, na.rm=T)/mean(x, na.rm=T) ) # calculate the relative bias in the mean of the raw observations with all biascorrected neighbour stations over the past period between 'comparestartrow[t]' and 't'. To do so, 'Nint', which is based on 'Ndataset' times 'defaultbiascorrection', should be divided by 'defaultbiascorrection'.
} # end of r-loop
SOflag <- rep(0, times=length(Nint))
SOflag[which(apply(cortable, 1, function(x) median(x, na.rm=T)) < corthres)] <- 1	# rows where the median of all neighbour correlations was below 'corthres' are flagged as Station Outlier
SOflag[which(apply(cortable, 1, function(x) length(which(is.na(x)==F))) < nstat)] <- -1	# 'cortable' will have more or equal amounts of NA values than 'biastable'. Where 'SOflag' is 0, enough stations were included in the calculation to consider that row in 'biastable'.
} # end of ifelse loop
if(exists("SO_flags")==F){ SO_flags <- SOflag
}else{ SO_flags <- cbind(SO_flags, SOflag) }
biascorrectiontimeline <- rep(defaultbiascorrection, times=length(Nint)) # start out with the 'defaultbiascorrection' for all stations at all time intervals
if(length(which(SOflag == 0)) > 0){
biasmed <- apply(biastable, 1, function(x) median(x, na.rm=T))	# the median of the bias between raw timeserie 'Nint' with neighbouring stations that are multiplied with the 'defaultbiascorrection'
for(brow in which(SOflag == 0)){
biasprev <- biascorrectiontimeline[brow]
biasnew <- 1 / (1+biasmed[brow])
if( abs(log(biasnew / biasprev)) > log(1+biasthres) ){	# means that if [1/(1+biasthres) > BCFnew/BCFprev > 1+biasthres], change it for the remainder of the timeline
biascorrectiontimeline[(brow+1):length(biascorrectiontimeline)] <- biasnew }
} # end of brow-loop
} # end of if-loop
if(exists("biascorrectiontable")==F){ biascorrectiontable <- biascorrectiontimeline
}else{ biascorrectiontable <- cbind(biascorrectiontable, biascorrectiontimeline) }
} # end of i loop
cortable
NeighbourVal
NeighbourValselec
length(columnselec) < nstat
cortable
cortable[t, columnselec]
apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs'))
class(apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs')))
length(apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs')))
length(columnselec)
length(apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs')))
apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs'))
apply(NeighbourValselec[,columnselec[-4]], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs'))
apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs'))
cortable[t,columnselec]
cortable[t,columnselec] <- apply(NeighbourValselec[,columnselec], 2, function(x) cor(x, Nint[comparestartrow[t]:t], use='complete.obs'))
cortable[t,columnselec]
apply(NeighbourValselec[,columnselec], 2, function(x) mean(Nint[comparestartrow[t]:t]/defaultbiascorrection - x, na.rm=T)/mean(x, na.rm=T) )
rm(list=(ls()))
url <- "https://data.4tu.nl/ndownloader/files/24056603"
temp <- tempfile()  # make temporary folder to download zip-file
download.file(url, temp, mode="wb")
AMS_dataset <- read.csv(unz(temp, "AMS_dataset.csv"))   # read AMS_dataset from downloaded zip-file
AMS_metadata <- read.csv(unz(temp, "AMS_metadata.csv"))   # read AMS_metadata from downloaded zip-file
unlink(temp)
rm(url, temp)
rm(list=(ls()))
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
load("OutputFolder/timerange_AMSdata.RData")
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/FZ_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/HI_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/SO_flags_Filtersettings",Filtersettings,".RData"))
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
Time <- seq(starttime, endtime, by= "5 min")
Time
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
Nselec <- Ndataset[,StationSelec]
Nselec <- Ndataset[,StationSelec]
NselecAccum <- apply(Nselec, 1, function(x) cumsum(ifelse(is.na(raincolumn), 0, raincolumn)) + raincolumn*0)	# count the number of neighbours with measurements at each interval
NselecAccum <- apply(Nselec, 1, function(x) cumsum(ifelse(is.na(x), 0, x)) + x*0)	# count the number of neighbours with measurements at each interval
NselecAccum
dim(NselecAccum)
NselecAccum[,1]
NselecAccum[1,]
plot(NselecAccum[1,])
NselecAccum <- apply(Nselec, 1, function(x) cumsum(ifelse(is.na(x), 0, x)) + x*0)	# count the number of neighbours with measurements at each interval
plot(NselecAccum[1,])
s <- 1
s <- 2
StationSelec[s]
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
plot(N_ac)
s
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
ymax <- 1600 # mm upper limit of the graph
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
N_ac
points(Time, N_ac, col=rainbow(length(StationSelec))[s])
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
s
rainbow(length(StationSelec))[s]
points(Time, N_ac, col=rainbow(length(StationSelec))[s])
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
dev.new()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
points(Time[100], 500)
graphics.off()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
dev.off()
dev.off()
dev.off()
graphics.off()
dev.off()
dev.new()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
getwd()
setwd("C:/Users/vosl/OPENSENSE_sandbox/PWSQC Rnotebook")
source("InputFiles/Filtersettings.txt")	# obtain 'range', 'nstat', 'nint', 'HIthresA', 'HIthresB', 'compareint', 'rainyint', 'matchint', 'corthres' and 'Filtersettings'
load("OutputFolder/timerange_AMSdata.RData")
Meta <- read.table("InputFiles/metatableAms.txt", header=T, sep=",")
load("OutputFolder/Ndataset.RData")
load(paste0("OutputFolder/FZ_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/HI_flags_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/neighbourlist_Filtersettings",Filtersettings,".RData"))
load(paste0("OutputFolder/SO_flags_Filtersettings",Filtersettings,".RData"))
Time <- seq(starttime, endtime, by= "5 min")
FilterSetting <- "Strict" # can be Flex: intervals with Flag = 0 and Flag = -1, or Strict: only intervals with Flag = 0.
ymax <- 1600 # mm upper limit of the graph
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
s <- 2
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
dev.new()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
dev.new()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
}
ymax <- 1800 # mm upper limit of the graph
legend("topleft", legend=StationSelec, pch=1, col=rainbow(length(StationSelec))[1:length(StationSelec)])
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)")
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
lines(Time, N_ac, col=rainbow(length(StationSelec))[s])
}
legend("topleft", legend=StationSelec, pch=1, col=rainbow(length(StationSelec))[1:length(StationSelec)])
```{r}
FilterSetting <- "Strict" # can be Flex: intervals with Flag = 0 and Flag = -1, or Strict: only intervals with Flag = 0.
FilterSetting <- "Strict" # can be Flex: intervals with Flag = 0 and Flag = -1, or Strict: only intervals with Flag = 0.
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
lines(Time, N_ac, col=rainbow(length(StationSelec))[s]) # makes cumulative plots even if values are missing.
}
legend("topleft", legend=StationSelec, pch=1, col=rainbow(length(StationSelec))[1:length(StationSelec)])
rainaccum <- cumsum(ifelse(is.na(raincolumn), 0, raincolumn)) + raincolumn*0
FilterSetting <- "Strict" # can be Flex: intervals with Flag = 0 and Flag = -1, or Strict: only intervals with Flag = 0.
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset[,StationSelec[s]]), 0, Ndataset[,StationSelec[s]])) + Ndataset[,StationSelec[s]]*0
lines(Time, N_ac, col=rainbow(length(StationSelec))[s]) # makes cumulative plots even if values are missing.
}
legend("topleft", legend=StationSelec, lty=1, col=rainbow(length(StationSelec))[1:length(StationSelec)])
Ndataset2 <- Ndataset
ymax <- 1800 # mm upper limit of the graph
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
Ndataset2 <- Ndataset*biascorrectiontable
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset2[,StationSelec[s]]), 0, Ndataset2[,StationSelec[s]])) + Ndataset2[,StationSelec[s]]*0
lines(Time, N_ac, col="darkgrey") # makes cumulative plots even if values are missing.
points(Time[which(FZ_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(FZ_flags[,StationSelec[s]]==1)], "red") # plot intervals with FZ flag
points(Time[which(HI_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(HI_flags[,StationSelec[s]]==1)], "orange") # plot intervals with HI flag
points(Time[which(SO_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(SO_flags[,StationSelec[s]]==1)], "green") # plot intervals with HI flag
}
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
#Ndataset2 <- Ndataset*biascorrectiontable
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset2[,StationSelec[s]]), 0, Ndataset2[,StationSelec[s]])) + Ndataset2[,StationSelec[s]]*0
lines(Time, N_ac, col="darkgrey") # makes cumulative plots even if values are missing.
points(Time[which(FZ_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(FZ_flags[,StationSelec[s]]==1)], "red") # plot intervals with FZ flag
points(Time[which(HI_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(HI_flags[,StationSelec[s]]==1)], "orange") # plot intervals with HI flag
points(Time[which(SO_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(SO_flags[,StationSelec[s]]==1)], "green") # plot intervals with HI flag
}
dev.off()
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
#Ndataset2 <- Ndataset*biascorrectiontable
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset2[,StationSelec[s]]), 0, Ndataset2[,StationSelec[s]])) + Ndataset2[,StationSelec[s]]*0
lines(Time, N_ac, col="darkgrey") # makes cumulative plots even if values are missing.
points(Time[which(FZ_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(FZ_flags[,StationSelec[s]]==1)], "red") # plot intervals with FZ flag
points(Time[which(HI_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(HI_flags[,StationSelec[s]]==1)], "orange") # plot intervals with HI flag
points(Time[which(SO_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(SO_flags[,StationSelec[s]]==1)], "green") # plot intervals with HI flag
}
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
#Ndataset2 <- Ndataset*biascorrectiontable
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset2[,StationSelec[s]]), 0, Ndataset2[,StationSelec[s]])) + Ndataset2[,StationSelec[s]]*0
lines(Time, N_ac, col="darkgrey") # makes cumulative plots even if values are missing.
points(Time[which(FZ_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(FZ_flags[,StationSelec[s]]==1)], "red") # plot intervals with FZ flag
points(Time[which(HI_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(HI_flags[,StationSelec[s]]==1)], "orange") # plot intervals with HI flag
points(Time[which(SO_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(SO_flags[,StationSelec[s]]==1)], "green") # plot intervals with HI flag
}
ymax <- 1800 # mm upper limit of the graph
#StationSelec <- c("...")  # fill in the id's in Meta that you want plotted. Must be at least one but can be several.
# example: #
StationSelec <- c(1, 30, 44, 100, 120, 134)
plot(Time, rep(NA, times=length(Time)), ylim=c(0,ymax), ylab="Cummulative rainfall (mm)", main="Raw observations")
for(s in 1:length(StationSelec)){
N_ac <-  cumsum(ifelse(is.na(Ndataset2[,StationSelec[s]]), 0, Ndataset2[,StationSelec[s]])) + Ndataset2[,StationSelec[s]]*0
lines(Time, N_ac, col="darkgrey") # makes cumulative plots even if values are missing.
points(Time[which(FZ_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(FZ_flags[,StationSelec[s]]==1)], "red") # plot intervals with FZ flag
points(Time[which(HI_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(HI_flags[,StationSelec[s]]==1)], "orange") # plot intervals with HI flag
points(Time[which(SO_flags[,StationSelec[s]]==1)], Ndataset2[,StationSelec[s]][which(SO_flags[,StationSelec[s]]==1)], "green") # plot intervals with HI flag
}
graphics.off()
FZ_flags
plotstart <- strptime("2016-08-20 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
plotstart
plotend <- strptime("2016-08-21 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
(plotstart < starttime | plotend > endtime)
plotstart <- strptime("2015-08-20 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
plotend <- strptime("2016-08-21 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
if(plotstart < starttime | plotend > endtime)
if(plotstart < starttime | plotend > endtime){stop("The selected plot period is out of bounds.")}
plotstart <- strptime("2016-08-20 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
plotend <- strptime("2016-08-21 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
if(plotstart < starttime | plotend > endtime){stop("The selected plot period is out of bounds.")}
dim(Ndataset)
Ntimeselec <- Ndataset2[which(Time >= plotstart & Time <= plotend),]
colSums(Ntimeselec, na.rm=T)
apply(Ntimeselec, 1, function(x) length(which(is.na(x)==F)))
apply(Ntimeselec, 2, function(x) length(which(is.na(x)==F)))
apply(Ntimeselec, 2, function(x) length(which(is.na(x)==F))/length(x))	# counts number of
Raintimeselec <- colSums(Ntimeselec, na.rm=T)
Raintimeselec
DataAvailabilityThreshold <- 0.8
colorRampPalette(c("black", "white"))
Raintimeselec
(all(is.na(Raintimeselec)))
Raintimeselec
10 * Raintimeselec/max(Raintimeselec)
max(10 * Raintimeselec/max(Raintimeselec))
round(100 * Raintimeselec/max(Raintimeselec))
round(100 * (Raintimeselec/max(Raintimeselec)+1) # rewrite each cumulative value to the shade closest to it's 1/100 value.
round(100 * (Raintimeselec/max(Raintimeselec)+1)) # rewrite each cumulative value to the shade closest to it's 1/100 value.
round(100 * (1+ Raintimeselec/max(Raintimeselec))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
round(100 * (Raintimeselec/max(Raintimeselec)))+1 # rewrite each cumulative value to the shade closest to it's 1/100 value and add one so no rain also gets nonzero value.
plotstart <- strptime("2016-08-20 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
plotend <- strptime("2016-08-21 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
if(plotstart < starttime | plotend > endtime){stop("The selected plot period is out of bounds.")}
Ndataset3 <- Ndataset*biascorrectiontable
DataAvailabilityThreshold <- 0.8 # needs to have at least 80 % of data availability for a cumulative sum to be calculated.
Raintimeselec <- colSums(Ntimeselec, na.rm=T)
fracdata <- apply(Ntimeselec, 2, function(x) length(which(is.na(x)==F))/length(x))	# counts fraction of time intervals with non-NA observations.
Raintimeselec[which(fracdata < DataAvailabilityThreshold)] <- NA
if(all(is.na(Raintimeselec))){stop("No station with sufficient data in selected plot period.")}
# make a color indication based on the cumulative amounts
colfunc <- colorRampPalette(c("white", "purple"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec)))+1 # rewrite each cumulative value to the shade closest to it's 1/100 value and add one so no rain also gets nonzero value.
plot(Meta$lon, Meta$lat, col=colfunc(101)[colval], pch=20, cex=2)
colval
Ndataset3 <- Ndataset
Ntimeselec <- Ndataset3[which(Time >= plotstart & Time <= plotend),]
DataAvailabilityThreshold <- 0.8 # needs to have at least 80 % of data availability for a cumulative sum to be calculated.
Raintimeselec <- colSums(Ntimeselec, na.rm=T)
fracdata <- apply(Ntimeselec, 2, function(x) length(which(is.na(x)==F))/length(x))	# counts fraction of time intervals with non-NA observations.
Raintimeselec[which(fracdata < DataAvailabilityThreshold)] <- NA
if(all(is.na(Raintimeselec))){stop("No station with sufficient data in selected plot period.")}
# make a color indication based on the cumulative amounts
colfunc <- colorRampPalette(c("white", "purple"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec)))+1 # rewrite each cumulative value to the shade closest to it's 1/100 value and add one so no rain also gets nonzero value.
plot(Meta$lon, Meta$lat, col=colfunc(101)[colval], pch=20, cex=2)
colfunc(101)[colval]
colfunc(101)
colcal
colval
round(100 * (Raintimeselec/max(Raintimeselec)))
fracdata
Raintimeselec
Raintimeselec/max(Raintimeselec, na.rm=T
)
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T)))+1 # rewrite each cumulative value to the shade closest to it's 1/100 value and add one so no rain also gets nonzero value.
colval
max(colval, na.rm=T)
colfunc(101)[colval]
plot(Meta$lon, Meta$lat, col=colfunc(101)[colval], pch=20, cex=2)
# make a color indication based on the cumulative amounts
colfunc <- colorRampPalette(c("blue", "green", "purple"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T)))+1 # rewrite each cumulative value to the shade closest to it's 1/100 value and add one so no rain also gets nonzero value.
plot(Meta$lon, Meta$lat, pch=1, cex=2)
points(Meta$lon, Meta$lat, col=colfunc(101)[colval], pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=1, cex=2)
points(Meta$lon, Meta$lat, col=colfunc(101)[colval], pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=1, cex=1.5)
plot(Meta$lon, Meta$lat, pch=1, cex=1.5, xlab="lon", ylab="lat")
100*max(Raintimeselec)
100*max(Raintimeselec, na.rm=T)
round(100*max(Raintimeselec, na.rm=T))
round(100*max(Raintimeselec, na.rm=T))*c((1, 26, 51, 76, 101)
round(100*max(Raintimeselec, na.rm=T))*c(1, 26, 51, 76, 101)
round(max(Raintimeselec, na.rm=T))*c(1, 26, 51, 76, 101)
round(max(Raintimeselec, na.rm=T)*c(1, 26, 51, 76, 101))
colfunc <- colorRampPalette(c("blue", "green", "purple"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "white"  # plot zero values as white circles.
plot(Meta$lon, Meta$lat, pch=20, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=21, bg="lightgrey", cex=1.5, xlab="lon", ylab="lat")
plot(Meta$lon, Meta$lat, pch=21, bg="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=21, bg="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, bg=colorlist, pch=21, cex=2)
plot(Meta$lon, Meta$lat, pch=20, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, bg=colorlist, pch=21, cex=2)
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "white"  # plot zero values as white circles.
plot(Meta$lon, Meta$lat, pch=20, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, bg=colorlist, pch=21, cex=2)
round(max(Raintimeselec, na.rm=T)*c(1, 26, 51, 76, 101))
colorlist
plot(Meta$lon, Meta$lat, pch=20, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=1.5, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1)
c("yellow", colfunc[100][c(25, 50, 75, 100)])
c("yellow", colfunc(100)[c(25, 50, 75, 100)])
colfunc <- colorRampPalette(c("blue", "red"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]))
plotstart <- strptime("2016-08-20 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
plotend <- strptime("2016-08-21 00:00:00", format="%Y-%m-%d %H:%M:%S", tz="GMT")
if(plotstart < starttime | plotend > endtime){stop("The selected plot period is out of bounds.")}
Ndataset3 <- Ndataset*biascorrectiontable
# make a color indication based on the cumulative amounts
colfunc <- colorRampPalette(c("blue", "red"))
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", legend=round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]))
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", legend=round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]), pch=20)
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", legend=round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]), pch=20, cex=2)
legend("topleft", legend=round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]), pch=20, cex.pt=2)
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", legend=paste(round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), "mm""), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]), pch=20)
paste(round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), "mm")
colval <- round(100 * (Raintimeselec/max(Raintimeselec, na.rm=T))) # rewrite each cumulative value to the shade closest to it's 1/100 value.
colorlist <- colfunc(100)[colval]
colorlist[which(colval == 0)] <- "yellow"  # plot zero values as yellow circles.
plot(Meta$lon, Meta$lat, pch=1, col="lightgrey", cex=2, xlab="lon", ylab="lat")
points(Meta$lon, Meta$lat, col=colorlist, pch=20, cex=2)
legend("topleft", legend=paste(round(seq(0, max(Raintimeselec, na.rm=T), length=5),digits=1), "mm"), col=c("yellow", colfunc(100)[c(25, 50, 75, 100)]), pch=20)