ThesiR:

# install.packages('rgeos', type='source')
# install.packages('rgdal', type='source')
rm(list=ls())
require(rgeos)
require(sp)
require(maptools)
require(rgdal)
require(spatstat)
require(kriging)
require(tidyverse)
#devtools::load_all()
#"UTF-8","latin1"

path_estimation <- "estimations_predicts/estimations_v2_1.rds"
results <- readRDS(path_estimation)
full_join(x = data.frame(a=c(1,2)),y=data.frame(a=c(3,4)))
# ajustando o model lm para o log(altitude)
data("data_1")
dados <- data_1
alt_dados <- dados %>% filter(Estação==1) %>% select(Estações,Altitude,LATITUDE,LONGITUDE) %>% unique.data.frame()
model.alt <- lm(log(Altitude)~as.numeric(LONGITUDE)+as.numeric(LATITUDE),data = alt_dados)


am <- readOGR("AM",encoding ="UTF-8")
plot(am)
coordinates(am)
ID = cut(coordinates(am)[,1],quantile(coordinates(am)[,1]),include.lowest=T)
isTRUE(gpclibPermitStatus())
novo_mapa = unionSpatialPolygons(am,ID)
mapa_dis = gUnaryUnion(novo_mapa)
LONGITUDE <- mapa_dis@polygons[[1]]@Polygons[[1]]@coords[,1]
LATITUDE <- mapa_dis@polygons[[1]]@Polygons[[1]]@coords[,2]
pontos_interpolcao <- cbind(LONGITUDE,LATITUDE) %>% data.frame()
alt_est <- predict(model.alt,newdata=pontos_interpolcao) %>% data.frame()
poly <- list(data.frame(pontos_interpolcao))
x11()
plot(am)
points(pontos_interpolcao,cex=2,col=2)
text(pontos_interpolcao,cex=.9)
p = c(2210,2100,831:1292,675,1300:1700,2249:2304)
grupo <- NULL
dist_min <-distt<- NULL
u1 = "Nova Olinda do Norte";u2="Maraã";u3="Itamarati"
u = c(u1,u2,u3)
long.lat.u <- am@data %>% filter(NOME %in%u) %>% select(NOME,LATITUDE,LONGITUDE)
for(j in 1:nrow(pontos_interpolcao)){
  for(i in 1:nrow(long.lat.u)){
    dist_min[i] <- dist(rbind(long.lat.u[i,-1],pontos_interpolcao[j,]))
  }
  grupo[j]<- paste0("Grupo ",which.min(dist_min))
  distt <- rbind(distt,dist_min[which.min(dist_min)]%>%data.frame())

}
length(grupo)
nrow(distt)
omega <- function(x){
  tau_m = median(x)#.5#max(distt)
  C_n = exp((-0.5*(x)^2)/tau_m^2)
  omega = C_n/sum(C_n)
  return(omega)
}
TT=252
n <- length(LATITUDE)
om <-data.frame(grupo,distt) %>% rename(d=".")
head(om)
nrow(om)
om1 <- om %>% group_by(grupo) %>% summarise(peso=omega(d))
len.gr <- om %>% group_by(grupo) %>% summarise(n=length(grupo))%>% select(n)
om1 <- om1 %>% group_by(grupo) %>% summarise(peso=rep(peso,each=TT))

results <- readRDS("dados_resultados\\estimations.rds")
database <- results$data
y.min <- database %>% group_by(Grupos,Ano,Mês) %>% summarise(y.min=mean(UR))


res = NULL
for(j in 1:3){
  g<-paste0("Grupo ",j)
  aux1 <- y.min %>% filter(Grupos==g)
  aux2 <- apply(aux1,2,rep,len.gr$n[j])
  res <- rbind(res,aux2 %>% data.frame())
}


cov_kl <- data.frame(a=rep(1,n*TT),rep(alt_est$.,each=TT),rep(sin(1:12),21*n),rep(sin(1:12),21*n)) %>% as.matrix()
colnames(cov_kl)<-colnames(results$cov_kl)
cov_delta <- data.frame(a=rep(1,n*TT),rep(sin(1:12),21*n),rep(sin(1:12),21*n))
colnames(cov_delta)<-colnames(results$cov_delta)
ncx <- ncol(cov_kl)
ncv <- ncol(cov_delta)
grupos<-om1$grupo
w <- om1$peso
y<-res$y.min %>% as.numeric()
quantil<-results$quantil
grupos= as.factor(grupos);gr = levels((grupos))
DADOS = data.frame(y,cov_kl,cov_delta,grupos)
colnames(DADOS) = c("y",
                    colnames(cov_kl),
                    paste0(colnames(cov_delta),"_delta"),
                    "grupos")

