alan navarro

##################################################################
#### Importar una base de datos de temperatura maxima diaria
#### Fuente Estacion Observatorio Hermosillo CONAGUA
#### Tabla de Excel, primer columna son "year-month"
#### columnas 2-32 dias del mes
library(lubridate)
library(chron)
library(gdata)
library(hydroGOF)
library(hydroTSM)
library(zoo)
library(openair)
#### Hay que tener cuidado pues al parecer muchas funciones se 
#### enmascaran
### Fijamos el directorio de trabajo
setwd('   ')
###################################
urlfile<-'https://alanphd.com/archivos/temp_max_obs_2000_2017.csv'
### Leer base de datos
tmaxs<- read.csv(urlfile, 
                  header=T, na.strings= c('NA', " ", ""), 
                  quote="\"", 
                  dec=".",
                  stringsAsFactors = FALSE)
### usamos estos valores para verificar que se conserva la 
### estructura de los datos 
tmaxs[1,3] # NA
tmaxs[5,3] # 40.5
names(tmaxs) # "fecha" "X1" ... "X31"
dim(tmaxs) # 32 columnas
216/12 ### 216 hileras entre 12 meses = 18 years
### Valores INAP se usa como NA 
which(tmaxs == "INAP")
### Aseguramos que sean variables numericas
for (i in 2:32) {if (is.numeric(tmaxs[,i])) 
  {tmaxs[,i] <- tmaxs[,i]} else
  {tmaxs[,i]<-as.numeric(tmaxs[,i])}}
## Adicionamos dos variables al inicio de la tabla
tmaxs = cbind(year = 0, mes= 0, tmaxs); names(tmaxs)
### Llenamos las primeras dos columnas con la variable "fecha"
class(tmaxs$fecha) # character
tmaxs$year = as.numeric(substr(tmaxs$fecha, 1, 4))
tmaxs$mes = as.numeric(substr(tmaxs$fecha, 6, 7))
tmaxs$fecha <- NULL; head(tmaxs)
#####
month31<-c(1,3,5,7,8,10,12) # meses que terminan en 31
# First loop
############# let's get rid of those days w/o 31 days
for (i in 1:nrow(tmaxs)) 
  {if (any(tmaxs$mes[i] == month31)) 
    {tmaxs$X31[i] = tmaxs$X31[i]} else
   {tmaxs$X31[i] = 333} }
####################################################
tmaxs$X31
# Los dias 30 de Febrero
####################################################
for (i in 1:nrow(tmaxs)) 
  {if (tmaxs$mes[i] == 2) {tmaxs$X30[i] = 333}
  else {tmaxs$X30[i] = tmaxs$X30[i]}}
####################################################
### Years bisiestos
fechas<-min(tmaxs$year):max(tmaxs$year)
bis <- as.character(leap.year(fechas))
### rufles is just a dog's name 
rufles = data.frame(years=fechas,bisiesto=bis); rufles
rufles.s<-subset(rufles, rufles$bisiesto == "TRUE"); rufles.s
dim(rufles.s)
bisiestos.y<-rufles.s$years; bisiestos.y
# I am sure there must be a shorter way to do this; I enjoy R 
# I had fun ...
##############################################################
for (i in 1:nrow(tmaxs)) { 
  if ((tmaxs$mes[i] == 2) & (any(tmaxs$year[i] == bisiestos.y))) 
    {tmaxs$X29[i] = tmaxs$X29[i]}
  else {if (tmaxs$mes[i] == 2) {tmaxs$X29[i] = 333}
    else {tmaxs$X29[i] = tmaxs$X29[i]}}
}
##############################################################
head(tmaxs)
### usamos estos valores para verificar que se conserva la 
### estructura de los datos 
tmaxs[5,4] # 40.5
tmaxs[1,4] # NA
### por seguridad salvamos el archivo
getwd() # en nuestro directorio de trabajo
save(tmaxs,file="tmaxs.rda")
##############################################################
### EMPEZAMOS A INTEGRAR UNA SERIE DE TIEMPO
tmaxs_m <- data.matrix(tmaxs[ ,3:33]) # creamos una matrix
#Vectorizar la matriz función del paquete gdata
tmaxs_m<-unmatrix(tmaxs_m, byrow=T) 
#esta es una mejor forma de convertir una matrix en un vector
#Quitar los factores al vector
names(tmaxs_m) <- NULL
tmaxs_v = tmaxs_m
# Hacer vector de fechas
ini = min(tmaxs$year)
final = max(tmaxs$year)
ini_d <- as.Date(paste(ini,"01","01",sep="-")); ini_d
final_d <- as.Date (paste (final,"12","31",sep="-")); final_d
datesE <- dip(ini_d, final_d)
length(datesE) ####
#### numero total de dias 6,575
### en este caso son dias calendario
### Estimamos manualmente el numero de dias
abs(min(tmaxs$year)- max(tmaxs$year))
### 17 en realidad son 18 years desde 2000 al 2017
dias=(365*(abs(min(tmaxs$year)- max(tmaxs$year)) + 1))+(length(bisiestos.y))
if(length(datesE)==(dias)) {paste("correcto :-)")} else { paste("mal :-(")}
##############################
# Esta es la prueba de fuego, los valores con 333 deben de ser igual a la 
# diferencia de los dias calendario y los dias de la matrix vectorizada
length(tmaxs_v) #[1] 6696
length(datesE) #[1] 6575
daysoff<-which(tmaxs_v==333)
length(daysoff) # [1] 121
length(tmaxs_v) - length(datesE) #[1] 121
# ya estamos mas cerca de poder crear nuestro "zoo" object
# ahora hay que eliminar los valores 333, antes 
### hacemos una copia
tmaxs_v1 = tmaxs_v
#### eliminamos los dias que NO existen
tmaxs_v1 <-tmaxs_v1[-daysoff]
length(tmaxs_v1) #[1] 6575 
tmaxs_zoo<-vector2zoo(tmaxs_v1, datesE, date.fmt = "%Y-%m-%d")
tmaxs_zoo[160]
tmaxs[tmaxs$year==2000 & tmaxs$mes== 6,]
# saves as an R object
save(tmaxs_zoo,file="tmaxs_zoo.rda") 
#### tomado de http://oscarperpinan.github.io/R/zoo.pdf
Year <- function(x)as.numeric(format(x, "%Y"))
Year=Year(index(tmaxs_zoo))
### termina tomado de http://oscarperpinan.github.io/R/zoo.pdf
### Funcion para contar dias temp > a 
daysAbove_X = function (x){
  length (which (x > 44))
}
#### cambiar el valor de la función de arriba 38, 40, 42, 44
Gthan38 =aggregate(tmaxs_zoo, by=Year,FUN=daysAbove_X)
Gthan40 =aggregate(tmaxs_zoo, by=Year,FUN=daysAbove_X)
Gthan42 =aggregate(tmaxs_zoo, by=Year,FUN=daysAbove_X)
Gthan44 =aggregate(tmaxs_zoo, by=Year,FUN=daysAbove_X)
min_v= c(min(Gthan38), min(Gthan40), min(Gthan42), min(Gthan44))
max_v= c(max(Gthan38), max(Gthan40), max(Gthan42), max(Gthan44))
max(tmaxs_zoo, na.rm = T) ### 49.5 grados
### 120 dias temp > 38 grados C. (101-146)
### 78 dias temp > 40 grados C. (57-95)
### 39 dias temp > 42 grados C. (25-51)
### 10 dias temp > 44 grados C. (2-16)
Umbral=c(38,40,42,44)
NumD=c(mean(Gthan38), mean(Gthan40), mean(Gthan42), mean(Gthan44))
#### adaptada de 
# https://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars
####### Primera gráfica ##############################################
plot(Umbral,NumD,
     ylim=range(c(0, 150)),
     pch=19, 
     xlab="Umbral temperatura máxima (grados Celsius)", ylab="Días (rango)",
     main="Estación Observatorio Hermosillo (2000-2017) \n Días al año con temperatura máxima > Umbral",
     xaxt="n",
     frame.plot=F)
