tutorial - mapas de calor en r

Cómo configurar usar ggplot2 para mapear un raster

Question

Me gustaría hacer un gráfico con R studio similar al de abajo (creado en Arc Map)

He intentado el siguiente código:

# data processing
library(ggplot2)
# spatial
library(raster)
library(rasterVis)
library(rgdal)

#
test <- raster(paste(datafold,''oregon_masked_tmean_2013_12.tif'',sep="")) # read the temperature raster
OR<-readOGR(dsn=ORpath, layer="Oregon_10N") # read the Oregon state boundary shapefile

gplot(test) +  
  geom_tile(aes(fill=factor(value),alpha=0.8)) + 
  geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
               fill=NA,color="grey50", size=1)+
  coord_equal()

La salida de ese código se ve así:

Un par de cosas a anotar. Primero, faltan los shapefiles de la cuenca en la versión R. eso está bien.

En segundo lugar, el fondo gris más oscuro en el gráfico R es valores Sin datos. En Arc, no se muestran, pero en R se muestran con gplot. No se muestran cuando uso ''plot'' del paquete raster:

plot(test)

Mis preguntas son las siguientes:

1. ¿Cómo me deshago del relleno gris oscuro NoData en el ejemplo ''gplot''?
2. ¿Cómo configuro la leyenda (barra de colores) para que sea razonable (como en las leyendas de ''trazado'' de ArcMap y ráster?)
3. ¿Cómo controlo el mapa de colores?

A tener en cuenta, he probado muchas versiones diferentes de

scale_fill_brewer
scale_fill_manual
scale_fill_gradient

y así sucesivamente pero me sale errores, por ejemplo

br <- seq(minValue(test), maxValue(test), len=8)

gplot(test)+
geom_tile(aes(fill=factor(value),alpha=0.8)) +

scale_fill_gradient(breaks = br,labels=sprintf("%.02f", br)) +

geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
             fill=NA,color="grey50", size=1)+
coord_equal()

Regions defined for each Polygons
Error: Discrete value supplied to continuous scale

Finalmente, una vez que tenga una solución para trazar uno de estos mapas, me gustaría trazar varios mapas en una figura y crear una única barra de colores para todo el panel (es decir, una barra de colores para todos los mapas) y me gustaría poder controla dónde se encuentra la barra de colores y el tamaño de la barra de colores. Este es un ejemplo de lo que puedo hacer con grid.arrange, pero no puedo entender cómo configurar una sola barra de colores:

r1 <- test
r2 <- test
r3 <- test
r4 <- test

colr <- colorRampPalette(rev(brewer.pal(11, ''RdBu'')))

l1 <- levelplot(r1, 
          margin=FALSE,                       
          colorkey=list(
             space=''bottom'',                   
             labels=list(at=-5:5, font=4),
             axis.line=list(col=''black'')       
          ),    
          par.settings=list(
             axis.line=list(col=''transparent'') 
          ),
          scales=list(draw=FALSE),            
          col.regions=viridis,                   
          at=seq(-5, 5, len=101)) +           
   layer(sp.polygons(oregon, lwd=3))

l2 <- levelplot(r2, 
                margin=FALSE,                       
                colorkey=list(
                   space=''bottom'',                   
                   labels=list(at=-5:5, font=4),
                   axis.line=list(col=''black'')       
                ),    
                par.settings=list(
                   axis.line=list(col=''transparent'') 
                ),
                scales=list(draw=FALSE),            
                col.regions=viridis,                   
                at=seq(-5, 5, len=101)) +           
   layer(sp.polygons(oregon, lwd=3))

l3 <- levelplot(r3, 
                margin=FALSE,                       
                colorkey=list(
                   space=''bottom'',                   
                   labels=list(at=-5:5, font=4),
                   axis.line=list(col=''black'')       
                ),    
                par.settings=list(
                   axis.line=list(col=''transparent'') 
                ),
                scales=list(draw=FALSE),            
                col.regions=viridis,                   
                at=seq(-5, 5, len=101)) +           
   layer(sp.polygons(oregon, lwd=3))

l4 <- levelplot(r4, 
                margin=FALSE,                       
                colorkey=list(
                   space=''bottom'',                   
                   labels=list(at=-5:5, font=4),
                   axis.line=list(col=''black'')       
                ),    
                par.settings=list(
                   axis.line=list(col=''transparent'') 
                ),
                scales=list(draw=FALSE),            
                col.regions=viridis,                   
                at=seq(-5, 5, len=101)) +           
   layer(sp.polygons(oregon, lwd=3))

grid.arrange(l1, l2, l3, l4,nrow=2,ncol=2) #use package gridExtra   

La salida es esta:

El archivo shape y el archivo raster están disponibles en el siguiente enlace:

https://drive.google.com/open?id=0B5PPm9lBBGbDTjBjeFNzMHZYWEU

Muchas gracias de antemano.

devtools :: session_info () Información de sesión ------------------------------------------ -------------------------------------------------- ------------------------- fijando el valor
versión R versión 3.1.1 (2014-07-10) sistema x86_64, darwin10.8.0
ui RStudio (0.98.1103)
idioma (EN)
cotejo en_US.UTF-8
tz America / Los_Angeles

Paquetes ------------------------------------------------- -------------------------------------------------- ---------------------- paquete * fuente fecha fecha
bitops 1.0-6 2013-08-17 CRAN (R 3.1.0) colorspace 1.2-6 2015-03-11 CRAN (R 3.1.3) devtools 1.8.0 2015-05-09 CRAN (R 3.1.3) digerido 0.6 .4 2013-12-03 CRAN (R 3.1.0) ggplot2 * 1.0.1 2015-03-17 CRAN (R 3.1.3) ggthemes * 2.1.2 2015-03-02 CRAN (R 3.1.3) git2r 0.10 .1 2015-05-07 CRAN (R 3.1.3) gridExtra 0.9.1 2012-08-09 CRAN (R 3.1.0) mesa 0.1.2 2012-12-05 CRAN (R 3.1.0) hexbin * 1.26. 3 2013-12-10 CRAN (R 3.1.0) celosía * 0.20-29 2014-04-04 CRAN (R 3.1.1) latticeExtra * 0.6-26 2013-08-15 CRAN (R 3.1.0) magrittr 1.5 2014 -11-22 CRAN (R 3.1.2) MASS 7.3-33 2014-05-05 CRAN (R 3.1.1) memoise 0.2.1 2014-04-22 CRAN (R 3.1.0) munsell 0.4.2 2013-07 -11 CRAN (R 3.1.0) plyr 1.8.2 2015-04-21 CRAN (R 3.1.3) proto 0.3-10 2012-12-22 CRAN (R 3.1.0) raster * 2.2-31 2014-03- 07 CRAN (R 3.1.0) rasterVis * 0.28 2014-03-25 CRAN (R 3.1.0) RColorBrewer * 1.0-5 2011-06-17 CRAN (R 3.1.0) Rcpp 0.11.2 2014-06-08 CRAN (R 3.1.0) RCurl 1.95-4.6 2015-04-24 CRAN (R 3.1.3) reshape2 1.4.1 2014-12-06 CRAN (R 3.1.2) rgdal * 0.8-16 2014-02-07 CRAN (R 3.1.0) versiones 1.0.0 2015-04-22 CRAN (R 3.1.3) escalas * 0.2.4 2014-04-22 CRAN (R 3.1.0) sp * 1.0-15 2014-04-09 CRAN (R 3.1.0) stringi 0.4-1 2014-12-14 CRAN (R 3.1.2) stringr 1.0.0 2015-04-30 CRAN (R 3.1.3) viridis * 0.3.1 2015 -10-11 CRAN (R 3.2.0) XML 3.98-1.1 2013-06-20 CRAN (R 3.1.0) zoo 1.7-11 2014-02-27 CRAN (R 3.1.0)

Answer 1

Así es como lo haría, con rasterVis::levelplot :

Cargar cosas

library(rgdal)
library(rasterVis)
library(RColorBrewer)

Leer cosas:

oregon <- readOGR(''.'', ''Oregon_10N'')
r <- raster(''oregon_masked_tmean_2013_12.tif'')

Defina una paleta de rampa de color (o un vector de colores con una longitud 1 más corta que el número de saltos para la rampa de color definida con el siguiente argumento).

colr <- colorRampPalette(brewer.pal(11, ''RdYlBu''))

Trazar cosas:

levelplot(r, 
          margin=FALSE,                       # suppress marginal graphics
          colorkey=list(
            space=''bottom'',                   # plot legend at bottom
            labels=list(at=-5:5, font=4)      # legend ticks and labels 
          ),    
          par.settings=list(
            axis.line=list(col=''transparent'') # suppress axes and legend outline
          ),
          scales=list(draw=FALSE),            # suppress axis labels
          col.regions=colr,                   # colour ramp
          at=seq(-5, 5, len=101)) +           # colour ramp breaks
  layer(sp.polygons(oregon, lwd=3))           # add oregon SPDF with latticeExtra::layer

Es posible que realmente desee trazar el contorno de la leyenda (incluidas sus marcas), en cuyo caso agregue axis.line=list(col=''black'') a la lista de colorkey de colorkey . Esto es necesario para anular la supresión general de cuadros causada por par.settings=list(axis.line=list(col=''transparent'')) :

