Introducción al ambiente R

Esta página muestra algunos ejercicios básicos para empezar a familiarizarse con el ambiente de trabajo R. Las funciones que veremos a continuación son las que todo principiante en R debería conocer y tratar de memorizar puesto que le servirán por siempre. A medida que el estudiante se familiarice con estas funciones podrá continuar con funciones y códigos un poco más complejos. Tengan paciencia y mucha persistencia.

Primero vamos a descargar los datos de este enlace (solo estará disponible estos dos días del taller) https://drive.google.com/drive/folders/1za3kjoyaOev-44yvafRJIRSle4zI_sUx?usp=sharing

Vamos a instalar los paquetes que trae R y que usaremos. # install.packages(“sp”, “raster”, “rgdal”, “dismo”, “envirem”)

Direccionar un directorio-folder de trabajo

setwd("/home2/Dropbox/Curso/TallerPINCC")

Algunas operaciones básicas

a <- 100
b <- 1:10
c <- "Curso"
d <- c("Curso", "UNAM", "Clima")

a
## [1] 100
b
##  [1]  1  2  3  4  5  6  7  8  9 10
c
## [1] "Curso"
d
## [1] "Curso" "UNAM"  "Clima"
# crear una matriz
e <- matrix(1:9, nrow = 3, ncol = 3)
e
##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]    3    6    9

Ahora vamos a jugar con datos. Subiremos un conjunto de datos en formato CSV (comma-separated-value). Estos archivos se pueden abrir con Excel o LibreOffice.

mis_datos <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/distancias_islas_caribe.csv", header = TRUE)

Funciones básicas para empezar a explorar

# head me muestra el encabezado del archivo que acabo de subir (las primeras filas)
head(mis_datos)
##      ID                   Island species distance
## 1 11342             Not assigned  aeneus    6.075
## 2 11359 Southern Lesser Antilles  aeneus    8.194
## 3 11360             Not assigned  aeneus    4.618
## 4 11361             Not assigned  aeneus    4.382
## 5 11362             Not assigned  aeneus    4.611
## 6 11363             Not assigned  aeneus    2.855
# tail me muestra la parte final del archivo
tail(mis_datos)
##            ID                   Island   species distance
## 129312 919813 Southern Lesser Antilles whitemani    9.946
## 129313 919816 Southern Lesser Antilles whitemani     8.92
## 129314 919817 Southern Lesser Antilles whitemani    8.755
## 129315 919818 Southern Lesser Antilles whitemani    8.831
## 129316 919831 Southern Lesser Antilles whitemani     8.73
## 129317 919833             Not assigned whitemani    8.848
# class: indica que tipo de archivo estoy usando, e.g., un vector, una lista, un dataframe, una matriz, etc.
class(mis_datos)
## [1] "data.frame"
# str: me muestra la estructura de los datos que estoy visualizando. Es útil para ver donde esta algo que luego necesito extraer o indexar.
str(mis_datos)
## 'data.frame':    129317 obs. of  4 variables:
##  $ ID      : int  11342 11359 11360 11361 11362 11363 11364 11366 11375 11377 ...
##  $ Island  : chr  "Not assigned" "Southern Lesser Antilles" "Not assigned" "Not assigned" ...
##  $ species : chr  "aeneus" "aeneus" "aeneus" "aeneus" ...
##  $ distance: chr  "6.075" "8.194" "4.618" "4.382" ...
# dim: permite ver las dimensiones de mi objeto ya sea un dataframe o una matriz.
dim(mis_datos)
## [1] 129317      4
# length: es igual que el anterior pero solo para objetos tipo lista
length(mis_datos)
## [1] 4
#names: permite ver los nombres  de las columnas
names(mis_datos)
## [1] "ID"       "Island"   "species"  "distance"

Ahora carguemos el mismo archivo pero que cada fila tenga un nombre único usando la primera columna (ID).

mis_datos <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/distancias_islas_caribe.csv", header=TRUE, row.names=1)