EMVteta2 <- c(results$ar$Estimate,
              results$kl$Estimate,
              results$delta$Estimate,
              results$alpha$Estimate)

betas <- EMVteta2[-1][(1:ncx)]
lambda <- EMVteta2[-1][(ncx+1):(ncx+ncv)]
phi_ar <- EMVteta2[1]
alfa.est <- EMVteta2[-1][(ncx+ncv+1):length(EMVteta2[-1])]
#alfa.est <- c(.5,.1,.9)
link.kl <- results$link.kl
link.delta <- results$link.delta
TT=252
delta <- link.delta$linkinv(as.matrix(DADOS[c(paste0(colnames(cov_delta),"_delta"))]) %*% lambda)
id_eta_first <- seq(1,nrow(DADOS),TT)
id_eta_last <- seq(TT,nrow(DADOS),TT)
eta_1 <- cov_kl[,-1]%*%betas[-1]
eta_2 <- phi_ar*(link.kl$linkfun(DADOS$y) - eta_1)
eta_kl <- betas[1]+eta_1[-c(id_eta_first)] + eta_2[-c(id_eta_last)]
klq <- link.kl$linkinv(eta_kl)


est <- data.frame(Grupos=DADOS$grupos[-id_eta_first],
                  klq,
                  delta=delta[-id_eta_first],
                  UR=DADOS$y[-id_eta_first],
                  Mês=res$Mês[-id_eta_first],
                  Ano=res$Ano[-id_eta_first],
                  peso=om1$peso[-id_eta_first],
                  t=rep(1:252,n)[-id_eta_first]
                 ) %>%
#  unite("t",Mês:Ano,sep = "/") %>%
  mutate(bq_t=1/delta,
         phi_ts = log(klq)/log(1-(1-quantil)^(delta)),
         aq_t = 1/phi_ts) %>%
         #t=unite("t",Mês:Ano)),
         #UR=database$UR[-id_eta_first],
         #Mês=res$Mês[-id_eta_first],
         #Longitude=database$LONGITUDE[-id_eta_first]) %>%
  arrange(Grupos)
head(est,13)

#alfa.est <- c(.9,.9,.9)
rstable_pos <- function(n,alfa){
  U <- runif(n,0,pi)
  v <- rexp(n,1)
  a_U <- ((sin(alfa*U)/sin(U))^(1/(1-alfa)))*((sin((1-alfa)*U))/sin(alfa*U))
  E <- (a_U/v)^((1-alfa)/alfa)
  return(E)

}
#alfa.est = c(.1,.1,.1)
nn=500
B=1
estt <-matrix(0,nrow(est),B)
for(b in 1:B){
  j=1
  esti<-NULL
  while(j <= length(alfa.est)){
    #qest=NULL
    G <- paste("Grupo",j)
    w.aux <- w[[j]];L<-length(w.aux)
    alfa.aux <- alfa.est[j]
    #set.seed(10)
    e <- est %>% dplyr::filter(Grupos==G) %>%
      mutate(#peso=w.aux %>% rep(each=251),
             zm = rstable_pos(nn,alfa.aux) %>% mean(trim=.30),
             H=zm*peso,
             betaa=rbeta(1,H,1),
             y = (1-betaa^(delta))^(phi_ts),
             yqcond = (1-(1-quantil)^(1/(bq_t*H)))^(1/aq_t),
             yqmarg = (1-exp(-(delta/peso)*(-log(1-quantil))^(1/alfa.aux)))^phi_ts,
             #ysim= rkumar(n,aq_t,bq_t*zm*peso) %>% median()
      )

    esti=rbind(esti,e)
    j=j+1
  }
  #estt[,b]<-as.matrix(esti)
}



fig2 = esti %>% slice(1:251)%>%dplyr::group_split(Grupos) %>%
  purrr::map(
    ~ggplot(.)+geom_line(aes(x=t,y=y))+#geom_point(aes(x=t,y=UR,colour=factor(d)))+
      #geom_line(aes(x=t,y=q50),col=2)+
      #geom_line(aes(x=t,y=q50,colour="red"))+
     # geom_ribbon(aes(x=t,ymin=q25,ymax=q75),alpha = 0.5)+
      ylab("Relative humidity")+ #scale_y_continuous(labels=scales::percent)+
      #scale_x_discrete(labels = c(month.abb))+
    #  facet_grid(~Estações,scale="free")+
      theme_bw()+
      #geom_hline(data = d,aes(yintercept = kp))+
      #geom_hline(yintercept = q)+
      theme(legend.position = "none",axis.text.x=element_text(angle=90,size=8)))%>%
  cowplot::plot_grid(plotlist = .,nrow=1)
x11()
fig2



