R help es - Jan 2014 - Aumentar resolución raster

*Ejecutando este código creo un mapa de España en el que el color viene
dado por la altura.*

library(ggplot2)
library(raster)

# orografia de españa
dem <- getData('alt', country='ESP')[[1]]
slope <- terrain(dem, opt='slope')
aspect <- terrain(dem, opt='aspect')
hill <- hillShade(slope, aspect, 40, 270)
#fronteras de España
polyEsp<-getData('GADM',country='ESP',level=1)

#Cortando los valores

espHill<-mask(hill,polyEsp)
espHill<-crop(monHill,polyEsp)
espDem<-mask(dem,polyEsp)
espDem<-crop(espDem,polyEsp)

#transformando raster en points
hdf <- rasterToPoints(espHill); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")
ddf <- rasterToPoints(espDem); ddf <- data.frame(ddf)
colnames(ddf) <- c("X","Y","DEM")

#creando vector para color breaks
b.hs <- seq(min(hdf$Hill),max(hdf$Hill),length.out=100)
b.dem <- seq(min(ddf$DEM),max(ddf$DEM),length.out=100)

#
ggplot(NULL, aes(X, Y)) +
  geom_raster(data = ddf, aes(fill = DEM)) +
  geom_raster(data = hdf, aes(alpha = Hill)) +
  scale_fill_gradientn(name="Altitude",colours = rainbow(20))+
  guides(fill = guide_colorbar()) +
  scale_alpha(range = c(0, 0.5), guide = "none") +
  scale_x_continuous(name=expression(paste("Longitud
(",degree,")")),
expand=c(0,0)) +
  scale_y_continuous(name=expression(paste("Latitud
(",degree,")")),
expand=c(0,0)) +
  coord_equal()



*Tengo muchos datos de la altura. Si quiero representar un valor del que
tengo menos medidas, varío un poco el código. *

#creo los datos
datos <- replicate(10, 10*rnorm(10))
datosRaster <- raster(datos)
e <- extent(c(-10, 10, 35,43))   # limites latitud y longitud de españa
extent(datosRaster) <- e

# orografia de españa
dem <- getData('alt', country='ESP')[[1]]
slope <- terrain(dem, opt='slope')
aspect <- terrain(dem, opt='aspect')
hill <- hillShade(slope, aspect, 40, 270)

dem <- datosRaster  # ahora quiero representar esto

#fronteras de España
polyEsp<-getData('GADM',country='ESP',level=1)

#Cortando los valores

espHill<-mask(hill,polyEsp)
espHill<-crop(monHill,polyEsp)
espDem<-mask(dem,polyEsp)
espDem<-crop(espDem,polyEsp)

#transformando raster en points
hdf <- rasterToPoints(espHill); hdf <- data.frame(hdf)
colnames(hdf) <- c("X","Y","Hill")
ddf <- rasterToPoints(espDem); ddf <- data.frame(ddf)
colnames(ddf) <- c("X","Y","DEM")

#creando vector para color breaks
b.hs <- seq(min(hdf$Hill),max(hdf$Hill),length.out=100)
b.dem <- seq(min(ddf$DEM),max(ddf$DEM),length.out=100)

#
ggplot(NULL, aes(X, Y)) +
  geom_raster(data = ddf, aes(fill = DEM)) +
  geom_raster(data = hdf, aes(alpha = Hill)) +
  scale_fill_gradientn(name="Altitude",colours = rainbow(20))+
  guides(fill = guide_colorbar()) +
  scale_alpha(range = c(0, 0.5), guide = "none") +
  scale_x_continuous(name=expression(paste("Longitud
(",degree,")")),
expand=c(0,0)) +
  scale_y_continuous(name=expression(paste("Latitud
(",degree,")")),
expand=c(0,0)) +
  coord_equal()




*Me sale con unos píxeles muy grandes. *


*¿Cual sería la manera de que esto no ocurra?*

*si hago*

plot(datosRaster)
plot(datosRaster, interpolate=T)



*Consigo el efecto que quiero, pero no puedo unirlo al mapa de españa. *

*¿Hay algún método para interpolar un objeto Raster y darle más definición?*

*o alternativamente,*



*¿Puedo extrapolar la matriz datos para que tenga más filas y columnas?
Parecido a lo que hace bkde2D o hist2D cuando tengo puntos. *

*Gracias*

	[[alternative HTML version deleted]]

Hola,

ggplot necesita un data.frame para poder funcionar. Lo estás
construyendo con rasterToPoints. Es mejor que uses sampleRegular
porque no necesitas pintar todas las celdas del Raster. Por otra
parte, la interpolación gráfica que necesitas con datos de baja
resolución se consigue gracias a la función rasterGrob del paquete
grid, en el que ggplot2 está basado. La función geom_raster tienen la
opción interpolate=TRUE, y teóricamente la pasa a rasterGrob para que
haga el interpolado. Sin embargo, no funciona, pero existe un atajo:
https://github.com/hadley/ggplot2/issues/794.

Todo esto, aplicado a los datos de baja resolución y dicho en código,
va a continuación.

Saludos.

Oscar.

PD. Te puede ser útil la función gplot del paquete rasterVis. Es un
wrapper muy sencillo para que ggplot funcione "directamente" con un
objeto Raster*.
PD2. Si estás dispuesto a probar herramientas basadas en lattice en
lugar de ggplot2, te aconsejo que pruebes la función spplot del
paquete raster, y levelplot del paquete rasterVis. Precisamente he
publicado un ejemplo con un DEM de España y "hill shading" construido
con rasterVis::levelplot aquí:
http://oscarperpinan.github.io/spacetime-vis/images/thumbs/hillShading.png
y su código aquí
https://github.com/oscarperpinan/spacetime-vis/blob/master/raster.R#L60

-------------------------------------------------------------

datos <- replicate(10, 10*rnorm(10))
dem <- raster(datos)
e <- extent(c(-10, 10, 35,43))
extent(dem) <- e
## Hay que añadir la proyección para que terrain funcione
projection(dem) <- '+proj=longlat +datum=WGS84'

slope <- terrain(dem, opt='slope')
aspect <- terrain(dem, opt='aspect')
hill <- hillShade(slope, aspect, 40, 270)

## sampleRegular para extraer puntos.
## En este caso concreto no sirve de nada porque sólo hay 100 celdas,
## pero sí es útil para Raster normales con millones de celdas.
s <- stack(dem, hill)
names(s) <- c('dem', 'hill')
slr <- sampleRegular(s, 1e5, asRaster = TRUE)

## Construimos el data.frame completo incorporando coordenadas
coords <- xyFromCell(slr, seq_len(ncell(slr)))
dat <- as.data.frame(getValues(slr))
dat <- cbind(coords, dat)

## Dibujamos...
 ggplot(aes(x, y, fill=dem, alpha=hill), data=dat) +
    ## geom_raster(aes(fill = dem, alpha=hill), interpolate=TRUE) +
    layer(geom="raster", geom_params=list(interpolate=TRUE)) +
    scale_fill_gradientn(name="Altitude",colours = rainbow(20))+
    guides(fill = guide_colorbar()) +
    scale_alpha(range = c(0, 0.5), guide = "none") +
    scale_x_continuous(name=expression(paste("Longitud
(",degree,")")),
                       expand=c(0,0)) +
    scale_y_continuous(name=expression(paste("Latitud
(",degree,")")),
                       expand=c(0,0)) +
    coord_equal()