Ahora carguemos otro archivo con datos para unas especies

datos <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/areas_nichos.csv", header = TRUE, row.names = 1)

Hasta ahora hemos aprendido a cargar archivos en formato CSV y a visualizar sus dimensiones y estructura. Ahora pasaremos a una de las cosas más importante de trabajar con R que es la INDEXACIÓN. La indexación es clave para aprender a dominar R.

Primero, vamos a generar un vector númerico a partir del segundo archivo que subimos (datos)

var1 <- datos[,1]
head(var1)
## [1]  5.373085  4.209303 12.610739  6.441746  2.915383  9.306049
dim(var1)
## NULL
class(var1)
## [1] "numeric"
length(var1)
## [1] 131

Queremos asignarle los nombres de las especies a cada valor del vector númerico

names(var1)<-row.names(datos)

Ahora vamos a indexar cada una de las variables del primer archivo (mis_datos):

# voy  a indexar la columna 2 y 3. 
variables1 <- mis_datos[,c(2,3)] 

# indexar la columna 1 y 3
variables2 <- mis_datos[,c(1,3)]

# indexar de la columna 2 a la 3 de forma continua.
variables3 <- mis_datos[,c(2:3)]

Voy a convertir el objeto lista en un dataframe

dataframe1 <- as.data.frame(variables1)
class(dataframe1)
## [1] "data.frame"

Ahora de dataframe a lista

data.list1 <- as.list(variables1)
class(data.list1)
## [1] "list"

De dataframe a matriz

vars.matrix <- as.matrix(variables1)
vars.matrix[1:100] # voy a visualizar los primeros 100 datos de esta matriz.
##   [1] "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"  
##   [7] "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"  
##  [13] "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"  
##  [19] "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"   "aeneus"  
##  [25] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [31] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [37] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [43] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [49] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [55] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [61] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [67] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "ahli"    
##  [73] "ahli"     "ahli"     "ahli"     "ahli"     "ahli"     "aliniger"
##  [79] "aliniger" "aliniger" "aliniger" "aliniger" "aliniger" "aliniger"
##  [85] "aliniger" "aliniger" "aliniger" "aliniger" "aliniger" "aliniger"
##  [91] "aliniger" "aliniger" "aliniger" "aliniger" "aliniger" "aliniger"
##  [97] "aliniger" "aliniger" "aliniger" "aliniger"

Ahora empezaremos a jugar con las tablas. Vamos a UNIR tres tablas diferentes con base en una columna en común.

species <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/datos_especies_islas.csv")
areas <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/areas_islas.csv")
distancias <- read.csv("/home2/Dropbox/Curso/TallerPINCC/Datos/Tablas/distancias.csv")

Concatenar dos tablas - dos dataframes-

species.areas <- merge(species, areas, by="Islands")

Concatenar esta tabla nueva con la de distancias

spp.areas.dist <- merge(species.areas, distancias, by="species")

Ahora vamos a indexar por subconjuntos (subsets)

dactyloa <- subset(spp.areas.dist, clades=="Dactyloa")
homolechis <- subset(spp.areas.dist, species=="homolechis")

Miremos como podemos indexar una columna de dos formas diferentes

dist.homo1 <- homolechis[,7]
dist.homo2 <- homolechis$distancias

Para evitar que se convierta en factor usar lo siguiente:

dist.homo <- as.numeric(as.vector(homolechis$distancias))

Algunas estadísticas y gráficas básicas en R

# Histograma de frecuencias
hist(dist.homo)

mean(dist.homo)
## [1] 3.917335
min(dist.homo)
## [1] 1.233
max(dist.homo)
## [1] 29.516
sd(dist.homo)
## [1] 1.909522
var(dist.homo)
## [1] 3.646274

Vamos a indexar el promedio

a <- mean(dist.homo)

Guardar y cargar un objeto de R (RData)

save.image("misDatos.RData")