pontos.interpolcao.2 <- apply(pontos_interpolcao,2,rep,each=251)

esti$LATITUDE <- pontos.interpolcao.2[,2]
esti$LONGITUDE <- pontos.interpolcao.2[,1]
esti.1 <- esti %>% select(Ano,Mês,y,LATITUDE,LONGITUDE)
esti.2 <-  readRDS("dados_resultados\\valores_preditos.rds")
esti.2 <- esti.2 %>% select(Ano,Mês,y,LONGITUDE,LATITUDE)




ano <- unique(esti.1$Ano)[-1][1]
mes <- month.abb[1:2]

for(j in ano){
  for(k in mes){
    z1<-esti.2%>%filter(Ano==j & Mês==k) %>% mutate(LATITUDE=as.numeric(LATITUDE),LONGITUDE=as.numeric(LONGITUDE))
    z = esti.1 %>%filter(Ano==j & Mês==k) %>% mutate(LATITUDE=as.numeric(LATITUDE),
                                                              LONGITUDE=as.numeric(LONGITUDE)) %>% slice(c(p,5)) %>%
      full_join(y=z1)

    kriged <- kriging(z$LONGITUDE,z$LATITUDE,
                      z$y, polygons=poly, pixels=300,lags = 10)

    name.file <- paste0("dados_resultados\\",k,"_t",j,".rds")
    valores <- kriged$map
    valores$Ano <- rep(j,nrow(valores))
    valores$Mês <- rep(k,nrow(valores))

    saveRDS(valores,name.file)

  }
}

g# Analysis of results

## Packages
require(ggmap)
require(scales)
require(metR)
library(gganimate)
require(transformr)
years <- c(2000:2020)[-1][1]
month <- factor(month.abb,levels=c(month.abb[1:12]))[1:2]
data.base.pred <- NULL
for(j in years){
  for(k in month){
    name.file <- paste0("dados_resultados\\",k,"_t",j,".rds")
    data.base.aux <- readRDS(name.file)
    data.base.pred <- rbind(data.base.pred,data.base.aux)
  }
}
# coordinates(data.base.pred)<- ~x+y
# gridded(data.base.pred) = TRUE
#teste <- cbind(data.base.pred$x,data.base.pred$y) %>% as.polygonal()

library(RColorBrewer) # mix das cores
lm.palette <- colorRampPalette(c("lightblue","blue", "darkblue"), space = "rgb")
lm.palette<-colorRampPalette(c("#000099", "#00FEFF", "#45FE4F",
                   "#FCFF00", "#FF9400", "#FF3100"))

#lm.palette <- colorRampPalette(c("blues"), space = "rgb")

maps <- data.base.pred %>%mutate(Mês=factor(Mês,levels = c(month.abb)))  %>% ggplot(aes(x=x,y=y,z=pred))+
  metR::geom_contour_fill(color="white")+
  #geom_contour(color="white")+
  scale_x_continuous(expression(LONGITUDE),expand = c(0,0))+
  scale_y_continuous(expression(LATITUDE),expand = c(0,0))+
  # scale_fill_gradient((as.expression(paste("l(a,b;y)")))
  #                     ,low="black", high="gray",oob=squish,na.value = "white")

  scale_fill_gradientn("UR",
                      colours =lm.palette(30) %>% sort(decreasing = T) ,na.value = "white")+
  facet_wrap(.~Mês)+
  theme_bw()+
  theme(legend.key.height = unit(2,"cm"))
x11()
maps
# Video

library(RColorBrewer) # mix das cores
lm.palette <- colorRampPalette(c("lightblue","blue", "darkblue"), space = "rgb")
maps <- data.base.pred %>%mutate(Mês=factor(Mês,levels = c(month.abb)),
                                 Ano=as.character(Ano)) %>% arrange(Ano)%>%
  ggplot(aes(x=x,y=y,z=pred))+
  metR::geom_contour_fill()+
  geom_contour(color="white")+
  scale_x_continuous(expression(LONGITUDE),expand = c(0,0))+
  scale_y_continuous(expression(LATITUDE),expand = c(0,0))+
  # scale_fill_gradient((as.expression(paste("l(a,b;y)")))
  #                     ,low="black", high="gray",oob=squish,na.value = "white")

  scale_fill_gradientn("UR",
                       colours =lm.palette(20) ,na.value = "white")+
  facet_wrap(.~Mês)+
  theme_bw()+
  theme(legend.key.height = unit(2,"cm"))+
  transition_manual(Ano)+
  labs(title = "Year =  {current_frame}")

animate(maps,duration = 63,width = 800, height = 600,renderer = av_renderer()) #renderer = gifski_renderer() gif
anim_save("maps.mp4")



21*3