TALLERES CLIMATOLOGÍA

INTRODUCCION

El objetivo de esta guía es tener acceso rápido a los contenidos prácticos que trataremos en el curso de Climatología 1 y 2. Comenzaremos desde la base de que los estudiantes dóminan y aplican los fundamentos de la programación en lenguaje R, sin embargo, existirá un pequeño recordatorio para refrescar la memoria.

En esta guía podrán encontrar un manual paso a paso de las actividades que realizaremos en los talleres, que pueden o no ser calificados.

Los archivos necesarios para el desarrollo de esta guía práctica los podran encontrar en el siguiente GitHub.

Comenzamos recordando que para instalar un paquete de datos usaremos la función install.packages() y dentro del paréntesis escribiremos el nombre del paquete entre comillas:

install.packages("tidyverse")

Solo se deben instalar una vez los paquetes. Luego, si los necesitamos solo debemos “llamarlos” o activarlos. Para ello usaremos la función library() y dentro del paréntesis escribiremos el nombre del paquete, pero esta vez sin comillas:

library(tidyverse)

## Warning: package 'lubridate' was built under R version 4.4.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

REPASO CREACION DE GRAFICOS

GRAFICOS SIMPLES

Para recordar como trabajar con los gráficos vamos a utilizar una tabla que contiene la precipitación registrada por una estación cercana a la ciudad de Frutillar (frutillar_pps.xls). Esta tabla contiene la precipitación en milímetros (mm) expresada en valores mensuales. Contiene 5 columnas: la primera nos indica el año (year), la segunda el mes (month), el día (day), la fecha completa (date) y la cifra de precipitación (valor):

El primer ejercicio será leer una tabla utilizando la funcion read_excel() de la libreria readxl:

frutillar <- readxl::read_excel("TuDirectorioPC/frutillar_pps.xls")

Luego revisaremos el contenido en R usando la función View():

## Warning: package 'DT' was built under R version 4.4.3

a. Gráficos de Barra

Comenzaremos manipulando nuestra tabla reduciendo el número de datos. Para ello vamos a filtrar usando la función filter() del paquete “dplyr”. Puede ser el año que deseen.

a2015 <- filter(frutillar, year == 2015)
# Nota: si la columna year es numerica escribimos el numero,
# si es de clase caracter lo escribimos entre comillas.

Luego, usaremos la función barplot() dónde indicaremos que columna vamos a plotear, recuerden pueden utilizar los paréntesis de corchetes (a2015[,5]) o la opción de subset con el signo peso (a2015$valor):

barplot(a2015$valor)

Este es el plot mas sencillo de todos, dónde viene todo por defecto. Ahora lo que haremos será cambiar el color de las barras a un tono azul:

barplot(a2015$valor, col = "blue")

Siguiente paso será editar los límites del eje y. Actualmente están por defecto ajustándose a los datos, yendo de 0 a 300. Nosotros lo cambiaremos de 0 a 350.

barplot(a2015$valor, col = "blue", ylim = c(0, 350))

En estos momentos nuestro gráfico cuenta con un eje y definido dónde aparecen los valores de precipitación, pero en el eje x no aparece nada entonces, los agregaremos.

barplot(a2015$valor, col = "blue", ylim = c(0, 350), names = a2015$month)

Si somos más exigentes podríamos cambiar el eje x para que en lugar de aparecer el numero del mes nos aparezca al menos la inicial del mes. Para ello crearemos un vector con 12 caracteres con la inicial de cada mes y luego lo indicaremos en el argumento names= de nuestra función barplot():

