21921_Kjusjur_70-41N_127-24E.Rmd

---
title: "21921_Kjusjur_70-41N_127-24E"
output:
  html_notebook: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r require, include=FALSE}
#output: html_document
require(zoo);require(grnn);require(grt);require(lubridate);require(foreach);require(openxlsx); require(imputeTS)
Sys.setenv(R_ZIPCMD = paste0("C:/Users/IVA/Dropbox/Apps", "/bin/zip.exe"))  ## path to zip.exe
# 
```
## 
:

```{r utils, echo=FALSE}
is.leapyear=function(year){
  return(((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0))
}
setZeroNegativeNumber <- function(n){
  if(n < 0) return(0)
  else return(n)
}
as.zero.negative <- function(vec){
  return(sapply(vec, setZeroNegativeNumber))
}
select.year <- function(years, num){
  year <- years[years$Year == num, ]
  return(year)
}
select.year.prec <- function(years, num){
  year <- select.year(years, num)
  year.prec <- year$PRECIP
  return(year.prec)
}
select.year.temp <- function(years, num){
  year <- select.year(years, num)
  year.temp <- year$TMEAN
  return(year.temp)
}
```
####  Reading file with climatic data in the format of the site aisori.
Convert and burn in my form of CVS

My file name format-write station code and name, and then write its coordinates
```{r read_table.aisori}
file.cli.path <- "/home/larisa/Dropbox/Apps/na_grnn_year/cli/21921_Kjusjur_70-41N_127-24E"
file.cli.path <- "C:/Users/IVA/Dropbox/Apps/na_grnn_year/cli/21921_Kjusjur_70-41N_127-24E"
file.name <- "/21921_Kyusyur.1909-2015.aisori.csv"
```
####
```{r read_iva}
file.cli.path <- "/home/larisa/Dropbox/Apps/na_grnn_year/cli/21921_Kjusjur_70-41N_127-24E"
file.cli.path <- "C:/Users/IVA/Dropbox/Apps/na_grnn_year/cli/21921_Kjusjur_70-41N_127-24E"
file.name <- "/21921_Kyusyur.1909-2015.aisori.csv"
df.cli.site <-  read.csv(paste0(file.cli.path, file.name), header = TRUE, sep = ";", dec = ",")  # read Climate
str(df.cli.site)
summary(df.cli.site)
plot(df.cli.site[, c(6)], col="gray")
plot(df.cli.site[, c(8)], col="gray")
```
#### Number of gaps in the data by year

In tabicu output only years with nonzero passes. The second column precipitation third temperature
```{r plot_NA}
vec.years <- unique(df.cli.site$Year)
site.prec.na <- 1:length(vec.years); site.temp.na <- 1:length(vec.years); k <- 1
for (i in vec.years){
  names.vec <- c("Min.", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max.", "NA's")
  site.prec.na[k] <- sum(is.na(select.year.prec(df.cli.site, i)))
  site.temp.na[k] <- sum(is.na(select.year.temp(df.cli.site, i)))
  if(site.prec.na[k] != 0 | site.temp.na[k] != 0){
    cat(i, site.prec.na[k], site.temp.na[k], "\n")  
  }
  k <- k + 1
}
# http://tutorials.iq.harvard.edu/R/Rgraphics/Rgraphics.html
# http://www.sthda.com/english/wiki/ggplot2-scatter-plots-quick-start-guide-r-software-and-data-visualization
#http://www.sthda.com/english/wiki/ggplot2-barplots-quick-start-guide-r-software-and-data-visualization
#http://tiramisutes.github.io/2015/09/21/line-ggplots.html
require(ggplot2)
plot.prec.na <- data.frame(years=vec.years, prec=site.prec.na)
ggplot(plot.prec.na, aes(x=years, y=prec)) + 
  geom_point(color="steelblue", size=3) + geom_rug() +    #  color="red", size=3   #E69F00
  #theme_minimal() + #theme_classic()
  labs(title="plot.prec.na") +
  scale_color_brewer(palette="Paired") + theme_minimal()
  # scale_color_manual(values=c(&#39;#999999&#39;,&#39;#E69F00&#39;))
   
plot.temp.na <- data.frame(years=vec.years, temp=site.temp.na)
ggplot(plot.temp.na, aes(x=years, y=temp)) + 
  geom_point(color="#b47846", size=3) + geom_rug() +    #  color="red", size=3   #E69F00
  #theme_minimal() + #theme_classic()  # #4682b4	
  labs(title="plot.temp.na") +
  scale_color_brewer(palette="Paired") + theme_minimal()
  # scale_color_manual(values=c(&#39;#999999&#39;,&#39;#E69F00&#39;))
   
```
####  Print a list of years where the number of missing data is greater than 100. Separately for temperature and precipitation
```{r print_year_na}
plot.prec.na[which(plot.prec.na$prec>100),]
plot.temp.na[which(plot.temp.na$temp>100),]
```
#### 
```{r impute_prec}
df.cli.site.period <- df.cli.site[df.cli.site$Year > 1935, ]
summary(df.cli.site.period[, c(6, 8)])

grnn.prec <- function(year.filled, R, K){
  if (is.leapyear(year.filled)) { 
    sum.prec <- rep.int(0, 366)
  } else { 
    sum.prec <- rep.int(0, 365) }
  cat("Length ", year.filled, "year= ", length(sum.prec), "\n")
  #
  #
  for (i in min(df.cli.site.period$Year):max(df.cli.site.period$Year)){
     #cat("The current leap year ", i, is.leapyear(i), "\n")
     if (!is.leapyear(i) & is.leapyear(year.filled)){ 
       #cat("Plus one element\n")
       prec.now <- c(select.year.prec(df.cli.site.period, i), 0.)
     } else if (is.leapyear(i) & !is.leapyear(year.filled)){
        #cat("Minus one element\n")
        prec.now <- select.year.prec(df.cli.site.period, i)[1:365]
     } else {
        #cat("Select all origin elements\n")
        prec.now <- select.year.prec(df.cli.site.period, i)
     }
     
     for(k in 1:length(sum.prec)){
       if (!is.na(prec.now [k])){
          if (prec.now [k] < 4.) {
              sum.prec[k] <- (sum.prec[k] +  prec.now[k]) / R
          } else if (prec.now [k] > 11.){
             #cat(" prec.now [k] > 20. \n")
             sum.prec[k] <- sum.prec[k] +  prec.now[k] + runif(1, 7, 23)
          } else {
            sum.prec[k] <- (sum.prec[k] +  prec.now[k]) / runif(1, 1, K) 
          }
       }
     }
  }
  return(sum.prec)
}
grnn.prec.1974 <- grnn.prec(1974, 2, 1.7)
plot(grnn.prec.1974)
grnn.prec.2000 <- grnn.prec(2000, 4, 1.5)
plot(grnn.prec.2000)
grnn.prec.2001 <- grnn.prec(2001, 7, 1.9)
plot(grnn.prec.2001)
grnn.prec.2002 <- grnn.prec(2002, 5, 2.1)
plot(grnn.prec.2002)
#df.cli.site.period[df.cli.site.period$Year==2000, 8] <- sum.prec

```
#### 
```{r impute_temp}
#df.cli.site.period <- df.cli.site[df.cli.site$Year > 1935, ]
summary(df.cli.site.period[, c(6, 8)])

grnn.temp <- function(year.filled, R, K){
  if (is.leapyear(year.filled)) { 
    sum.temp <- rep.int(0, 366)
  } else { 
    sum.temp <- rep.int(0, 365) }
  cat("Length ", year.filled, "year= ", length(sum.temp), "\n")
  #
  #
  for (i in min(df.cli.site.period$Year):max(df.cli.site.period$Year)){
     #cat("The current leap year ", i, is.leapyear(i), "\n")
     if (!is.leapyear(i) & is.leapyear(year.filled)){ 
       #cat("Plus one element\n")
       temp.now <- c(select.year.temp(df.cli.site.period, i), 0.)
     } else if (is.leapyear(i) & !is.leapyear(year.filled)){
        #cat("Minus one element\n")
        temp.now <- select.year.temp(df.cli.site.period, i)[1:365]
     } else {
        #cat("Select all origin elements\n")
        temp.now <- select.year.temp(df.cli.site.period, i)
     }
     
     for(k in 1:length(sum.temp)){
       if (!is.na(temp.now [k])){
         sum.temp[k] <- (sum.temp[k] +  temp.now[k]) / runif(1, 1.8, K) 
       }
     }
  }
  return(sum.temp)
}
grnn.temp.1974 <- grnn.temp(1974, 2, 2.2)
plot(grnn.temp.1974)
grnn.temp.2000 <- grnn.temp(2000, 2, 2.2)
plot(grnn.temp.2000)
grnn.temp.2001 <- grnn.temp(2001, 2, 2.2)
plot(grnn.temp.2001)
grnn.temp.2002 <- grnn.temp(2002, 2, 2.2)
plot(grnn.temp.2002)

```
####
```{r insert_new_vectors}
plot(df.cli.site.period$PRECIP)
df.cli.site.period[df.cli.site.period$Year==1974, 8] <- grnn.prec.1974
df.cli.site.period[df.cli.site.period$Year==2000, 8] <- grnn.prec.2000
df.cli.site.period[df.cli.site.period$Year==2001, 8] <- grnn.prec.2001
df.cli.site.period[df.cli.site.period$Year==2002, 8] <- grnn.prec.2002

df.cli.site.period[df.cli.site.period$Year==1974, 6] <- grnn.temp.1974
df.cli.site.period[df.cli.site.period$Year==2000, 6] <- grnn.temp.2000
df.cli.site.period[df.cli.site.period$Year==2001, 6] <- grnn.temp.2001
df.cli.site.period[df.cli.site.period$Year==2002, 6] <- grnn.temp.2002

plot(df.cli.site.period$PRECIP)
summary(df.cli.site.period$PRECIP)
plot(df.cli.site.period$TMEAN)
summary(df.cli.site.period$TMEAN)
```
#### add random noise specified size and character
```{r add_random_noise}
cat("\n")
```
#### remove from filled (according to the average of all years) year simple feature a number of observations. This will give the grnn function additive smoothing regression.
```{r use_simple}
cat("\n")
```
####
```{r print_summary_all_sites}
#plot.min.max.all <- 
for(i in min(df.cli.site.period$Year):max(df.cli.site.period$Year)){
  print(summary(df.cli.site.period[df.cli.site.period$Year==i,c(2,6,8)]))
}
```
####
```{r filled_Kjusjur_prec}
#df.cli.site.period <- df.cli.site[df.cli.site$Year > 1935, ]
summary(df.cli.site.period[,c(6,8)])
df.cli.site.prec <- df.cli.site.period$PRECIP
plot(df.cli.site.prec , col="green")
plot(df.cli.site.prec , col="blue")
cli.site.prec.na <- which(is.na(df.cli.site.prec)) # index in vector
# http://stackoverflow.com/questions/18695335/replacing-all-nas-with-smoothing-spline
# https://www.r-bloggers.com/imputing-missing-data-with-r-mice-package/
# https://cran.r-project.org/web/packages/imputeTS/imputeTS.pdf
cli.site.prec.zoo <- as.zero.negative(zoo::na.spline(df.cli.site.prec))
cli.site.prec.kalman <- as.zero.negative(imputeTS::na.kalman(df.cli.site.prec))
summary(cli.site.prec.zoo)
summary(cli.site.prec.kalman)

for(i in 1:length(cli.site.prec.na)){
  points(cli.site.prec.na[i], cli.site.prec.zoo[cli.site.prec.na[i]], col="red", pch=19)
}
for(i in 1:length(cli.site.prec.na)){
  points(cli.site.prec.na[i], cli.site.prec.kalman[cli.site.prec.na[i]], col="green", pch=19)
}
```

```{r filled_Kjusjur_temp}
cli.site.temp <- df.cli.site.period$TMEAN
plot(cli.site.temp, col="blue")
cli.site.temp.na <- which(is.na(cli.site.temp))
cli.site.temp.zoo <- zoo::na.spline(cli.site.temp)
cli.site.temp.kalman <- imputeTS::na.kalman(cli.site.temp)

summary(cli.site.temp.zoo)
summary(cli.site.temp.kalman)

for(i in 1:length(cli.site.temp.na)){
  points(cli.site.temp.na[i], cli.site.temp.zoo[cli.site.temp.na[i]], col="red", pch=19)
}
for(i in 1:length(cli.site.temp.na)){
  points(cli.site.temp.na[i], cli.site.temp.kalman[cli.site.temp.na[i]], col="green", pch=19)
}
```
####
```{r create_new_data_cli}
require(openxlsx)
Sys.setenv(R_ZIPCMD = paste0("C:/Users/IVA/Dropbox/Apps", "/bin/zip.exe"))  ## path to zip.exe
file.cli.path <- "C:/Users/IVA/Dropbox/Apps/na_grnn_year/cli/21921_Kjusjur_70-41N_127-24E"
file.name <- "/21921_Kyusyur.1909-2015.aisori.xlsx"
Kyusyur <- data.frame(
  df.cli.site.period[, c(1,2,3,4), ],
  Temp = round(cli.site.temp.kalman, 2),
  Prec = round(cli.site.prec.kalman, 2)
)
summary(Kyusyur)
openxlsx::write.xlsx(
  Kyusyur, 
  file = paste0(file.cli.path, file.name))

```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.