############################################################
#R como un SIG: extracción de datos climáticos de WorldClim#    
############################################################

#Paula N Fergnani 
#Ecología Austral

#Introducción a R como un SIG y 
#Extracción de datos climáticos de WorldClim a partir de un archivo de coordenadas de sitios de muestreo.
#Esta rutina se ha realizado usando la versión de R 4.0.3. 

#Parte 1: descargar los paquetes necesarios, instalarlos y cargarlos a R. Establecer el directorio de trabajo.

#Descargar los paquetes necesarios e instalarlos. Esta acción se realiza por única vez.
install.packages("sp")
install.packages("rgdal") 
install.packages("maptools")
install.packages("raster")
install.packages("rgeos")
install.packages("viridis")

#Cargar a R los paquetes. 
library(sp)
library(rgdal) 
library(rgeos)
library(maptools)
library(raster)
library(viridis)

#Establecer el directorio de trabajo en la carpeta que contiene el archivo con la lista de coordenadas
#Respetar el uso de la doble barra \\ en la dirección
#Este directorio es un ejemplo, modificarlo  manualmente poniendo la ruta correspondiente
setwd("H:\\R como un SIG")


#Parte 2: Cargar las coordenadas y crear un objeto espacial de puntos 

#Cargar y guardar las coordenadas en un marco de datos llamado "coords"
coords=read.table("Localidades.Abundancia.txt",header=T)

#Ver el marco de datos
coords #en la consola
View(coords) #en un apartado

#Convertir las coordenadas a un objeto espacial de puntos con la función "coordinates"
coordinates(coords)=~x+y

#Ver un sumario del objeto espacial de coordenadas "coords"
summary(coords)

#Ver el marco de datos asociado al objeto espacial
#Es necesario usar @data para acceder 
#Observar que el marco de datos solo contiene los datos de Abundancia
coords@data # Muestra el marco de datos asociado en la consola
View(coords@data) #Muestra el marco de datos asociado en un apartado

#Ver las coordenadas del objeto espacial
#Se utiliza @coords para ver esta información
coords@coords

#Especificar el sistema de coordenadas del objeto espacial
#Estas coordenadas fueron tomadas con el sistema WGS84 "+proj=longlat +datum=WGS84"
proj4string(coords)=CRS("+proj=longlat +datum=WGS84")

#Ver las propiedades del objeto espacial
summary(coords)

#Si se quiere conocer la estructura interna del objeto espacial de coordenadas
str(coords) 

#Si se necesita  agregar datos "manualmente" a la tabla de atributos
#Se crea el vector "Ambiente" dentro del marco de datos accediendo con @data, y con la función "c", se especifica el tipo de ambiente en un vector
coords@data$Ambiente=c("bosque","bosque","bosque","bosque","bosque","bosque","bosque","bosque","bosque","bosque","bosque", "bosque","bosque","bosque","bosque","bosque","bosque","matorral","matorral","matorral","matorral","matorral","matorral","matorral",
                       "matorral","matorral", "estepa", "estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa","estepa")

#ver nuevamente el marco de datos y observar que se creó el vector "Ambiente"
View(coords@data)


#Parte 3: Graficar un mapa para verificar que las coordenadas están correctamente situadas

#Graficar el objeto espacial de coordenadas, el argumento axes=TRUE muestra los valores de los ejes X e Y
plot(coords, axes=TRUE, main="Mis sitios")

#Graficar nuevamente cambiando argumentos del gráfico (pch=tipo de punto, col= color, cex=tamaño del punto)
plot(coords, pch=20, col="red", cex=4, main="Mis sitios", axes=TRUE)

#Graficar solo los sitios de bosque coords@data$Ambiente=="bosque"
#con el argumento add = TRUE, el mapa se grafica sobre el anterior
plot(coords[ coords@data$Ambiente=="bosque", ], col = "turquoise", add = TRUE, axes=TRUE) 

#Ver la ubicación de los puntos respecto del mapa del mundo
#cargar el mapa del mundo que está en el paquete maptools
data(wrld_simpl) 
#Graficarlo
plot(wrld_simpl, axes=TRUE)
#Agregar las coordenadas
plot(coords, pch=20, col="red",  main="Mis sitios", add=TRUE)

#Graficar solo Argentina

#Identificar en la tabla de atributos de wrld_simpl de que forma se consignan los países
#Para ello utilizar la función head() que muestra las primeras líneas de la tabla de atributos
head(wrld_simpl@data) #Ver que en la variable NAME están los paises

#Graficar Argentina 
plot(wrld_simpl[ wrld_simpl@data$NAME=="Argentina", ])

#Agregar las coordenadas de muestreo
plot(coords, pch=20, col="red",  main="Mis sitios", add=TRUE)