mes <- c("E", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
barplot(a2015$valor, col = "blue", ylim = c(0, 350), names = mes)

Y para finalizar, vamos a agregarle algo muy importante: las etiquetas de los ejes x e y, que son vitales para poder interpretar un gráfico que no está contextualizado (sin su explicación por escrito o expuesto).

barplot(a2015$valor,
  col = "blue", ylim = c(0, 350), names = mes,
  ylab = "Precipitacion mm", xlab = "Meses"
)

Recuerden que pueden revisar más argumentos consultado desde la consola con ?barplot, o directo desde “Help”:

# ?barplot # barplot es función

b. Gráficos de puntos

Como estamos trabajando con una variable con una connotación temporal es necesario revisar si cuenta con la columna de tiempo (date) en formato de tiempo (fecha/Date). Para hacer esa consulta ocuparemos la función class() que ya hemos visto antes, que nos indicará que tipo o clase es nuestro objeto, o nuestra columna (si se la indicamos), para eso escribimos en nuestra consola:

class(frutillar$date)

## [1] "character"

Para transformar esa columna en formato fecha vamos a usar la función as.Date() como se indica:

frutillar$date <- as.Date(frutillar$date, format = "%Y-%m-%d")
# Nota: el argumento format va a depender de como vengan nuestros datos de fechas en la tabla,
# ej: yyyy/mm/dd sería %Y/%m/%d

class(frutillar$date)

## [1] "Date"

Ahora estamos seguros de que es formato fecha podemos plotear. En este caso debemos indicar explícitamente que corresponde al eje x y que al eje y utilizando la función plot():

plot(x = frutillar$date, y = frutillar$valor)

Lo siguiente que haremos será cambiar el formato de los puntos:

plot(x = frutillar$date, y = frutillar$valor, pch = 20)

Nota 1

Por último, cambiaremos el color y las etiquetas de los ejes x e y, de la misma forma que en el gráfico de puntos:

plot(
  x = frutillar$date, y = frutillar$valor,
  pch = 20, col = "cyan", ylab = "Precipitacion (mm)", xlab = "Años"
)

c. Gráficos de línea

Igual que en el gráfico de puntos, nos debemos fijar si nuestras fechas tienen la clase fecha dentro de nuestra tabla, sino debemos cambiarlo. Posterior a ello, usaremos la misma función plot, pero incorporaremos el argumento type= que nos permite definir si queremos que nuestro gráfico sea de línea.

plot(x = frutillar$date, y = frutillar$valor, type = "l") # en type= “l” (ele) de line

El siguiente paso será cambiarle el tipo de línea, es decir continua, punteada, con guiones, etc.

plot(x = frutillar$date, y = frutillar$valor, type = "l", lty = 2)

Nota 2

Por último, cambiaremos el color y las etiquetas de los ejes x e y, de la misma forma que en el gráfico de lineas:

plot(
  x = frutillar$date, y = frutillar$valor, type = "l",
  lty = 2, col = "blue", ylab = "Precipitacion (mm)", xlab = "Años"
)

d. Gráficos de lineas y puntos

Aquí lo que debemos hacer es primero plotear el gráfico de puntos, y luego para agregar el siguiente plot (de líneas) usaremos la función lines() cuyos argumentos son los mismos, pero debemos considerar no repetir los argumentos que sean globales, como el nombre de los ejes. Además, el argumento type= ya no es necesario porque la función lines() plotea líneas que se agregan a otro gráfico, en este caso se agregan al gráfico de puntos. Debemos correr una línea después de la otra.

plot(
  x = frutillar$date, y = frutillar$valor,
  pch = 20, col = "brown", ylab = "Precipitacion (mm)", xlab = "Años"
)
lines(x = frutillar$date, y = frutillar$valor, col = "blue")

Si queremos cambiar algún argumento sobre la línea, como por ejemplo el tipo de línea, debemos volver a correr la línea del plot de puntos, ya que lines() plotea encima.

Por ejemplo, si después de haber corrido las líneas anteriores corremos nuevamente la línea de lines() con otros argumentos, se pondrá encima, por ejemplo:

lines(x = frutillar$date, y = frutillar$valor, col = "green", lty = 2)

Ejemplo de superposición de plots

e. Histogramas o Gráficos de frecuencia

Al igual que en el gráfico de barras solo debemos indicarle la columna dónde se encuentran los datos, porque automáticamente el gráfico calcula el eje y, que recordemos que es frecuencia.

hist(frutillar$valor)

De manera automática, el eje x determinó un intervalo de menor a mayor de precipitaciones, y dividió las barras cada 50 mm de precipitación, o sea son 10 barras con 9 cortes. Esto quiere decir que la primera barrita representa la frecuencia de los datos cuya precipitación esté entre los 0 y los 50 mm. Si nosotros queremos cambiar el número de barras, podemos hacerlo con el argumento breaks=:

hist(frutillar$valor, breaks = 20)

Ahora tenemos 21 barras con 20 cortes. Este argumento tiene que ver con los datos, por lo que puede que muchas veces indiquemos un numero de cortes, pero no cambie porque los datos no lo permiten; no se le puede exigir más o menos de lo que datos representan.

Ejemplo, si tenemos los valores 5, 20, 35 y 40 cada uno tiene una frecuencia de 4, 14, 9 y 20, tendríamos un gráfico así:

Ejemplo de resultado

No podemos pedirle 20 cortes porque los datos no lo permiten.

Por último, incorporaremos el color, etiquetas del eje x e y, y además le cambiaremos el titulo con el argumento main=.

hist(frutillar$valor,
  breaks = 20, col = "red", xlab = "Precipitaciones (mm)", ylab = "Frecuencia",
  main = "Histograma de Precipitaciones; breaks=20"
)

f. Boxplot o Gráfico de caja con bigotes

Cuando queremos representar solo una caja debemos utilizar la función boxplot() e indicar la columna de los datos, igual que con las barras y el histograma.

boxplot(frutillar$valor)

Pero si queremos representar variabilidad por grupos, debemos indicarles cual columna contiene esos grupos que serán representados en el eje x, y la cual columna tiene los valores que serán representados en el eje y, usando la formula con la expresión: boxplot(y~x)

boxplot(frutillar$valor ~ frutillar$month)

Podemos cambiar el color de todas las cajas indicando un único color, o podemos hacer una concatenación de colores usando la función concatenar c() cuya extensión (cantidad) sea la misma que el número de cajas.

boxplot(frutillar$valor ~ frutillar$month, col = "blue")

Concatenando:

boxplot(frutillar$valor ~ frutillar$month, col = c(
  "blue", "yellow", "green", "lightblue",
  "brown", "pink", "purple", "orange",
  "cyan", "darkgreen", "grey", "gold"
))

También existen funciones que entregan colores en paletas de colores según nosotros le indiquemos. Luego finalizamos:

boxplot(frutillar$valor ~ frutillar$month,
  col = rainbow(12)
) # el valor dentro del parentesís es el número de cajas

boxplot(frutillar$valor ~ frutillar$month,
  col = rainbow(12),
  names = mes, xlab = "Meses", ylab = "Precipitacion (mm)"
)

Nota 3

g. Varios plots en un mismo gráfico

Para este ejemplo filtraremos igual que para el gráfico de barras 4 años diferentes, y los guardaremos en nuevos objetos.

library(dplyr) # activa los paquetes o librerías para acceder a sus funciones
a2011 <- filter(frutillar, year == 2011)
a2012 <- filter(frutillar, year == 2012)
a2013 <- filter(frutillar, year == 2013)
a2014 <- filter(frutillar, year == 2014)

Luego, bajo la misma lógica que cuando ploteamos puntos y líneas, debemos crear un primer plot de línea usando la función plot() y el argumento type=. Luego, usamos la función lines() para los otros 3 años restantes y al final legend() para agregar la leyenda.

plot(
  x = a2011$month, y = a2011$valor, type = "l", col = "red",
  xlab = "Meses", ylab = "Precipitacion (mm)", ylim = c(0, 400)
)
lines(a2012$valor, col = "blue")
lines(a2013$valor, col = "green")
lines(a2014$valor, col = "orange")
legend(
  x = 1, y = 400, legend = c("Año 2011", "Año 2012", "Año 2013", "Año 2014"),
  col = c("red", "blue", "green", "orange"), lty = c(1, 1, 1, 1), cex = 0.8
)

En legend(), x= e y= corresponden a la posición; hay que probar distintas combinaciones hasta que se ajuste dónde queremos, legend= corresponde al texto de nuestra leyenda, concatenamos 4 textos correspondiente a los años seleccionados. En col= también debemos concatenar los colores de nuestras líneas, en el mismo orden en que ploteamos. El argumento lty= nos dice que tipo de línea (punteada, guiones, etc.) escogimos en nuestro plot, nuevamente concatenamos tantos como líneas tengamos, y finalmente cex= es el tamaño del cuadro de la leyenda, aquí puesto para que quede al 80% de su tamaño real.

h. Varios plots en una misma lámina pero por separado

Igual que en el punto anterior, filtraremos igual que para el gráfico de barras 4 años diferentes, y los guardaremos en nuevos objetos.

library(dplyr) # activa los paquetes o librerías para acceder a sus funciones
a2011 <- filter(frutillar, year == 2011)
a2012 <- filter(frutillar, year == 2012)
a2013 <- filter(frutillar, year == 2013)
a2014 <- filter(frutillar, year == 2014)

Ahora, antes de plotear, usaremos la función par() que nos permiten dividir nuestro plot en áreas de ploteo según cuantas filas y columnas indiquemos. Como plotearemos 4 gráficos pondremos que sea 2 filas y 2 columnas de la siguiente forma:

par(mfrow = c(2, 2))

Luego debemos usar la función plot() indicando el tipo con el argumento type= para cada uno con sus otros argumentos que ajustan, debido a que cada plot será un gráfico por serado.

plot(
  x = a2011$month, y = a2011$valor, type = "l", col = "red",
  xlab = "", ylab = "Precipitacion (mm)", ylim = c(0, 400), main = 2011
)
plot(
  x = a2012$month, y = a2012$valor, type = "l", col = "blue",
  xlab = "", ylab = "", ylim = c(0, 400), main = 2012
)
plot(
  x = a2013$month, y = a2013$valor, type = "l", col = "green",
  xlab = "Meses", ylab = "Precipitacion (mm)", ylim = c(0, 400), main = 2013
)
plot(
  x = a2014$month, y = a2014$valor, type = "l", col = "orange",
  xlab = "Meses", ylab = "", ylim = c(0, 400), main = 2014
)

Si nos equivocamos en el número de filas o columnas debemos correr nuevamente par() o si queremos dejar de plotear varios gráficos por lámina debemos escribir en la consola dev.off().

GRAFICOS AVANZADOS CON GGPLOT2

Ya aprendidos como hacer gráficos simples utilizando paquetes nativos de R. Ahora daremos un paso más adelante en dificultad y trabajaremos con un paquete de datos que está especializado en crear gráficos a partir de tablas. Este paquete es uno de los más utilizados por su gran potencial de creación y edición de gráficos, su nombre es GGPLOT2.

Vamos a trabajar los mismos gráficos y con la misma base de datos para que se puedan hacer comparaciones de como se crean los mismos gráficos con códigos diferentes. Al ir avanzando en cada gráfico, se hará énfasis en líneas de códigos nuevas que vayan apareciendo.

Para crear gráficos avanzados con esta base de datos instalaremos el paquete GGPLOT2:

install.packages("ggplot2")

library(ggplot2)

a. Gráficos de barra

Para comenzar, filtraremos igual que en el práctico anterior usando la función filter() del paquete “dplyr”. Puede ser el año que deseen.

library(dplyr)
a2017 <- filter(frutillar, year == 2017)

Luego escribiremos las primeras líneas que nos permiten crear nuestro gráfico usando la librería ggplot2. Todos los gráficos de GGPLOT2 comenzarán utilizando la función ggplot, en la que designaremos cual es la base de datos que vamos a utilizar:

ggplot(a2017) +

Algo muy importante, es que ahora cada línea de código va a finalizar con un signo más (+). Esto quiere decir que la siguiente línea va a corresponder a un elemento de diseño y se le va a agregar (por eso el +) a la línea anterior.

El siguiente elemento va a ser utilizar la función que designa que es un gráfico de barras, esta es la función geom_bar(). Dentro de ella, definiremos algunos parámetros del gráfico, como el eje x e y, para ello usaremos la función aes():

ggplot(a2017) +
  geom_col(aes(x = month, y = valor))

Lo primero que podemos advertir de este código es que la segunda línea posee sangría, eso también es un indicador que esa línea es la continuación de la primera.

ggplot(a2017) +
  geom_col(aes(x = month, y = valor))

Ahora, agregaremos algunos argumentos que nos permiten cambiar el color del borde de las barras (color=) y el color de relleno (fill=):

ggplot(a2017) +
  geom_col(aes(x = month, y = valor), color = "blue", fill = "lightblue")

El siguiente paso será cambiarle los límites a nuestro eje x, utilizando la función ylim().

ggplot(a2017) +
  geom_col(aes(x = month, y = valor), color = "blue", fill = "lightblue") +
  ylim(0, 600)

Ahora cambiaremos las etiquetas de los ejes x e y usando las funciones xlab() e ylab() respectivamente:

ggplot(a2017) +
  geom_col(aes(x = month, y = valor), color = "blue", fill = "lightblue") +
  ylim(0, 600) +
  xlab("Meses") +
  ylab("Precipitación (mm)")

Si se fijan, en la tercera línea podemos seguir incorporando elementos con el +. También lo podemos hacer en la siguiente línea, pero debemos tener presente que debe estar con sangría. De ambas maneras es válido y funciona.

Lo último será un poco más avanzado, y corresponde a cambiar el texto de cada columna (eje x) y su orientación a 90° para que no se superpongan una sobre las otras.

ggplot(a2017) +
  geom_col(aes(x = month, y = valor), color = "blue", fill = "lightblue") +
  ylim(0, 600) +
  xlab("Meses") +
  ylab("Precipitación (mm)") +
  scale_x_discrete(labels = c("Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic")) +
  theme(axis.text.x = element_text(angle = 90))

b. Gráficos de Lineas y puntos

Recordemos que cuando estamos trabajando con una variable con una connotación temporal es necesario revisar si cuenta con la columna de tiempo (date) en formato de tiempo (fecha/Date). Para hacer esa consulta ocuparemos la función class() que ya hemos visto antes, que nos indicará que tipo o clase es nuestro objeto, o nuestra columna (si se la indicamos), para eso escribimos en nuestra consola:

class(frutillar$date)

## [1] "Date"

Para transformar esa columna en formato fecha vamos a usar la función as.Date() como se indica:

frutillar$date <- as.Date(frutillar$date, format = "%Y-%m-%d")

class(frutillar$date)

## [1] "Date"

Ahora podemos comenzar con el gráfico y lo haremos igual que la vez pasada. Utilizando como primera función ggplot() y para crear las líneas usaremos la función geom_line() y para los puntos geom_point().

ggplot(frutillar, aes(x = date, y = valor)) +
  geom_line() +
  geom_point()

Ahora cambiaremos el tipo de línea y el tipo de punto con algunos parámetros dentro de cada función respectiva:

ggplot(frutillar, aes(x = date, y = valor)) +
  geom_line(linetype = "dashed") +
  geom_point(shape = 18)

Finalmente cambiaremos el color de cada elemento y los textos en los ejes.

ggplot(frutillar, aes(x = date, y = valor)) +
  geom_line(linetype = "dashed", colour = "darkblue") +
  geom_point(shape = 18, colour = "cyan") +
  xlab("Años") +
  ylab("Precipitación (mm)")

c. Histogramas o Gráficos de frecuencia

Para el histograma usaremos la función geom_histogram() y solo debemos indicarle el eje x, ya que el eje y siempre será frecuencia:

ggplot(frutillar) +
  geom_histogram(aes(valor))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Luego, usaremos el argumento binwidth= para definir el tamaño de las barras (el ancho del intervalo de cada barra), fill= para cambiar el color, y además agregaremos las etiquetas de los ejes.

ggplot(frutillar) +
  geom_histogram(aes(valor), binwidth = 10, fill = "red") +
  xlab("Precipitación (mm)") +
  ylab("Frecuencia")

d. Boxplot o Gráfico de caja con bigotes

Comenzamos igual que los otros gráficos, y esta vez usaremos la función geom_boxplot(). Adicionalmente usaremos el argumento Group= en nuestra función aes() cuando queramos incorporar varios boxplot en un mismo gráfico.

ggplot(frutillar) +
  geom_boxplot(aes(y = valor))

ggplot(frutillar) +
  geom_boxplot(aes(x = month, y = valor, group = month))

Ahora agregaremos color y cambiaremos las etiquetas:

ggplot(frutillar) +
  geom_boxplot(aes(x = month, y = valor, group = month), outlier.color = "red", fill = topo.colors(12)) +
  xlab("Meses del año") +
  ylab("Precipitación (mm)")

e. Guía rápida de colores

La siguiente imagen nos presenta una serie de gamas y tonalidades de colores con sus respectivos nombres para poder usarlos a la hora de definir un color.

Recurso: http://sape.inf.usi.ch/sites/default/files/ggplot2-colour-names.png

En lugar de escribir el nombre de un color, tambien lo podemos usar utilizando su código html, el cual podemos consultarlo en la siguiente web: https://htmlcolorcodes.com/es/

Ejemplo de código html de un color

Tambien existen paletas de colores ya definidas como “topo.cor” pero en otro paquete. Este se llama “R Color Brewer”

Recurso: https://www.r-graph-gallery.com/38-rcolorbrewers-palettes_files/figure-html/thecode-1.png

Para utilizarla, debemos instalar el paquete:

install.packages("RColorBrewer")

ANALISIS EXPLORATORIO DE DATOS

Para comenzar vamos a activar un paquete que agrupa una serie de librerias que permiten realizar analisis exploratorio de datos:

library(tidyverse)

El primer ejercicio será leer una tabla utilizando la funcion read_csv()

tabla <- read_csv("TuDirectorioPC/MET_VALPO.csv")

## Rows: 744 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (6): Year, Month, PP, Tmin, Tmed, Tmax
## date (1): Date
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Para conocer nuestros datos, aplicaremos funciones que nos permitan describirlos a través de la estadística:

Estadística descriptiva

mean(tabla$PP) # promedio

## [1] 34.80914

median(tabla$PP) # mediana

## [1] 10

sd(tabla$PP) # desviacion estandar

## [1] 57.71894

summary(tabla$PP) # resumen

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    3.00   10.00   34.81   43.00  458.00

Distribución de frecuencia

ggplot(tabla, aes(x = PP)) + # definimos variable
  geom_histogram() + # definimos el tipo de gráfico
  ylab("Frecuencia") + # etiqueta eje y
  xlab("Precipitación (mm)") + # etiqueta eje x
  theme_bw() + # cambio de tema
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Variabilidad datos (boxplot)

ggplot(tabla, aes(y = PP)) + # definimos variable
  geom_boxplot() + # definimos el tipo de gráfico
  xlab("Precipitación (mm)") + # etiqueta eje x
  theme_bw() + # cambio de tema
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

Variabilidad datos por año (boxplot)

ggplot(tabla, aes(y = PP, x = Year, group = Year)) + # definimos variable
  geom_boxplot() + # definimos el tipo de gráfico
  ylab("Precipitación (mm)") + # etiqueta eje x
  xlab("") +
  theme_bw() + # cambio de tema
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

Variabilidad datos por mes (boxplot)

tabla$Month <- as.factor(tabla$Month) # transformamos variable numerica a categorica
tabla$month_label <- month(tabla$Date, label = T)
ggplot(tabla, aes(y = PP, x = Month, group = Month)) + # definimos variable
  geom_boxplot() + # definimos el tipo de gráfico
  ylab("Precipitación (mm)") + # etiqueta eje x
  xlab("") +
  theme_bw() + # cambio de tema
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

Variabilidad datos por mes en el tiempo

ggplot(tabla, aes(y = PP, x = Year)) + # definimos variable
  geom_line() + # definimos el tipo de gráfico
  geom_smooth(color = "red") +
  # stat_spline(color = "red") +
  ylab("Precipitación (mm)") + # etiqueta eje x
  xlab("") +
  facet_wrap(~month_label, nrow = 1, ncol = 12) +
  theme_bw() + # cambio de tema
  # invertimos la posicion de las etiquetas del eje
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(tabla, aes(y = Tmed, x = Year)) + # definimos variable
  geom_line() + # definimos el tipo de gráfico
  geom_smooth(color = "red") +
  # stat_spline(color = "red") +
  ylab("Temperatura (°C)") + # etiqueta eje x
  xlab("") +
  facet_wrap(~month_label, nrow = 1, ncol = 12) +
  theme_bw() + # cambio de tema
  # invertimos la posicion de las etiquetas del eje
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  # editamos el tema
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

REMUESTREO TEMPORAL

El remuestreo temporal nos permite realizar agregaciones temporales de nuestros datos, es decir, permite agrupar nuestros valores de acuerdo a una condición temporal. Por ejemplo, si tenemos datos diarios (365 al año), podemos agregar nuestros datos en mensuales (12 al año), semestrales (2 al año), anuales (1 al año), etc. Esto nos permite poder observar nuestros datos desde otra perspectiva o tambien permite resumir el volumen de nuestros datos. Tambien es util cuando la naturaleza de nuestros datos tienen detalles temporales diferentes.

Datos Anuales: Precipitación

year_pp <- tabla %>%
  group_by(Year) %>% # agrupamos segun la categoria de "Año"
  summarise(PP_yearly = sum(PP)) # definimos que estadistico usaremos para agrupar nuestros datos

ggplot(year_pp, aes(x = Year, y = PP_yearly)) +
  geom_col(fill = "blue") +
  ylab("Precipitación (mm)") +
  xlab("") +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

Datos Anuales: Temperatura

year_tmed <- tabla %>%
  group_by(Year) %>% # agrupamos segun la categoria de "Año"
  summarise(Tmed_yearly = mean(Tmed)) # definimos que estadistico usaremos para agrupar nuestros datos

ggplot(year_tmed, aes(x = Year, y = Tmed_yearly)) +
  geom_line(color = "red") +
  geom_point(size = 1, color = "red") +
  ylab("Temperatura media (°C)") +
  xlab("") +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  )

