Cómo configurar el uso de ggplot2 para mapear un ráster
Me gustaría hacer una gráfica usando R studio similar a la siguiente (creada 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 tener en cuenta. En primer lugar, los shapefiles de la cuenca faltan en la versión R. eso está bien.
En segundo lugar, El fondo gris más oscuro en el gráfico R no es ningún valor de datos. En Arc, no se muestran, pero en R mostrar con gplot. No se muestran cuando uso 'plot'del paquete raster:
plot(test)
Mis preguntas son las siguientes:
- ¿Cómo puedo eliminar el relleno de NoData gris oscuro en el ejemplo 'gplot'?
- ¿Cómo puedo configurar la leyenda (barra de color) para que sea razonable (como en las leyendas de 'parcela' de ArcMap y raster?)
- ¿Cómo puedo controlar el mapa de colores?
Para tener en cuenta, he probado muchas versiones diferentes de
scale_fill_brewer
scale_fill_manual
scale_fill_gradient
Y así en y así sucesivamente, pero tengo 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 tengo una solución para conspirar uno de estos mapas, me gustaría parcela de varios mapas en una figura y crear un único colorbar para todo el panel (es decir, un colorbar para todos los mapas) y me gustaría ser capaz de controlar donde el colorbar se encuentra y el tamaño de la colorbar. Aquí hay un ejemplo de lo que puedo hacer con grid.organizar, pero no puedo averiguar cómo establecer una sola barra de color:
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
El la salida es la siguiente:
El shapefile y el archivo raster están disponibles en el siguiente enlace:
Https://drive.google.com/open?id=0B5PPm9lBBGbDTjBjeFNzMHZYWEU
Muchas gracias por adelantado.
Devtools:: session_info()
Información de la sesión ---------------------------------------------------------------------------------------------------------------------
valor de ajuste
versión R versión 3.1.1 (2014-07-10)
sistema x86_64, darwin10.8.0
ui RStudio (0.98.1103)
language (EN)
intercalar en_US.UTF-8
tz America / Los_Angeles
Paquetes -------------------------------------------------------------------------------------------------------------------------
paquete * versión fecha fuente
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)
digerir 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)
gtable 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)
MASA 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)
rversions 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)
3 answers
Así es como lo haría, con rasterVis::levelplot:
Cargar cosas:
library(rgdal)
library(rasterVis)
library(RColorBrewer)
Lee las 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 roturas para la rampa de color definida con el argumento at a continuación).
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 desee trazar el contorno de la leyenda (incluidas sus marcas), en cuyo caso agregue axis.line=list(col='black') a la lista de colorkey args. Esto es necesario para anular la supresión general de cajas 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 es a menudo una mejor rampa para usar, a pesar de ser being en mi opinión.un poco feo. Los principales beneficios sobre algo como ColorBrewer RdYlBu son que los colores siguen siendo distintos cuando están desaturados, y las diferencias de color reflejan mejor las diferencias en los valores. Creo que RdYlBu es perfectamente accesible para Deuteranopía/Protanopía / Tritanopía daltonismo, sin embargo.
Aquí está la versión de 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í está cómo trazar varios rásters según lo solicitado.
Suponiendo que todos los rásters tienen la misma extensión, resolución, proyecciones, etc., puede apilarlos en un RasterStack, y luego usar levelplot en la pila. Puede pasar width como un elemento de la lista pasada a colorkey para controlar la altura de la leyenda ("width" es un poco contra-intuitivo, pero por defecto las leyendas son verticales). Si desea suprimir las etiquetas de tira por encima de cada panel (como he hecho a continuación - por defecto se etiquetan con los nombres de capa de la pila [ver names(s)]), puede agregar strip.border y strip.background a la lista pasada 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))
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2015-10-24 00:58:32
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 incluiste lo que estabas usando para hacer {[4] } así que hice esto:
test_spdf <- as(test, "SpatialPixelsDataFrame")
test_df <- as.data.frame(test_spdf)
colnames(test_df) <- c("value", "x", "y")
Y, entonces es solo cuestión de enviar eso + el shapefile 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, tho. Viridis es solo uno de los mejores para venir alrededor en un tiempo muy largo.
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"))
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2015-10-21 11:09:08
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.valor = NA.
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2018-06-30 11:10:00