Estadística Espacio Temporal

Datos espaciales y análisis exploratorio




logo

Maestría en Estadística Aplicada

Roberto Trespalacios

Temas

  • Datos espaciales y análisis exploratorio
    • Datos georeferenciados
      • Creación y tratamiento de datos tipo raster en R (librería raster)
      • Gráficos exploratorios
      • Aplicaciones

Creación y tratamiento de datos tipo raster en R (librería raster)

  • La librería raster soporta el análisis de datos espaciales en R.
  • raster define un conjunto de clases para representar datos espaciales.
library(raster)

#crear el raster
r <- raster()

hasValues(r) #tiene valores el raster?
ncol(r) <- nrow(r) <- 10 #num de col
dim(r) # dim del raster
xmax(r) # max en el eje horizontal(ext)
projection(r) <- "+proj=utm +zone=48 +datum=WGS84"

#agregar valores al raster 1:ncell(r)
values(r) <- 1:ncell(r)

#ejemplo
set.seed(0)
values(r) <- runif(ncell(r))

plot(r, main='Raster with 1000 cells')

Manejo de datos tipo raster

Vamos el ejemplo siguiente:

#filename <- "/home/rober/Downloads/geog5330-master/week8/Data/ext5.grd"
# filename <- "/home/rober/Downloads/ggmapTemp.png"
filename <- system.file("external/test.grd", package="raster")
filename

# Creamos el raster
r <- raster(filename)
filename(r)

hasValues(r)

inMemory(r)

Grafiquemos el raster

plot(r, main='RasterLayer del file')

Cambiemos su resolución

aggregate.r <- aggregate(r, fact=2)
plot(aggregate.r, main='RasterLayer del file')

Otras funciones: crear raster de un raster múltiple

filename <- system.file("external/rlogo.grd", package="raster")
filename
b <- brick(filename)
r <- raster(b, layer=2)
plot(r)

locs = read.csv(file="/home/rober/Downloads/geog5330-master/week8/Data/trees2.csv", header=T)
# Generate psuedo-absence tree locations as background data
coordinates(locs)=c('lngitude', 'latitude')
#proj4string(locs) = CRS('+init=epsg:4326')
x = circles(coordinates(locs), d=50000, lonlat=T)
#bg = spsample(x@polygons, 1000, type='random', iter=1000)
# load the climate conditions raster datasets
clim = brick('/home/rober/Downloads/geog5330-master/week8/Data/ext5.grd')
plot(clim, 1, cex=0.5, legend=T, mar=par("mar"), xaxt="n", yaxt="n", main="Annual mean temperature")
map("state", xlim=c(-119, -113), ylim=c(33.5, 38), fill=F, col="black", add=T)
# presence of trees
points(locs, col='red')
# psuedo-absense
#points(bg, col='blue')

require(maps)
require(sp)
require(maptools)
map.county <- map("county", plot = FALSE, fill = TRUE, res = 0)
IDs <- sapply(strsplit(q$names, ":"), function(x) x[1])
q <- map2SpatialPolygons(q, IDs=q$names, proj4string=CRS("+proj=longlat +datum=WGS84"))


# map("state", xlim=c(-119, -113), ylim=c(33.5, 38), fill=F, col="black", add=T)

Álgebra de rasters

A continuación vemos algunas operaciones algebraicas que se pueden realizar con los raster.

>, >=, <, ==, ! # operadores logicos
abs, round, ceiling, floor, trunc, sqrt
log, log10, exp # logaritmos y exponentes
cos, sin, atan, tan # funciones trigonometricas
max, min # maximo y minimo
range, prod, sum # algebra basica
any, all # consultas

#raster vacio
r <- raster(ncol=10, nrow=10)
values(r) <- 1:ncell(r)

#algunas operaciones
s <- r + 10
s <- sqrt(s)
s <- s * r + 5
r[] <- runif(ncell(r))
r <- round(r)
r <- r == 1

#reemplazo
s[r] <- -0.5
s[!r] <- 5
s[s == 5] <- 15

Observaciones importantes

  • Si usa múltiples objetos Raster (en funciones donde esto es relevante, como el rango), estos deben tener la misma resolución y origen.
  • El origen de un objeto Raster es el punto más cercano a (0, 0) que podría obtener si se movía desde las esquinas de un objeto Raster hacia ese punto en pasos de la resolución x e y.
  • Normalmente, estos objetos también tendrían la misma extensión, pero si no lo hacen, el objeto devuelto cubre la intersección espacial de los objetos utilizados.
  • Cuando utiliza múltiples objetos multicapa con diferentes números o capas, los objetos 'más cortos' son 'reciclados'. Por ejemplo, si multiplicas un objeto de 4 capas (a1, a2, a3, a4) con un objeto de 2 capas (b1, b2), el resultado es un objeto de cuatro capas (a1b1, a2b2, a3b1, a3b2).
r <- raster(ncol=5, nrow=5)
#r <- raster(ncol=2, nrow=2)
r[] <- 1
s <- stack(r, r+1)
q <- stack(r, r+2, r+4, r+6)
x <- r + s + q