Datos Mensuales: Creación de Climograma

# a. Tabla mensual de precipitacion
climo_pp <- tabla %>%
  group_by(Month) %>% # agrupamos segun la categoria de "Mes"
  summarise(PP_clim = mean(PP)) # definimos que estadistico usaremos para agrupar nuestros datos

# b. Tabla mensual de temperatura
climo_tmed <- tabla %>%
  group_by(Month) %>%
  summarise(Tmed_clim = mean(Tmed))

Ahora juntamos ambos objetos de tablas mensuales en una misma tabla:

tabla_clim <- climo_pp
tabla_clim$Tmed_clim <- climo_tmed$Tmed_clim

Con esta información, ya podemos generar nuestra climograma con GGPLOT2

ggplot(tabla_clim, mapping = aes(group = 1)) +
  # Definimos el eje X e Y de la precipitación (eje Y izquierdo)
  geom_col(aes(x = Month, y = PP_clim), fill = "blue") +
  # Definimos el eje X e Y de la precipitación (eje Y derecho)
  # Debemos reescalar nuestros datos de temperatura para que sea equivalente con
  # la escala de la PP
  geom_line(mapping = aes(x = Month, y = Tmed_clim * 10), color = "red") +
  geom_point(mapping = aes(x = Month, y = Tmed_clim * 10), color = "red", size = 1.5) +
  # En esta sección volvemos a reescalar los datos de Temperatura a tu valor original
  # por eso realizamos operación inversa
  scale_y_continuous("Precipitation [mm]",
    sec.axis = sec_axis(~ . / 10, name = "Temperature [°C]")
  ) +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  ) +
  xlab("")