#Parte 4: explorar las variables bio de la base de datos climática. Obtener los datos de WorldClim mundiales y visualizarlos usando una escala de celda grande (10 minutos, aproximadamente 340 Km2)

#Establecer el url  de descarga en el link que provee la base de datos de Worldclim en su página web 
#Obvervar en el url que la versión a descargar es la "v2.1" y la resolución 10m "wc2.1_10m_bio". 
#Es importante notar que cuando cambie la versión de la base (cuando se actualice), cambiará también el nombre y el link, por lo que será necesario corregir ambas cosas en todos los lugares en que se mencione el nombre en el script
W_url_10m <- "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_bio.zip"

#Descargar en el directorio de trabajo las variables bio a 10 m de resolución, esto puede demorar unos minutos
download.file(W_url_10m, destfile="wc2.1_10m_bio.zip") 

#Descomprimir el archivo zip
unzip("wc2.1_10m_bio.zip")

#Aclaración es posible descargar las mismas variables en otras 3 resoluciones: 5m, 2.5m, y 30 segundos. 
#Para trabajar con otras resoluciones, solo es necesario cambiar el término "10m" por el que corresponda
#Se deja un ejemplo para la descarga a una resolución de 2.5m
#Si se desea descargar, es necesaio borrar el "#" en cada línea
#Más adelante en el script se descargará la base en la resolución de 30 segundos

#W_url_2.5m <- "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_2.5m_bio.zip"
#download.file(W_url_2.5m, destfile="wc2.1_2.5m_bio.zip")
#unzip("wc2.1_2.5m_bio.zip") #Pesa aproximadamente 643MB y la descarga puede demorar

#Cargar las variables descargadas a R
bio1_10 = raster("wc2.1_10m_bio_1.tif")
bio2_10 = raster("wc2.1_10m_bio_2.tif")
bio3_10 = raster("wc2.1_10m_bio_3.tif")
bio4_10 = raster("wc2.1_10m_bio_4.tif")
bio5_10 = raster("wc2.1_10m_bio_5.tif")
bio6_10 = raster("wc2.1_10m_bio_6.tif")
bio7_10 = raster("wc2.1_10m_bio_7.tif")
bio8_10 = raster("wc2.1_10m_bio_8.tif")
bio9_10 = raster("wc2.1_10m_bio_9.tif")
bio10_10 = raster("wc2.1_10m_bio_10.tif")
bio11_10 = raster("wc2.1_10m_bio_11.tif")
bio12_10 = raster("wc2.1_10m_bio_12.tif")
bio13_10 = raster("wc2.1_10m_bio_13.tif")
bio14_10 = raster("wc2.1_10m_bio_14.tif")
bio15_10 = raster("wc2.1_10m_bio_15.tif")
bio16_10 = raster("wc2.1_10m_bio_16.tif")
bio17_10 = raster("wc2.1_10m_bio_17.tif")
bio18_10 = raster("wc2.1_10m_bio_18.tif")
bio19_10 = raster("wc2.1_10m_bio_19.tif")

#Graficar la variable bio1 (temperatura media anual)
plot(bio1_10, main="Temperatura media anual")

#Realizar un histograma en base a los valores del raster
#Utiliza un porcentaje de las celdas, para realizar el gráfico más rápido
hist(bio1_10, main="Distribución de frecuencias de temperatura")
#En base a todas las celdas
hist(bio1_10, maxpixels=ncell(bio1_10),main="Distribución de frecuencias de temperatura")

#Crear un raster "RasterStack" (agrupación de rasters) que contenga todas las variables
#Esto es útil para realizar una misma acción sobre todas las variables juntas (cómo graficar o recortar etc.)
Bio_10_Stack <- stack(bio1_10, bio2_10, bio3_10, bio4_10, bio5_10, bio6_10, bio7_10, bio8_10, bio9_10, bio10_10, bio11_10, bio12_10, bio13_10, bio14_10, bio15_10, bio16_10, bio17_10, bio18_10, bio19_10)

#Graficar todas las variables obtenidas de worldclim
plot(Bio_10_Stack)

#Las variables de WorldClim obtenidas son las "bio":
#BIO1 = Annual Mean Temperature
#BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
#BIO3 = Isothermality (BIO2/BIO7) (* 100)
#BIO4 = Temperature Seasonality (standard deviation *100)
#BIO5 = Max Temperature of Warmest Month
#BIO6 = Min Temperature of Coldest Month
#BIO7 = Temperature Annual Range (BIO5-BIO6)
#BIO8 = Mean Temperature of Wettest Quarter
#BIO9 = Mean Temperature of Driest Quarter
#BIO10 = Mean Temperature of Warmest Quarter
#BIO11 = Mean Temperature of Coldest Quarter
#BIO12 = Annual Precipitation
#BIO13 = Precipitation of Wettest Month
#BIO14 = Precipitation of Driest Month
#BIO15 = Precipitation Seasonality (Coefficient of Variation)
#BIO16 = Precipitation of Wettest Quarter
#BIO17 = Precipitation of Driest Quarter
#BIO18 = Precipitation of Warmest Quarter
#BIO19 = Precipitation of Coldest Quarter