Otras funciones(min, max, mean,...)

  • Las funciones de resumen (min, max, mean, prod, sum, Median, cv, range, any, all) siempre devuelven un objeto RasterLayer.
  • Quizás esto no sea obvio cuando se usan funciones como min, sum o mean.
a <- mean(r,s,10)
b <- sum(r,s)
st <- stack(r, s, a, b)
sst <- sum(st)
sst

La funciones cellStats, Aggregate y disaggregate, Crop y merge

  • Use cellStats si en lugar de un RasterLayer desea un solo número que resuma los valores de todas las celdas en cada capa.
cellStats(st, 'sum')
cellStats(sst, 'sum')
cellStats(st, 'mean')
  • Usamos Aggregate y disaggregate para cambiar la resolución del raster
r <- raster()
r[] <- 1:ncell(r)
ra <- aggregate(r, 20)
rd <- disaggregate(ra, 20)
  • Crop y merge para recortar y juntar raster.
r1 <- crop(r, extent(-50,0,0,30))
r2 <- crop(r, extent(-10,50,-20, 10))
m <- merge(r1, r2, filename='~/Desktop/test.grd', overwrite=TRUE)
plot(m)

plot of chunk unnamed-chunk-15

Reclasificar con cut o reclassify

cut o reclassify para reemplazar rangos de val. con val. únicos o subs para sust. val. indiv.

r <- raster(ncol=3, nrow=2)
r[] <- 1:ncell(r)
getValues(r)

Valores \( < 4 = \) NA

s <- calc(r, fun=function(x){ x[x < 4] <- NA; return(x)} )
as.matrix(s)

Realizar la operación r/(2*sqrt(y)) + 5

w <- overlay(r, s, fun=function(x, y){ x / (2 * sqrt(y)) + 5 } )
as.matrix(w)

Remover de r todos los valores que son NA en w

u <- mask(r, w)
as.matrix(u)

Identifica los valores de celda en u que son los mismos que en s.

v <- u==s
as.matrix(v)

Reemplazar NA's en w por los valores de r

cvr <- cover(w, r)
as.matrix(w)

Cambiar valores entre 0 y 2 por 0, etc.

x <- reclassify(w, c(0,2,1,  2,5,2, 4,10,3))
as.matrix(x)

Substituir 2 por 40 y por 50

y <- subs(x, data.frame(id=c(2,3), v=c(40,50)))
as.matrix(y)

Configuración espacial

La función clumpidentifica grupos de celdas que están conectadas.

  • edges, identifica bordes, es decir, transiciones entre valores de celda.
  • area, calcula el tamaño de cada celda de cuadrícula (para rásteres no proyectados),
    • Puede ser útil para, por ej. calcular el área cubierta por una determinada clasesobre una longitud/latitud en el raster.
  • zonal calcula estadísticas zonales, es decir, valores resumidos de un objeto Raster* para cada “zona” definida por un RasterLayer.
r <- raster(nrow=45, ncol=90)
r[] <- round(runif(ncell(r))*3)
a <- area(r)
zonal(a, r, 'sum')

Sumarización

r <- raster(ncol=36, nrow=18)
r[] <- runif(ncell(r))
cellStats(r, mean)

Estadística de zona

s <- r
s[] <- round(runif(ncell(r)) * 5)
zonal(r, s, 'mean')

Conteo de celdas

freq(s)
freq(s, value=3)

Celdas, filas, columnas y coordenadas 

library(raster)
r <- raster(ncol=36, nrow=18)
ncol(r) #n columnas

nrow(r) # n filas

ncell(r) # n celdas

rowFromCell(r, 100) # fila de la celda 100

colFromCell(r, 100) #col de la celda 100

cellFromRowCol(r,5,5) # celda de la fila y columna  

xyFromCell(r, 100) # coordenadas de la celda 100

cellFromXY(r, c(0,0)) #celda con la coordenada 0,0

colFromX(r, 0) #col de la coordenada x=0

rowFromY(r, 0) #fila de la coordenada y=0

Acceso a celdas

Se puede acceder a los valores de celda con varios métodos. Use getValues para obtener todos los valores o una sola fila; y getValuesBlock para leer un bloque (rectángulo) de valores de celda.

r <- raster(system.file("external/test.grd", package="raster"))
v <- getValues(r, 50)
v[35:39]
getValuesBlock(r, 50, 1, 35, 5)

También puede leer valores usando números de celda o coordenadas (x,y) usando el método de extracción.

cells <- cellFromRowCol(r, 50, 35:39)
cells
extract(r, cells)
xy <- xyFromCell(r, cells)
xy

Acceso a valores de celda

Raster tratado como vector

cells <- cellFromRowCol(r, 50, 35:39)
r[cells]
r[1:4]
r[2:3] <- 10
r[1:4]

Raster tratado como matriz

r[1]

r[2,2]

r[1,]

r[,2]

r[1:3,1:3]

# keep the matrix structure
r[1:3,1:3, drop=FALSE]