Climate Stripe

Para crear el climate stripe, utilizaremos la tabla de temperatura anual para poder visualizar de manera intuitiva, si la temperatura ha ido en aumento o ha disminuido. Para este ejercicio utilizaremos como referencia el periodo completo de nuestros datos, pero lo ideal es escoger un periodo acotado de referencia para determinar los cambios.

Primero calcularemos las fluctuaciones respecto a la media. Para ello utilizaremos la función scales.

year_tmed$scale <- scale(year_tmed$Tmed_yearly, scale = T)

Luego definiremos una paleta de colores que nos permitirá visualizar los años más fríos respecto a la media de tonos azules, y los años más calidos respecto a la media, de tonos rojizos. Para ello nos apoyaremos de la librería “RColorBrewer”.

library(RColorBrewer)
col_strip <- brewer.pal(11, "RdBu")

Ahora crearemos un “tema” para nuestro plot.

theme_strip <- theme_minimal() +
  theme(
    axis.text.y = element_blank(),
    axis.line.y = element_blank(),
    axis.title = element_blank(),
    panel.grid.major = element_blank(),
    legend.title = element_blank(),
    axis.text.x = element_text(vjust = 3),
    panel.grid.minor = element_blank(),
    plot.title = element_text(size = 14, face = "bold")
  )