load("misDatos.RData")

Gráficas sencillas y transformación de datos

hist(dist.homo, main="histograma",xlab="distancias", ylab="frecuencias")

Transformación logarítmica

log.dist <- log(dist.homo)
hist(log.dist, main="histograma", xlab="distancias")

Correlaciones y regresión lineal

cor(datos$nicho, datos$area)
## [1] 0.09052859
# primero x y luego y
plot(datos$nicho, datos$area)

Modelo de regresión lineal

reg.l <- lm(datos$nicho~datos$area)

Hacer un histograma de los residuales Usaré la indexación de los residuales del objeto reg.l con el simbolo de $

hist(reg.l$residuals)

Gŕafica de las dos variables

plot(datos$area, datos$nicho)

Voy a crear una columna nueva que contenga los valores de la variable “nicho” pero transformados con log10

datos["logNicho"] <- log(datos$nicho)
plot(datos$area, datos$logNicho)

Vuelvo a hacer el modelo de regresión lineal, pero con la nueva variable

reg.l2 <- lm(datos$logNicho~datos$area)
reg.l2
## 
## Call:
## lm(formula = datos$logNicho ~ datos$area)
## 
## Coefficients:
## (Intercept)   datos$area  
##     1.75181      0.03249

Extraigo los valores de la línea de tendencia (intercepto y la pendiente). Uso str para saber donde están.

str(reg.l2)
## List of 12
##  $ coefficients : Named num [1:2] 1.7518 0.0325
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "datos$area"
##  $ residuals    : Named num [1:131] -0.2301 -0.4484 0.6436 -0.0346 -0.8426 ...
##   ..- attr(*, "names")= chr [1:131] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:131] -21.1738 -0.4626 0.6935 0.0209 -0.7738 ...
##   ..- attr(*, "names")= chr [1:131] "(Intercept)" "datos$area" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:131] 1.91 1.89 1.89 1.9 1.91 ...
##   ..- attr(*, "names")= chr [1:131] "1" "2" "3" "4" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:131, 1:2] -11.4455 0.0874 0.0874 0.0874 0.0874 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:131] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "datos$area"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.09 1.07
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 129
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = datos$logNicho ~ datos$area)
##  $ terms        :Classes 'terms', 'formula'  language datos$logNicho ~ datos$area
##   .. ..- attr(*, "variables")= language list(datos$logNicho, datos$area)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "datos$logNicho" "datos$area"
##   .. .. .. ..$ : chr "datos$area"
##   .. ..- attr(*, "term.labels")= chr "datos$area"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(datos$logNicho, datos$area)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "datos$logNicho" "datos$area"
##  $ model        :'data.frame':   131 obs. of  2 variables:
##   ..$ datos$logNicho: num [1:131] 1.68 1.44 2.53 1.86 1.07 ...
##   ..$ datos$area    : num [1:131] 4.92 4.12 4.28 4.48 4.95 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language datos$logNicho ~ datos$area
##   .. .. ..- attr(*, "variables")= language list(datos$logNicho, datos$area)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "datos$logNicho" "datos$area"
##   .. .. .. .. ..$ : chr "datos$area"
##   .. .. ..- attr(*, "term.labels")= chr "datos$area"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(datos$logNicho, datos$area)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "datos$logNicho" "datos$area"
##  - attr(*, "class")= chr "lm"

Los voy a indexar…

intercepto <- reg.l2$coefficients[1] # intercepto en Y
pendiente <- reg.l2$coefficients[2] # pendiente

Hago la gráfica con la línea de ajuste

plot(datos$area, datos$logNicho)
abline(a=intercepto, b=pendiente)

La misma gráfica pero le pongo las etiquetas a los ejes

plot(datos$area, datos$logNicho, xlab="Area", ylab="Nicho")
abline(a=intercepto, b=pendiente)

EJERCICIO # 1 (van a usar lo aprendido)