abline(h=seq(0,150,10), v=seq(38,44, 2), lty=3, col="gray")
axis(1, at = seq(38, 44, by = 2))
arrows(Umbral, min_v, Umbral, max_v, length=0.05, angle=90, code=3)
text(Umbral,NumD, labels=round(NumD,0), cex= 0.8, pos=4)
text(Umbral,min_v, labels=round(min_v,0), cex= 0.7, pos=1, col="red")
text(Umbral,max_v, labels=round(max_v,0), cex= 0.7, pos=3, col="red")
###################################################################
#### De objeto zoo lo regresamos a dataframe
tmaxs_df= fortify.zoo(tmaxs_zoo)
head(tmaxs_df)
class(tmaxs_df$Index)
tmaxs_df$date <- as.POSIXct(tmaxs_df$Index)
#### Cambiar al idioma ya que mi Mac está en Inglés
# trace(calendarPlot, edit=TRUE)
# weekday.abb <- substr(c("Domingo", "Lunes", "Martes", "Miercoles", "Jueves", "Viernes", "Sabado"), 1, 1)[((6:12) + w.shift)%%7 + 1]
Sys.getlocale() 
# "en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8"
Sys.setlocale("LC_TIME", "es_ES.UTF-8")
##############################################
### Segunda Gráfica
calendarPlot(tmaxs_df, pollutant="tmaxs_zoo", year=2017, annotate = "value",
             main="Temperatura máxima diaria (2017) \n Estación Hermosillo Observatorio", 
             breaks = c(0, 38, 40, 42, 44, 50), 
             labels = c("<=38", "38", "40", "44", ">44"),
             key.header = "Grados C.")
##############################################
library(ggplot2)
tmaxs1_df <- data.frame(dias1=c(coredata(Gthan38), coredata(Gthan40), coredata(Gthan42), coredata(Gthan44)),
years1=rep(2000:2017,4), grados1 = c(rep(38,18), rep(40,18), rep(42,18), rep(44,18)),
colores= c(rep("#ffdb6c",18), rep("#ff831e",18), rep("#f50000",18), rep("#8b001e",18)))
tmaxs1_df
##############################################
# Funcion tomada de 
### https://stackoverflow.com/users/2372064/mrflick
### https://stackoverflow.com/questions/46327431/ggplot2-and-facet-wrap-add-geom-hline
StatMeanLine <- ggproto("StatMeanLine", Stat,
                        compute_group = function(data, scales) {
                          transform(data, yintercept=mean(y))
                        },
                        required_aes = c("x", "y")
)

stat_mean_line <- function(mapping = NULL, data = NULL, geom = "hline",
                           position = "identity", na.rm = FALSE, show.legend = NA, 
                           inherit.aes = TRUE, ...) {
  layer(
    stat = StatMeanLine, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}
##############################################
#### Tercera gráfica
ggplot(tmaxs1_df, aes(y=dias1, x=years1, fill=factor(grados1))) +
  geom_bar( stat="identity") + 
  facet_grid(grados1 ~ .) +
  scale_fill_manual(values = as.character(unique(tmaxs1_df$colores)),
  name="Umbral °C \n T. Max.") + xlab("Años") +
  ylab("Días con T. Max. >") +
  scale_x_continuous(breaks = round(seq(2000, 2017, by = 2),1))  +
  stat_mean_line(color="black")
#####################################################################