Y ahora podemos plotear con GGPLOT2:

ggplot(year_tmed, aes(x = Year, y = 1, fill = scale)) +
  geom_tile() +
  scale_y_continuous(expand = c(0, 0)) +
  # revertimos la paleta de colores para darle sentido usando rev()
  scale_fill_gradientn(colors = rev(col_strip)) +
  # definimos años que aparecerán en el plot
  scale_x_continuous(breaks = seq(1960, 2020, 5)) +
  # agregamos la fuente de nuestros datos
  labs(
    title = "TEMPERATURA REGION DE VALPARAISO 1960-2021",
    caption = "Datos: CR2MET"
  ) +
  theme_strip # indicamos el tema que creamos anteriormente

OPERACIONES SOBRE UN RASTER

Operaciones de consulta

library(terra)
library(tidyverse)
cr2met_pp <- rast("C:/TuDirectorio/CR2MET_PP_Mensual_1960-2021.tif")
plot(cr2met_pp[[1]]) # con doble parantesis indicamos el n° de la banda que queremos plotear

## Warning: package 'terra' was built under R version 4.4.2

## terra 1.8.5

## 
## Adjuntando el paquete: 'terra'

## The following object is masked from 'package:tidyr':
## 
##     extract

Clase de objeto

class(cr2met_pp)

## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"

Dimensiones

c(nrow(cr2met_pp), ncol(cr2met_pp), nlyr(cr2met_pp)) # filas, columnas, bandas

## [1]  70 101 744

Resolución

res(cr2met_pp)

## [1] 0.05 0.05

Extensión

ext(cr2met_pp)

##   xmin   xmax   ymin   ymax 
## -73.40 -68.35 -34.80 -31.30

Sistema de referencia y Coordenadas

crs(cr2met_pp)

## [[1]]
##  [1] "GEOGCRS[\"WGS 84\","                                      
##  [2] "    ENSEMBLE[\"World Geodetic System 1984 ensemble\","    
##  [3] "        MEMBER[\"World Geodetic System 1984 (Transit)\"],"
##  [4] "        MEMBER[\"World Geodetic System 1984 (G730)\"],"   
##  [5] "        MEMBER[\"World Geodetic System 1984 (G873)\"],"   
##  [6] "        MEMBER[\"World Geodetic System 1984 (G1150)\"],"  
##  [7] "        MEMBER[\"World Geodetic System 1984 (G1674)\"],"  
##  [8] "        MEMBER[\"World Geodetic System 1984 (G1762)\"],"  
##  [9] "        MEMBER[\"World Geodetic System 1984 (G2139)\"],"  
## [10] "        MEMBER[\"World Geodetic System 1984 (G2296)\"],"  
## [11] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"      
## [12] "            LENGTHUNIT[\"metre\",1]],"                    
## [13] "        ENSEMBLEACCURACY[2.0]],"                          
## [14] "    PRIMEM[\"Greenwich\",0,"                              
## [15] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"       
## [16] "    CS[ellipsoidal,2],"                                   
## [17] "        AXIS[\"geodetic latitude (Lat)\",north,"          
## [18] "            ORDER[1],"                                    
## [19] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"   
## [20] "        AXIS[\"geodetic longitude (Lon)\",east,"          
## [21] "            ORDER[2],"                                    
## [22] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"   
## [23] "    USAGE["                                               
## [24] "        SCOPE[\"Horizontal component of 3D system.\"],"   
## [25] "        AREA[\"World.\"],"                                
## [26] "        BBOX[-90,-180,90,180]],"                          
## [27] "    ID[\"EPSG\",4326]]"

Tipo de dato

datatype(cr2met_pp[[1]])

## [1] "INT2S"

Codificación

cr2met_pp

## class       : SpatRaster 
## dimensions  : 70, 101, 744  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -73.4, -68.35, -34.8, -31.3  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : CR2MET_PP_Mensual_1960-2021.tif 
## names       : 1960-01-15, 1960-02-15, 1960-03-15, 1960-04-15, 1960-05-15, 1960-06-15, ... 
## min values  :          0,          0,          0,          0,          0,          6, ... 
## max values  :        298,        168,        195,        105,        102,        631, ...

Operaciones de alteración

Recorte, Clip o Crop

El recorte signica sustraer parte de la matriz ráster según una extensión definida. Podemos crear una extensión ingresando las coordenas límite (xmin, xmax, ymin, ymax), usando la extensión de un geodato, o dibujando una extensión sobre un plot usando la función draw().