Hagamos lo mismo, pero con un conjunto de datos más grandes.

datos_grandes <- spp.areas.dist
datos_grandes["log_Nichos"]<- NA

Extraigo los datos de distancias

nichos <- as.numeric(as.vector(datos_grandes$distancias))
## Warning: NAs introduced by coercion

Transformación logarítmica

log_nichos<-log(nichos)

Ahora se los pego a mis datos

datos_grandes$log_Nichos <- log_nichos

TAREA:

Ahora cada uno debe hacer la regresión y gráficar con la línea de ajuste del modelo.

Cartografía básica y operaciones con datos espaciales y rasters

Aquí empezaremos con algunos ejercicios para trabajar con objetos espaciales y hacer algunos mapas básicos. Vamos a generar unos datos aleatorios de coordenadas geográficas y de precipitación.

name <- LETTERS[1:10]
longitud <- c(-116.7, -120.4, -116.7, -113.5, -115.5,
-120.8, -119.5, -113.7, -113.7, -110.7)
latitud <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9,
36.2, 39, 41.6, 36.9)
estaciones <- cbind(longitud, latitud)

# Simular datos de precipitación
set.seed(0)
# función runif: sirve para generar valores aleatorios siguiendo una distribución uniforme
precip <- (runif(length(latitud))*10)^3

Ahora vamos a gráficar estos datos y el tamaño de los puntos será proporcional a la cantidad de precipitación.

psize <- 1 + precip/500
plot(estaciones, cex=psize, pch=20, col='red', main='Precipitación')

# adicionarle los nombres a la gráfica
text(estaciones, name, pos=4)

# adicionar una leyenda
breaks <- c(100, 500, 1000, 2000)
legend("topright", legend=breaks, pch=20, pt.cex=psize, col='red', bg='gray')

Puntos espaciales:

Vamos a generar vectores con datos de longitud y latitud, y luego los combinamos.

longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)

Vamos a usar el paquete “sp”. D Debemos instalarlo usando la función install.packages(“sp”) y luego cargarlo con library(sp)

library(sp)
pts <- SpatialPoints(lonlat)

¿Qué clase de objeto es?

class(pts)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

Exploramos que tiene

showDefault(pts)
## An object of class "SpatialPoints"
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments: NA

El sistema de coordenadas de referencia de nuestro objeto espacial no está definido. Vamos a definirlo con la ayuda de la función CRS

crdref <- CRS('+proj=longlat +datum=WGS84')
pts <- SpatialPoints(lonlat, proj4string=crdref)

Cargamos el paquete raster

library(raster)
pts
## class       : SpatialPoints 
## features    : 10 
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

Vamos a generar (inventar) unos datos de precipitación asociados a esas coordenadas. Estos datos NO son reales, es solo para mostrar como podemos visualizarlos como objetos espaciales.

df <- data.frame(ID=1:nrow(lonlat), precip=(latitude-30)^3)

Combinamos los datos espaciales y los datos de precipitación

ptsdf <- SpatialPointsDataFrame(pts, data=df)
ptsdf
## class       : SpatialPointsDataFrame 
## features    : 10 
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 2
## names       : ID,   precip 
## min values  :  1,  185.193 
## max values  : 10, 3581.577

Vemos que tiene ese objeto espacial

str(ptsdf)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  10 obs. of  2 variables:
##   .. ..$ ID    : int [1:10] 1 2 3 4 5 6 7 8 9 10
##   .. ..$ precip: num [1:10] 3582 2000 705 1772 185 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:10, 1:2] -117 -120 -117 -114 -116 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "longitude" "latitude"
##   ..@ bbox       : num [1:2, 1:2] -120.8 35.7 -110.7 45.3
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "longitude" "latitude"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs"
##   .. .. ..$ comment: chr "GEOGCRS[\"unknown\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.25722"| __truncated__

También podemos verlos así:

