Introducción al ambiente R

En esta página se muestran 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: https://drive.google.com/drive/folders/1UOp6HdmGZuZCQAH8G8BqJmuOEW1hTWp0?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("/mnt/tb/DD/Dropbox/Curso_PEMBU")

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("./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("./Tablas/distancias_islas_caribe.csv", header=TRUE, row.names=1)

Ahora carguemos otro archivo con datos para unas especies

datos <- read.csv("./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.

especies <- read.csv("./Tablas/datos_especies_islas.csv")
areas <- read.csv("./Tablas/areas_islas.csv")
distancias <- read.csv("./Tablas/distancias.csv")

Concatenar dos tablas - dos dataframes-

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

Concatenar esta tabla nueva con la de distancias

spp.areas.dist <- merge(especies.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_Temp"]<- NA

Extraigo los datos de distancias

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

Transformación logarítmica

log_temp<-log(nichos)

Ahora se los pego a mis datos

datos_grandes$log_Temp <- log_temp

TAREA:

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

Análisis de datos meteorológicos

Existen varias librerias de R para el análisis de datos meteorológicos. Vamos a usar una muy común que se llama “climate”. Este paquete permite automatizar la descarga de datos de varios repositorios disponibles públicamente.

OGIMET (ogimet.com) University of Wyoming - atmospheric vertical profiling data (http://weather.uwyo.edu/upperair/).

National Oceanic & Atmospheric Administration - Earth System Research Laboratory - Global Monitoring Division (NOAA)

# install.packages(c("climate", "dplyr", "tidyr))
library(climate)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)

ns <- nearest_stations_ogimet(country ="Mexico",
                             point = c(-4, 56),
                             no_of_stations = 100, 
                             add_map = TRUE)
## [1] "http://ogimet.com/cgi-bin/gsynres?lang=en&state=Mexico&osum=no&fmt=html&ord=REV&ano=2022&mes=06&day=27&hora=06&ndays=1&Send=send"
## /tmp/RtmpBpIN9C/file127c66e3cf16

# Vemos el encabezado

head(ns)
##    wmo_id                    station_names       lon      lat alt distance [km]
## 43  76647                Valladolid, Yuc.  -88.21667 20.70001  22      10245.31
## 56  76698 Felipe Carrillo Puerto, Q. Roo.  -88.03333 19.58335  20      10275.51
## 61  76750                Chetumal, Q. Roo  -88.30001 18.48334   9      10352.53
## 36  76593                  Progreso, Yuc.  -89.65001 21.26667   2      10369.75
## 42  76644   Merida / lic Manuel Crecencio  -89.65001 20.93335  10      10383.86
## 55  76695                 Campeche, Camp.  -90.50001 19.83335   5      10519.17
PL = stations_ogimet(country = "Mexico", add_map = TRUE)
## [1] "http://ogimet.com/cgi-bin/gsynres?lang=en&state=Mexico&osum=no&fmt=html&ord=REV&ano=2022&mes=06&day=27&hora=06&ndays=1&Send=send"
## /tmp/RtmpBpIN9C/file127c67716d3fb

Descargar datos por hora en un año aleatorio de Ogimet para la Ciudad de México. Primero buscar el ID de la estación https://www1.ncdc.noaa.gov/pub/data/noaa/isd-history.txt

# Descarga de datos del servicio de la NOAA.

df = meteo_noaa_hourly(station = "766800-99999", 
                       year = sample(2000:2020, 1))
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2017/766800-99999-2017.gz"
## /tmp/RtmpBpIN9C/file127c659a5ae82
head(df)
##                  date year month day hour     lon  lat  alt  t2m dpt2m ws  wd
## 1 2017-01-01 00:00:00 2017     1   1    0 -99.183 19.4 2303 18.8  11.1  3 230
## 2 2017-01-01 03:00:00 2017     1   1    3 -99.183 19.4 2303 15.9  10.4  2  50
## 3 2017-01-01 09:00:00 2017     1   1    9 -99.183 19.4 2303 12.8   9.6  1 140
## 4 2017-01-01 12:00:00 2017     1   1   12 -99.183 19.4 2303 11.2   8.5  1 180
## 5 2017-01-01 18:00:00 2017     1   1   18 -99.183 19.4 2303 18.5  10.9  1  80
## 6 2017-01-01 21:00:00 2017     1   1   21 -99.183 19.4 2303 23.2  13.0  3 180
##      slp visibility
## 1 1008.7       5000
## 2 1009.5       4000
## 3 1011.1       2000
## 4 1013.0       1000
## 5 1015.0       1000
## 6 1010.0       2000
hist(df$t2m, main="CDMX 2020", xlab="Temperatura horaria", ylab="Frecuencia")

# Exploramos cómo están los datos 
sapply(df, class)
## $date
## [1] "POSIXct" "POSIXt" 
## 
## $year
## [1] "integer"
## 
## $month
## [1] "integer"
## 
## $day
## [1] "integer"
## 
## $hour
## [1] "integer"
## 
## $lon
## [1] "numeric"
## 
## $lat
## [1] "numeric"
## 
## $alt
## [1] "integer"
## 
## $t2m
## [1] "numeric"
## 
## $dpt2m
## [1] "numeric"
## 
## $ws
## [1] "numeric"
## 
## $wd
## [1] "integer"
## 
## $slp
## [1] "numeric"
## 
## $visibility
## [1] "integer"
# Indexamos la variable de temperatura y excluimos los datos faltantes (NAs)
temp <- na.omit(df$t2m)

hist(temp)

# Sacar la temperatura máxima horaria para el 2020
max(temp)
## [1] 32.4
# Temperatura mínima

min(temp)
## [1] 0.3
plot(temp)

# Ahora vamos a descargar 10 años de datos para la misma estación

df2 = meteo_noaa_hourly(station = "766800-99999", 
                       year = c(2010:2020))
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2010/766800-99999-2010.gz"
## /tmp/RtmpBpIN9C/file127c63131c708
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2011/766800-99999-2011.gz"
## /tmp/RtmpBpIN9C/file127c6d2f43c5
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2012/766800-99999-2012.gz"
## /tmp/RtmpBpIN9C/file127c66763536
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2013/766800-99999-2013.gz"
## /tmp/RtmpBpIN9C/file127c668c303ee
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2014/766800-99999-2014.gz"
## /tmp/RtmpBpIN9C/file127c62d2abc74
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2015/766800-99999-2015.gz"
## /tmp/RtmpBpIN9C/file127c6156fa9ce
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2016/766800-99999-2016.gz"
## /tmp/RtmpBpIN9C/file127c644b9d09e
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2017/766800-99999-2017.gz"
## /tmp/RtmpBpIN9C/file127c695181f8
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2018/766800-99999-2018.gz"
## /tmp/RtmpBpIN9C/file127c65ae98d0a
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2019/766800-99999-2019.gz"
## /tmp/RtmpBpIN9C/file127c64a41c758
## [1] "https://www1.ncdc.noaa.gov/pub/data/noaa/2020/766800-99999-2020.gz"
## /tmp/RtmpBpIN9C/file127c6491c501
head(df2)
##                  date year month day hour     lon  lat  alt  t2m dpt2m ws  wd
## 1 2010-01-01 00:00:00 2010     1   1    0 -99.183 19.4 2303 20.3   8.9  1 140
## 2 2010-01-01 03:00:00 2010     1   1    3 -99.183 19.4 2303 17.0   9.7  2  50
## 3 2010-01-02 06:00:00 2010     1   2    6 -99.183 19.4 2303 10.6  10.3  1 360
## 4 2010-01-02 09:00:00 2010     1   2    9 -99.183 19.4 2303  9.4   9.1  1  90
## 5 2010-01-02 12:00:00 2010     1   2   12 -99.183 19.4 2303  9.0   8.7  1 180
## 6 2010-01-04 06:00:00 2010     1   4    6 -99.183 19.4 2303 10.0   9.7  1 360
##      slp visibility
## 1 1008.6       8000
## 2 1011.9       8000
## 3 1012.6       8000
## 4 1012.5       6000
## 5 1011.9       6000
## 6 1013.2      10000
temp2 <- df2$t2m
hist(temp2)

temp_mensual = df2 %>% 
  group_by(month) %>% 
  summarise(tmax = max(t2m, na.rm = TRUE), 
            tmin = min(t2m, na.rm = TRUE),
            tavg = mean(t2m, na.rm = TRUE))

Ahora vamos a explorar los datos de la red de PEMBU

https://www.ruoa.unam.mx/pembu/index.php?page=historical_facts

datos = read.csv("./Datos_Estaciones/2018-CCA-L1/2018-01-CCA-L1.CSV", skip = 8, header = T)
datos <- datos[,c(1:2)]
colnames(datos) <- c("Fecha_hora", "Temp")
sapply(datos, class)
##  Fecha_hora        Temp 
## "character" "character"
datos$Temp <- as.numeric(datos$Temp)
## Warning: NAs introduced by coercion
sapply(datos, class)
##  Fecha_hora        Temp 
## "character"   "numeric"
head(datos)
##            Fecha_hora Temp
## 1 2018/01/01 00:00:00   NA
## 2 2018/01/01 00:30:00 12.7
## 3 2018/01/01 01:00:00 12.8
## 4 2018/01/01 01:30:00 12.3
## 5 2018/01/01 02:00:00 12.5
## 6 2018/01/01 02:30:00 12.6
# Los datos de las estaciones se toman cada media hora (48 observaciones por mes)

dim(datos)
## [1] 1488    2
t <- 1488/31

# agrupar por días
t2 <- split(seq_len(1488), rep(1:31, each=48))

# separar los días del mes de enero
datos2 <- lapply(1:31, function(x) datos[t2[[x]],])

# calcular el promedio diario
temp.mean.diaria <- lapply(1:31, function(x) mean(datos2[[x]][,2], na.rm=TRUE))
temp.mean.diaria <- unlist(temp.mean.diaria)
hist(temp.mean.diaria)

plot(temp.mean.diaria, xlab="días")

# Calcular el mínimo diario
temp.min.diaria <- lapply(1:31, function(x) min(datos2[[x]][,2], na.rm=TRUE))
temp.min.diaria <- unlist(temp.min.diaria)
hist(temp.min.diaria)

plot(temp.min.diaria, xlab="días")

# Calcular el máximo diario
temp.max.diaria <- lapply(1:31, function(x) max(datos2[[x]][,2], na.rm=TRUE))
temp.max.diaria <- unlist(temp.max.diaria)
hist(temp.max.diaria)

plot(temp.max.diaria, xlab="días")

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":
## Coordinate Reference System:
## Deprecated Proj.4 representation: NA
## Warning in wkt(x): CRS object has no comment

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)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
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":
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

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="./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) 
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.0, November 1st, 2021, [PJ_VERSION: 820]
## Path to PROJ shared files: /home/juvelas/.local/share/proj:/usr/share/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-7
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
ame1 <- readOGR(dsn="./Shapefile", layer="cont_america")
## OGR data source with driver: ESRI Shapefile 
## Source: "/mnt/tb/DD/Dropbox/Curso_PEMBU/Shapefile", layer: "cont_america"
## with 28 features
## It has 1 fields
ame2 <- readOGR("./Shapefile/cont_america.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/mnt/tb/DD/Dropbox/Curso_PEMBU/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)

!Eso 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) y https://cran.r-project.org/web/packages/climate/vignettes/getstarted.html (Czernecki, B.; Głogowski, A.; Nowosad, J. Climate: An R Package to Access Free In-Situ Meteorological and Hydrological Datasets for Environmental Assessment. Sustainability 2020, 12, 394. https://doi.org/10.3390/su12010394“)