vector_valpo <- vect("C:/TuDirectorio/valpo.gpkg")
recorte_PP <- cr2met_pp %>% crop(vector_valpo)
plot(recorte_PP[[1]])
plot(vector_valpo, add = T)

Máscara o Mask

La aplicación de una máscara o enmascaramiento implica el sustraer los valores de un ráster dentro de los límites de un geodato, reclasificanto los pixeles fuera de estos límites como “No Dato”, “not available” o “NA”. Se aconseja recortar y enmascarar al área de interes.

Un elemento importante a considerar es hacerse la pregunta, ¿Qué ocurre cuando el límite de mi área de interés toca parte del pixel pero no lo cubre por completo? ¿lo incluye o no?

Existe un parámetro de tipo lógico dentro la función mask() llamado “touches”. Cuando touches es T si toca el pixel, queda incluido en la máscara, si touches es F el pixel debe estar contenido en un 100% dentro del área de interes para ser incluido.

mascara_completa <- cr2met_pp %>% mask(vector_valpo) # mascara sin clip
mascara_crop <- recorte_PP %>% mask(vector_valpo) # mascara con clip

Operaciones de transformación

Vectorización: contornos o isolineas

nivel <- c(5, 50, 150, 300, 450, 600) # serie de intervalos de isolineas
pp_contorno <- as.contour(cr2met_pp[[6]], levels = nivel) # creamos vector de lineas
plot(cr2met_pp[[6]])
terra::contour(cr2met_pp[[6]], levels = nivel, add = T) # nos permite plotear encima con etiqueta

Extracción de valores de los pixeles

titulo <- paste("Histograma Precipitacion", names(mascara_crop[[6]])) # Titulo plot
hist(mascara_crop[[6]], main = titulo)

valores_elev <- values(mascara_crop) %>% na.omit() # omitimos las celdas NA
mean(valores_elev)

## [1] 34.7065

median(valores_elev)

## [1] 10

Estadística zonal: Precipitacion promedio comunas de la provincia de Valparaíso

comunas_valpo <- vect("TuDirectorio/comunas_valpo.gpkg")
provincia_valpo <- subset(comunas_valpo, comunas_valpo$NOM_PROV == "Valparaíso")
media_PP_comunas <- zonal(mascara_crop, provincia_valpo,
  fun = mean, weights = F, na.rm = T, touches = T, as.polygons = T
)
plot(media_PP_comunas, "1960-06-15", col = heat.colors(10))

tabla_resultado <- as.data.frame(media_PP_comunas)

Estadística zonal: Precipitacion promedio región de Valparaíso

media_PP_valpo <- zonal(mascara_crop, vector_valpo,
  fun = mean, weights = F, na.rm = T, touches = T, as.polygons = T
)
plot(media_PP_valpo, "1960-06-15", col = heat.colors(10))

tabla_valpo <- as.data.frame(media_PP_valpo)

df <- data.frame(PP = as.numeric(tabla_valpo[, 2:745]), fechas = names(mascara_crop) %>% as.Date())

ggplot(df, aes(x = fechas, y = PP)) +
  geom_col(col = "blue")

SERIES DE TIEMPO

setwd("C:/TuDirectorio/") # Tu directorio de trabajo (donde tienes tus archivos)
pp_valpo <- read.csv("PP_Valparaiso.csv")
# transformamos la columna llamada "date" de caracter (chr) a fecha (Date)
pp_valpo$date <- as.Date(pp_valpo$date)

Remuestreo temporal

Opción 1: usando tidyverse

pp_valpo$month <- month(pp_valpo$date) # creamos una nueva columna con los meses
pp_valpo$year <- year(pp_valpo$date) # creamos una nueva columna con los años

tabla_anual <- pp_valpo %>%
  group_by(year) %>% # agrupamos por año y mes
  summarise(PP_anual = sum(PP, na.rm = T)) # agrupamos usando la suma

## 
## Adjuntando el paquete: 'ggpubr'

## The following object is masked from 'package:terra':
## 
##     rotate

Opción 2: usando xts

library(xts) # la versión para raster es "rts"

## Cargando paquete requerido: zoo

## 
## Adjuntando el paquete: 'zoo'

## The following object is masked from 'package:terra':
## 
##     time<-

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################

## 
## Adjuntando el paquete: 'xts'

## The following objects are masked from 'package:dplyr':
## 
##     first, last

xts_PP <- xts(
  pp_valpo$PP, # vector numerico con la variable
  pp_valpo$date
) # vector de fechas
PP_Anual_xts <- apply.yearly(xts_PP, FUN = sum) # agregación mensual

Si deseamos extraer los valores del objeto xts calculado, debemos usar la función as.numeric() para transformarlo en un vector numérico. Las fechas las podemos crear con una secuencia:

fecha <- seq(as.Date("1960-01-15"), as.Date("2021-12-15"), "month")

Descomposición de una serie

ts_mensual <- ts(pp_valpo$PP, # vector con la variable
  start = c(1960, 1), # año y mes de inicio de la serie
  end = c(2021, 12), # año y mes de termino de la serie
  frequency = 12
) # frecuencia de datos al año

dec_PP <- decompose(ts_mensual) # descomposición

Tendencia (T - trend): corresponde al comportamiento de la media y la varianza en el tiempo.

Ciclos (C - cycle): oscilaciones respecto a la tendencia en largas o medias escalas temporales. Como se desprende de la tendencia, algunos autores no la consideran como un componente propiamente tal.

Estacionalidad (S - seasonality): comportamientos determinados en periodos fijos. Corresponden a cambios dentro de unidades de tiempo de manera regular (como cambios intra-anuales).