#Si se quiere seleccionar de un "Raster Stack" solo los rasters de algunas variables, se utilizan los [[]]
#De esta forma se identifican las primeras 4 y se grafican
plot(Bio_10_Stack[[1:4]])

#Graficar la temperatura con colores de fríos a cálidos, creando una paleta de colores personalizada
tempcol <- colorRampPalette(c("purple", "blue", "skyblue", "green", "lightgreen", "yellow", "orange", "red", "darkred"))
#El número entre paréntesis después de "tempcol" indica cuantos colores "crear", usando los colores definidos en el objeto "tempcol"  como paleta. Entonces, tempcol(100) crea 100 colores formando una escala "suave" de colores.
plot(Bio_10_Stack[[1]], main="Annual Mean Temperature", col=tempcol(100))

#También se puede usar una paleta de colores de fríos a cálidos llamada "turbo" contenida en el paquete de R viridis
plot(Bio_10_Stack[[1]], main="Annual Mean Temperature", col=turbo(256))
#Si se utilizan pocos colores de la paleta (4):
plot(Bio_10_Stack[[1]], main="Annual Mean Temperature", col=turbo(4))

#Superponer la ubicación de las coordenadas de muestreo
plot(Bio_10_Stack[[1]], main="Annual Mean Temperature", col=turbo(256))
points(coords) 

#Ver la información sumaria del RasterStack
Bio_10_Stack 

#Parte 5: Obtener y extraer los datos de WorldClim para las coordenadas de muestreo

#Para esto se muestra como obtener los datos de worldclim con resolución de mucho detalle, 30 segundos (aprox 1 km2)
#Sin embargo, esta descarga puede demorar horas, en el caso de no requerir tanta resolución, se pueden descargar los datos a una resolución de 2.5m o 5m, como se explicó en la parte 4 del script
#El URL se obtuvo de la página web de WorldClim
W_url_30s <- "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_30s_bio.zip"
download.file(W_url_30s, destfile="wc2.1_30s_bio.zip") #Puede demorar horas
unzip("wc2.1_30s_bio.zip") #Al deszipear surge un error si no hay espacio en el disco. 

#Cargar las variables descargadas a R

bio1_30s = raster("wc2.1_30s_bio_1.tif")
bio2_30s = raster("wc2.1_30s_bio_2.tif")
bio3_30s = raster("wc2.1_30s_bio_3.tif")
bio4_30s = raster("wc2.1_30s_bio_4.tif")
bio5_30s = raster("wc2.1_30s_bio_5.tif")
bio6_30s = raster("wc2.1_30s_bio_6.tif")
bio7_30s = raster("wc2.1_30s_bio_7.tif")
bio8_30s = raster("wc2.1_30s_bio_8.tif")
bio9_30s = raster("wc2.1_30s_bio_9.tif")
bio10_30s = raster("wc2.1_30s_bio_10.tif")
bio11_30s = raster("wc2.1_30s_bio_11.tif")
bio12_30s = raster("wc2.1_30s_bio_12.tif")
bio13_30s = raster("wc2.1_30s_bio_13.tif")
bio14_30s = raster("wc2.1_30s_bio_14.tif")
bio15_30s = raster("wc2.1_30s_bio_15.tif")
bio16_30s = raster("wc2.1_30s_bio_16.tif")
bio17_30s = raster("wc2.1_30s_bio_17.tif")
bio18_30s = raster("wc2.1_30s_bio_18.tif")
bio19_30s = raster("wc2.1_30s_bio_19.tif")

#crear un RasterStack
Bio_30s_Stack <- stack(bio1_30s, bio2_30s, bio3_30s, bio4_30s, bio5_30s, bio6_30s, bio7_30s, bio8_30s, bio9_30s, bio10_30s, bio11_30s, bio12_30s, bio13_30s, bio14_30s, bio15_30s, bio16_30s, bio17_30s, bio18_30s, bio19_30s)

#grafico la temperatura a la resolución de 30s
plot(Bio_30s_Stack[[1]], main="Annual Mean Temperature")
points(coords)

#Grafico el RasterStack usando la paleta previamente creada de colores cálidos-fríos y lo exporto a PDF, TIFF, EPS
plot(Bio_30s_Stack[[1]], main="Annual Mean Temperature", col=tempcol(100))