showDefault(ptsdf)
## An object of class "SpatialPointsDataFrame"
## Slot "data":
##    ID   precip
## 1   1 3581.577
## 2   2 2000.376
## 3   3  704.969
## 4   4 1771.561
## 5   5  185.193
## 6   6  704.969
## 7   7  238.328
## 8   8  729.000
## 9   9 1560.896
## 10 10  328.509
## 
## Slot "coords.nrs":
## numeric(0)
## 
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

Líneas y polígonos espaciales. Hacer puntos espaciales es bastante fácil. Hacer líneas y polígonos es un poco más complicado pero vamos a usar las funciones SpatialLines y SpatialPolygons para esto.

lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(lon, lat)
lns <- spLines(lonlat, crs=crdref)
lns
## class       : SpatialLines 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

Graficamos nuestro objeto de líneas

plot(lns)

Creamos un shapefile (polígono) con los puntos

pols <- spPolygons(lonlat, crs=crdref)
pols
## class       : SpatialPolygons 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

Graficamos el polígono

plot(pols, axes=TRUE, las=1)
plot(pols, border='blue', col='yellow', lwd=3, add=TRUE)
points(pts, col='red', pch=20, cex=3)

Ahora vamos a manipular algunos archivos tipo raster

El paquete sp permite trabajar con archivos de clase SpatialGridDataFrame y SpatialPixelsDataFrame. Estos archivos son archivos que usamos muy a menudo en operaciones espaciales y SIG. De aquí en adelante nos vamos a empezar a usar el paqute raster para trabajar con archivo tipo raster (celdas). Con este paquete podemos hacer operaciones matemáticas con rasters y un sinnúmero de otros procesos que se usan en programas dedicadas a SIG (e.g., ArcGIS, QGIS, etc.).

¿Qué es un capa de tipo raster? Un RasterLayer es un objeto de R que contiene datos de una sola variables. Como vimos en la clase este objeto debe tener información del número de celdas, hileras (rows) y columnas (columns), la extensión geográfica o espacial (x, y - min-, -max-) y lo más importante en SIG: el sistema de coordenadas de referencia.

library(raster)
r <- raster(ncol=10, nrow=10, xmx=-80, xmn=-150, ymn=20, ymx=60)
r
## class      : RasterLayer 
## dimensions : 10, 10, 100  (nrow, ncol, ncell)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs
values(r) <- runif(ncell(r))
r
## class      : RasterLayer 
## dimensions : 10, 10, 100  (nrow, ncol, ncell)
## resolution : 7, 4  (x, y)
## extent     : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.01339033, 0.9926841  (min, max)

Ahora podemos gráficar y ver como sale nuestro raster

plot(r)

Ahora vamos a subir unas capas climáticas que tenemos guardadas en nuestro computador.

capas <- list.files(path="/home2/Dropbox/Curso/TallerPINCC/Datos/Bioclimaticos/", pattern=".tif", full.names = T)
capas <- stack(capas)
capas
## class      : RasterStack 
## dimensions : 180, 360, 64800, 19  (nrow, ncol, ncell, nlayers)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : NA 
## names      :      bio_01,      bio_02,      bio_03,      bio_04,      bio_05,      bio_06,      bio_07,      bio_08,      bio_09,      bio_10,      bio_11,      bio_12,      bio_13,      bio_14,      bio_15, ... 
## min values : -234.128322,    9.500000,    9.396564,   68.000000,  -55.126503, -525.923628,   53.489130, -237.702800, -432.567114,  -91.281594, -471.184995,    0.000000,    0.000000,    0.000000,    0.000000, ... 
## max values :    313.0625,    195.2251,     95.0000,  22549.9589,    484.3726,    258.0000,    722.3212,    369.2097,    355.0000,    375.4345,    288.8571,   6995.0284,   1482.7247,    345.1818,    215.9236, ...

Ahora vamos a graficar una sola capa climática (bio1)

plot(capas[[1]])

Operaciones básicas con rasters