Aleatorio (R - random): comportamiento sin un patrón determinado. Es un comportamiento practicamente impredecible. Lo podemos encontrar en la literatura también como “erratico”.

plot(dec_PP)

Anomalías respecto a desviaciones de la media (Z-score) usando periodo completo

tabla_anual$Zscore <- scale(tabla_anual$PP_anual, scale = T)

ggplot(
  tabla_anual,
  aes(
    x = year, y = Zscore,
    fill = I(ifelse(tabla_anual$Zscore > 0, "blue", "brown")),
    color = I(ifelse(tabla_anual$Zscore > 0, "blue", "brown"))
  )
) +
  geom_col() +
  ylab("Anomalia PP") +
  xlab("") +
  theme_bw()

## Warning: Use of `tabla_anual$Zscore` is discouraged.
## ℹ Use `Zscore` instead.
## Use of `tabla_anual$Zscore` is discouraged.
## ℹ Use `Zscore` instead.

Taller 1

Descargue la información de temperatura media de la región de Valparaíso (disponible en el Aula Virtual) y desarrolle:

Cree un mapa conceptual, en el programa que desee (ej power point), que explique el proceso que usted debe realizar a esta tabla para poder descomponer la serie de tiempo.
Descomponga la serie de tiempo y describa que ocurre con la tendencia. En terminos genericos ¿aumenta? ¿disminuye? Recuerde construir una respuesta bien desarrollada.
Sin realizar el calculo de test responda: ¿La serie es o no estacionaria? Justifique su respuesta.
Realice un remuestreo temporal por año, y calcule las desviaciones respecto a la media. Plotee su resultado y realice una descripción de lo que observa.

Consideraciones

Puede realizar el taller en pareja o de manera individual. Pero esta decisión se debe mantener en los 3 talleres del curso.
Sea ordenada/o en la construcción de su script, ya que este será solicitado cuando deba entregar su evaluación.
Escriba sus respuestas en un archivo Word de manera ordenada. Incorpore gráficos cuando se le solicita o si su respuesta depende del análisis de uno.

CORRELACIONES Y REGRESION LINEAL

Correlaciones

library(tidyverse)
library(readxl)

emisiones <- read_excel("TuDirectorio/Co2_temp.xlsx")

cor(x = emisiones$co2, y = emisiones$temp, method = "spearman")

## [1] 0.9498891

cor(x = emisiones$co2, y = emisiones$temp, method = "kendall")

## [1] 0.8080259

cor(x = emisiones$co2, y = emisiones$temp, method = "pearson")

## [1] 0.9625997

ggplot(data = emisiones, aes(x = co2, y = temp)) +
  geom_point(color = "blue", size = 2) +
  geom_smooth(method = "lm", formula = y ~ x, se = F, color = "red") +
  theme(text = element_text(family = "Arial", size = 18)) +
  scale_x_continuous(expand = c(0.005, 0.005)) +
  theme(
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black")
  ) +
  ylab("Anomalias de temperatura global") +
  xlab("Emisiones globales de CO2")

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

cor.test(
  x = emisiones$co2, y = emisiones$temp, alternative = "two.side",
  conf.level = 0.95, method = "pearson"
)

## 
##  Pearson's product-moment correlation
## 
## data:  emisiones$co2 and emisiones$temp
## t = 28.201, df = 63, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9392044 0.9770984
## sample estimates:
##       cor 
## 0.9625997

arboles <- trees

Correlogramas

library(corrplot)

## corrplot 0.95 loaded

c_arboles <- cor(arboles)
corrplot(c_arboles)

library(PerformanceAnalytics)

## Warning: package 'PerformanceAnalytics' was built under R version 4.4.3

## 
## Adjuntando el paquete: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

chart.Correlation(arboles)

Regresión Lineal

# formula = Y ~ X
lm(formula = temp ~ co2, data = emisiones)

## 
## Call:
## lm(formula = temp ~ co2, data = emisiones)
## 
## Coefficients:
## (Intercept)          co2  
##    -3.30160      0.01029

El primer valor entregado (-3.3016) corresponde al intercepto, mientras que el segundo (0.01029) la pendiente.

Y = m * X + b

Y = 0.01029 * X + (-3.3016)

¿Qué ocurriria si las concentraciones de CO2 aumentaran a 500 y 600 kton?

0.01029 * 500 + (-3.3016)

## [1] 1.8434

0.01029 * 600 + (-3.3016)

## [1] 2.8724

… y en R?

modelo.lin <- lm(formula = temp ~ co2, data = emisiones)
predict(modelo.lin, newdata = data.frame(co2 = c(500, 600)))

##        1        2 
## 1.845397 2.874796

Taller 2

Descargue la información de Precipitación de Balmaceda y Coyhaique (disponible en el Aula Virtual) y desarrolle:

Cree un mapa conceptual, en el programa que desee (ej power point), que explique el proceso que usted debe realizar para gráficar la relación entre ambas precipitaciones y el calculo de su correlación y regresión lineal para predecir los datos ausentes.
Grafique ambas series de tiempo.
Gráfique la relación entre ambas precipitaciones, colocando en el eje correcto la estación patrón y la estación con valores ausentes.
Calcule la correlación entre ambas estaciones utilizando niveles de confianza de 0.90, 0.95 y 0.99 utilizando pearson.

¿De cuanto es la correlación? ¿Es suficiente?
¿Cómo varían los coeficientes de correlación?
¿Cómo varian los intervalos de confianza?

Utilizando ambas series de tiempo, realice una regresión lineal y complete todos los valores ausentes presentes en la serie correspondiente.