#Es útil exportar los gráficos porque se pueden visualizar mejor que en la ventana de gráficos de RStudio
#Si se está utilizando RStudio se puede exportar el gráfico de esta manera
#En la ventana donde se observan los gráficos "Plots" se puede hacer click en "Export" y "Save as PDF" o "Save as image" (seleccionar TIFF o EPS)
#Alternativamente se pueden exportar utilizando comandos
#Por ejemplo, exportar en pdf:
pdf("Temperatura_30seg.pdf", width = 8, height = 5)
plot(Bio_30s_Stack[[1]], main="Annual Mean Temperature", col=tempcol(100))
dev.off()

#Extraer las variables bio de WorldClim para mis coordenadas de muestreo
valores <- extract(Bio_30s_Stack, coords)

#Armar un marco de datos con las coordenadas, los valores extraidos y la Abundancia
df <- cbind.data.frame(coordinates(coords),valores, abundancia=coords@data$Abundancia)

View(df)

#Exportar los datos como un archivo de texto
write.table(df, "df.txt", row.names = FALSE)

#Crear un objeto espacial con los resultados y exportarlo como un shapefile para ser leido en plataformas SIG
AbundF=SpatialPointsDataFrame(coords,df) #crear el objeto espacial que contiene los valores extraidos
head(AbundF)
summary(AbundF)

writeOGR(AbundF,".", "Abund_bio", driver="ESRI Shapefile")  #exportar como un shapefile de puntos


#Parte 6: Graficar un mapa de los sitios de muestreo, utilizando como fondo los valores obtenidos de WorldClim

#Graficar para que quede una visualización
#Utilizar la función getData, para obtener el mapa de la Argentina y sus provincias
Argentina1 <- getData('GADM' , country="ARG", level=1) 
plot(Argentina1, axes=TRUE)

#Superponer las sitios de muestreo para identificar qué provincias seleccionar para el gráfico
plot(coords, pch=20, col="black", cex=1, add=TRUE)

#Observar la tabla de atributos del mapa de la Argentina e identificar que el nombre de la variable que contiene las provincias es NAME_1 y observar que las provincias están  identificadas como 'Río Negro' y 'Neuquén' 
View(Argentina1@data)

#Obtener un objeto espacial que solo contenga a las provincias de Río Negro y Neuquén
s <- subset(Argentina1, NAME_1 %in% c('Río Negro', 'Neuquén'))
plot(s, axes=TRUE)

#Cortar un rectángulo de la variable de temperatura de WorldClim de la extensión de mis provincias de interés
rect <- crop(Bio_30s_Stack[[1]], s) 
plot (rect)

#Hacer una máscara, para observar solo los valores de temperatura  que están contenidos en las provincias
mascara <- mask(rect, s)
plot(mascara, main="Sitios de muestreo y temperatura media anual")
plot(coords, pch=20, col="black", cex=1, add=TRUE)
plot(s, add=TRUE)
dev.off()  



#Parte 7: estadística elemental. Exploración de los datos de temperatura extraídos para los sitios de muestreo y la relación entre variables

cor(df)  #correlaciones entre todas las variables

dev.off() #cerrar cualquier venta gráfica que esté abierta

#Gráfico de cajas, variación de la temperatura media anual en los sitios de muestreo
boxplot(df$wc2.1_30s_bio_1)
#En el comando  anterior se mostró como hacer un gráfico de la manera más sencilla, solo poniendo la variable de estudio como argumento
#En el resto de los gráficos se agregarán comandos para mostrar como incorporar los títulos de los ejes y título del gráfico, pero estos no son necesarios para hacer una versión elemental del gráfico 
boxplot(df$wc2.1_30s_bio_1,  xlab="Mis sitios", ylab=expression(paste("Temperatura media anual (",degree,"C)")))  

#Histograma, frecuencia de los valores de temperatura media anual en los sitios de muestreo
hist(df$wc2.1_30s_bio_1, main="", ylab= "Número de sitios", xlab=expression(paste("Temperatura media anual (",degree,"C)")))  

#Gráfico de dispersión, relación entre la abundancia de hormigas y la temperatura
plot(df$wc2.1_30s_bio_1, df$abundancia,  ylab= "Abundancia de hormigas", xlab=expression(paste("Temperatura media anual (",degree,"C)"))) 

#Gráfico de dispersión, con la abundancia en escala logaritmica
plot(df$wc2.1_30s_bio_1, df$abundancia, log='y', ylab= "Abundancia de hormigas", xlab=expression(paste("Temperatura media anual (",degree,"C)"))) 


#Regresión lineal simple con distribución normal (logaritmo de la abundancia de hormigas como variable respuesta y temperatura media anual como variable predictora)
reg= lm(log10(df$abundancia) ~ df$wc2.1_30s_bio_1) 
plot(df$wc2.1_30s_bio_1, df$abundancia, log='y', ylab= "Abundancia de hormigas", xlab=expression(paste("Temperatura media anual (",degree,"C)"))) 
abline(reg, col="blue") #agregar la línea de tendencia
summary(reg)

#FIN#