promedio_bio1 <- cellStats(capas[[1]], 'mean')
promedio_bio12 <- cellStats(capas[[12]], 'mean')

Vamos a jugar un poco con shapefiles y rasters

library(rgdal) 
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
ame1 <- readOGR(dsn="/home2/Dropbox/Curso/TallerPINCC/Datos/Shapefile", layer="cont_america")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home2/Dropbox/Curso/TallerPINCC/Datos/Shapefile", layer: "cont_america"
## with 28 features
## It has 1 fields
ame2 <- readOGR("/home2/Dropbox/Curso/TallerPINCC/Datos/Shapefile/cont_america.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home2/Dropbox/Curso/TallerPINCC/Datos/Shapefile/cont_america.shp", layer: "cont_america"
## with 28 features
## It has 1 fields
plot(ame1)
plot(ame2)

Vamos a indexar un país y luego cortar las capas climáticas para ese país

head(ame1)
##                   PAÍS
## 0 Aruba (Paises Bajos)
## 1            Argentina
## 2               Belice
## 3              Bolivia
## 4               Brasil
## 5               Canadá
mx <- subset(ame1, PAÍS=="México")
plot(mx)

Cortamos las capas usando la función mask. También podemos usar la función crop

# Usando mask
mx.clima <- raster::mask(capas, mx)
plot(mx.clima[[1]])

# Usando crop
mx.clima.v2 <- raster::crop(capas, mx)
plot(mx.clima.v2[[1]])

Lo ideal es primero crop y luego mask

mx.clima.v3 <- crop(capas, mx)
mx.clima.v3 <- mask(mx.clima.v3, mx)
plot(mx.clima.v3[[1]])

Agregarle el mapa de México como contorno

plot(mx.clima.v3[[1]])
plot(mx, add=T)

¿Cómo se calculan los bioclimáticos?

Hoy vamos a hacer algunos ejercicios para trabajar con datos bioclimáticos de algunas bases de datos. Se instalarán algunos paquetes necesarios de CRAN y GitHub.

library(raster)
library(dismo)
library(envirem)
## Loading required package: palinsol
## Loading required package: gsl
## 
## Please see the vignette at https://envirem.github.io/ENVIREM_tutorial.html for a detailed walk-through of this package.

Descargae de datos directamente de Worldclim

https://www.worldclim.org/data/worldclim21.html

# worldclim (resolución 10 minutos)
clim <- getData("worldclim",var="bio",res=10) 
clim
## class      : RasterStack 
## dimensions : 900, 2160, 1944000, 19  (nrow, ncol, ncell, nlayers)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :  bio1,  bio2,  bio3,  bio4,  bio5,  bio6,  bio7,  bio8,  bio9, bio10, bio11, bio12, bio13, bio14, bio15, ... 
## min values :  -269,     9,     8,    72,   -59,  -547,    53,  -251,  -450,   -97,  -488,     0,     0,     0,     0, ... 
## max values :   314,   211,    95, 22673,   489,   258,   725,   375,   364,   380,   289,  9916,  2088,   652,   261, ...
names(clim)
##  [1] "bio1"  "bio2"  "bio3"  "bio4"  "bio5"  "bio6"  "bio7"  "bio8"  "bio9" 
## [10] "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16" "bio17" "bio18"
## [19] "bio19"

Extraigo solo la temperatura promedio anual y la precipitación anual

tpa <- clim[[1]]
pa <- clim[[12]]

Graficar

par(mfrow=c(2,1))
plot(tpa, main="Temperatura promedio anual") 
plot(pa, main="Precipitación anual")

tpa
## class      : RasterLayer 
## dimensions : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio1.bil 
## names      : bio1 
## values     : -269, 314  (min, max)
pa
## class      : RasterLayer 
## dimensions : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : bio12.bil 
## names      : bio12 
## values     : 0, 9916  (min, max)

Estadísticas a través de celdas (media, máximo, mínimo y desviación estándar)

mean.tpa <- cellStats(tpa, stat='mean')
mean.pa <- cellStats(pa, stat='mean')

max.tpa <- cellStats(tpa, stat='max')
max.pa <- cellStats(pa, stat='max')

min.tpa <- cellStats(tpa, stat='min')
min.pa <- cellStats(pa, stat='min')

sd.tpa <- cellStats(tpa, stat='sd')
sd.pa <- cellStats(pa, stat='sd')

Vamos a subir los datos de temperaturas y precipitación que bajaron con Oscar. Borrar primero el raster del 2019

# primero temperatura mínima mensual
tmin <- stack(list.files(path="/home2/Dropbox/Curso/TallerPINCC/Datos/WorldClim_Historical/WC_Tmin/", pattern='*.nc', full.names=TRUE))
## Loading required namespace: ncdf4
# luego temperatura máxima mensual
tmax <- stack(list.files(path="/home2/Dropbox/Curso/TallerPINCC/Datos/WorldClim_Historical/WC_Tmax/", pattern='*.nc', full.names=TRUE))

# precipitación mensual
prec <- stack(list.files(path="/home2/Dropbox/Curso/TallerPINCC/Datos/WorldClim_Historical/WC_Prec/", pattern='*.nc', full.names=TRUE))

Calcular los promedios mensuales para todo el periodo (1961-2018). 696 rasters son 58 años de datos.

t <- 696/12

# Separar los años
t1 <- split(seq_len(696), rep(1:58, each=12))

# Agrupar cada mes de todos los años
t2 <- split(seq_len(696), rep(1:12, length.out=696))

t.tmin <- lapply(1:12, function(x) tmin[[t2[[x]]]])
t.tmax <- lapply(1:12, function(x) tmax[[t2[[x]]]])
t.prec <- lapply(1:12, function(x) prec[[t2[[x]]]])

# calcular el promedio de cada mes de todos los años
t.tmin.mean <- lapply(1:12, function(x) mean(t.tmin[[x]]))
t.tmax.mean <- lapply(1:12, function(x) mean(t.tmax[[x]]))
t.prec.mean <- lapply(1:12, function(x) mean(t.prec[[x]]))

# convierto la lista de rasters en un stack
t.tmin.mean2 <- stack(t.tmin.mean)
t.tmax.mean2 <- stack(t.tmin.mean)
t.prec.mean2 <- stack(t.prec.mean)

names(t.tmin.mean2) <- c("ene","feb","mar","abr","may","jun","jul","ago","sep","oct","nov","dic")
names(t.tmax.mean2) <- c("ene","feb","mar","abr","may","jun","jul","ago","sep","oct","nov","dic")
names(t.prec.mean2) <- c("ene","feb","mar","abr","may","jun","jul","ago","sep","oct","nov","dic")

# indexado manual
# tmin.ene <- tmin[[c(1,13,25,37,49,61,73,85,97,109,121,133,...,)]]

Ahora vamos a generar los bioclimáticos (19) para los 58 años

library(dismo)
bios1961_2018 <- biovars(t.prec.mean2, t.tmin.mean2, t.tmax.mean2)

# generar el mapa para la bio1
plot(bios1961_2018[[1]])

# Generar el mapa para la bio12
plot(bios1961_2018[[12]])

Calcular las diferencias entre los bioclimáticos generados y los descargados en la primera parte.

# remuestrear al pixel más grande
clim2 <- raster::resample(clim, bios1961_2018[[1]], method="bilinear")

# diferencia
dif.tpa <- bios1961_2018[[1]] - clim2[[1]]/10
plot(dif.tpa)

# histograma

dif.tpa2 <- rasterToPoints(dif.tpa)
hist(dif.tpa2[,3])

Esto es todo por ahora!

Fuentes: Elaboración propia (Julián A. Velasco) y con rutinas adaptadas de: http://www.rspatial.org/intr/index.html (R. Hijmans, 2016)