levelplot(r, 
          margin=FALSE,                       
          colorkey=list(
            space=''bottom'',                   
            labels=list(at=-5:5, font=4),
            axis.line=list(col=''black'')       
          ),    
          par.settings=list(
            axis.line=list(col=''transparent'') 
          ),
          scales=list(draw=FALSE),            
          col.regions=colr,                   
          at=seq(-5, 5, len=101)) +           
  layer(sp.polygons(oregon, lwd=3))                  

Estoy de acuerdo con @hrbrmstr en que viridis suele ser una rampa mejor de usar , a pesar de ser, en mi opinión, un poco fea. Los principales beneficios sobre algo como RdYlBu de RdYlBu son que los colores aún son distintos cuando se desaturan, y las diferencias de color reflejan mejor las diferencias en los valores. Sin RdYlBu creo que RdYlBu es perfectamente accesible para el color ciego de Deuteranopia / Protanopia / Tritanopia.

Aquí está la versión viridis:

library(viridis)
levelplot(r, 
          margin=FALSE,                       
          colorkey=list(
            space=''bottom'',                   
            labels=list(at=-5:5, font=4),
            axis.line=list(col=''black'')       
          ),    
          par.settings=list(
            axis.line=list(col=''transparent'') 
          ),
          scales=list(draw=FALSE),            
          col.regions=viridis,                   
          at=seq(-5, 5, len=101)) +           
  layer(sp.polygons(oregon, lwd=3))

EDITAR

En respuesta a la pregunta adicional de OP, aquí se explica cómo trazar varios rásteres según se solicite.

Suponiendo que todos los rásteres tienen la misma extensión, resolución, proyecciones, etc., puede levelplot en un RasterStack y luego usar levelplot en la pila. Puede pasar el width como un elemento de la lista que se pasa a colorkey para controlar la altura de la leyenda ("ancho" es un poco contraintuitivo, pero por defecto las leyendas son verticales). Si desea suprimir las etiquetas de tira encima de cada panel (como lo he hecho a continuación, por defecto están etiquetadas con los nombres de capa de la pila [ver names(s) ]), puede agregar strip.border y strip.background a la lista par.settings a par.settings .

s <- stack(r, r*0.8, r*0.6, r*0.4)
levelplot(s, 
          margin=FALSE,                       
          colorkey=list(
            space=''bottom'',                   
            labels=list(at=-5:5, font=4),
            axis.line=list(col=''black''),
            width=0.75
          ),    
          par.settings=list(
            strip.border=list(col=''transparent''),
            strip.background=list(col=''transparent''),
            axis.line=list(col=''transparent'')
          ),
          scales=list(draw=FALSE),            
          col.regions=viridis,                   
          at=seq(-5, 5, len=101),
          names.attr=rep('''', nlayers(s))) +           
  layer(sp.polygons(oregon, lwd=3))

Answer 2

Esta es la solución fácil usando ggplot:

scale_fill_gradientn(colours = terrain.colors(4),limits=c(0,1),  
                 space = "Lab",name=paste("Probability /n"),na.value = NA)

En scale_fill_gradientn (también debería funcionar para scale_file_gradient), establezca na.value = NA.

Answer 3

library(ggplot2)
library(raster)
library(rasterVis)
library(rgdal)
library(grid)
library(scales)
library(viridis)  # better colors for everyone
library(ggthemes) # theme_map()

datafold <- "/path/to/oregon_masked_tmean_2013_12.tif"
ORpath <- "/path/to/Oregon_10N.shp"

test <- raster(datafold) 
OR <- readOGR(dsn=ORpath, layer="Oregon_10N") 

No incluyó lo que estaba usando para hacer la test así que hice esto:

test_spdf <- as(test, "SpatialPixelsDataFrame")
test_df <- as.data.frame(test_spdf)
colnames(test_df) <- c("value", "x", "y")

Y, luego, solo es cuestión de enviar ese + archivo de forma a ggplot2:

ggplot() +  
  geom_tile(data=test_df, aes(x=x, y=y, fill=value), alpha=0.8) + 
  geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
               fill=NA, color="grey50", size=0.25) +
  scale_fill_viridis() +
  coord_equal() +
  theme_map() +
  theme(legend.position="bottom") +
  theme(legend.key.width=unit(2, "cm"))

Funcionará con cualquier escala de temperatura continua ahora, aunque. Viridis es solo uno de los mejores para venir en mucho tiempo.

Puedes usar lo siguiente si tienes que usar gplot :

gplot(test) +  
  geom_tile(aes(x=x, y=y, fill=value), alpha=0.8) + 
  geom_polygon(data=OR, aes(x=long, y=lat, group=group), 
               fill=NA, color="grey50", size=0.25) +
  scale_fill_viridis(na.value="white") +
  coord_equal() +
  theme_map() +
  theme(legend.position="bottom") +
  theme(legend.key.width=unit(2, "cm"))