在R中的ggplot2中的地图上覆盖栅格图层? [英] Overlay raster layer on map in ggplot2 in R?
问题描述
我正在尝试将栅格图层叠加到ggplot中的地图上.栅格图层包含来自卫星标签的每个时间点的似然面.我还想在栅格图层上设置累积概率(95％，75％，50％).
我已经弄清楚了如何在ggplot地图上显示栅格图层，但是坐标彼此不对齐.我尝试使每个投影都具有相同的投影，但是似乎不起作用...我希望它们都适合我模型的边界(xmin = 149，xmax = 154，ymin = -14，ymax = -8.75
随附的是我的r代码和图形结果:
#load data
ncname <- "152724-13-GPE3"
ncfname <- paste(ncname, ".nc", sep = "")
ncin <- nc_open(ncfname)
StackedObject<-stack("152724-13-GPE3.nc", varname = "monthly_residency_distributions")
MergedObject<-overlay(StackedObject,fun=mean )
MergedObject[is.na(MergedObject)]<-0
Boundaries<-extent(c(149, 154, -14, -8.75))
ExtendedObject<-extend(MergedObject, Boundaries)
Raster.big<-raster(ncol=1200,nrow=900,ext=Boundaries)
Raster.HR<-resample(x=ExtendedObject, y=Raster.big, method="bilinear")
Raster.HR@data@values<- Raster.HR@data@values/sum(Raster.HR@data@values)
RasterVals<-sort(Raster.HR@data@values)
Raster.breaks <- c(RasterVals[max(which(cumsum(RasterVals)<= 0.05 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.25 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.50 ))], 1)
Raster.cols<-colorRampPalette(c("yellow","orange","red"))
RasterCols<- c(Raster.cols(3))
#Create Map
shape2 <- readOGR(dsn = "/Users/shannonmurphy/Desktop/PNG_adm/PNG_adm1.shp", layer = "PNG_adm1")
map<- crop(shape2, extent(149, 154, -14, -8.75))
projection(map)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
p <- ggplot() + geom_polygon(data = map, aes(x = long, y = lat, group = group), color = "black", size = 0.25) + coord_map()
projection(Raster.HR)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#plot raster and ggplot
par(mfrow=c(1,1))
plot(p)
par(mfrow=c(1,1), new = TRUE)
plot(Raster.HR, col=RasterCols, breaks=Raster.breaks, legend = NULL, bbox(map))
请让我知道我是否应该使用另外的程序包/代码行来执行此操作！感谢任何帮助
我知道了.您想在ggplot上绘制多个栅格图层，或者希望栅格对象位于背景多边形对象上方. rasterVis::gplot的问题在于它直接绘制栅格，不允许在上方或下方添加另一个栅格.您提醒我，我已经有了这种需要，并修改了功能gplot来检索数据作为小标题，这样您就可以根据需要依次使用dplyr和ggplot2进行尽可能多的播放.感谢您的提醒，我将其添加到当前的github库中供以后使用！
让我们使用一个可重现的示例来显示此功能:
创建数据集
- 创建世界地图作为
Raster用作背景栅格地图 - 创建数据栅格，这里是到点的距离(限制为最大距离)
代码:
library(raster)
# Get world map
library(maptools)
data(wrld_simpl)
# Transform World as raster
r <- raster(wrld_simpl, res = 1)
wrld_r <- rasterize(wrld_simpl, r)
# Lets create a raster of data
pt1 <- matrix(c(100,0), ncol = 2)
dist1 <- distanceFromPoints(r, pt1)
values(dist1)[values(dist1) > 5e6] <- NA
plot(dist1)
# Plot both
plot(wrld_r, col = "grey")
plot(dist1, add = TRUE)
提取栅格值(的一部分)并进行小波变换的功能
#' Transform raster as data.frame to be later used with ggplot
#' Modified from rasterVis::gplot
#'
#' @param x A Raster* object
#' @param maxpixels Maximum number of pixels to use
#'
#' @details rasterVis::gplot is nice to plot a raster in a ggplot but
#' if you want to plot different rasters on the same plot, you are stuck.
#' If you want to add other information or transform your raster as a
#' category raster, you can not do it. With `SDMSelect::gplot_data`, you retrieve your
#' raster as a tibble that can be modified as wanted using `dplyr` and
#' then plot in `ggplot` using `geom_tile`.
#' If Raster has levels, they will be joined to the final tibble.
#'
#' @export
gplot_data <- function(x, maxpixels = 50000) {
x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)
coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))
## Extract values
dat <- utils::stack(as.data.frame(raster::getValues(x)))
names(dat) <- c('value', 'variable')
dat <- dplyr::as.tbl(data.frame(coords, dat))
if (!is.null(levels(x))) {
dat <- dplyr::left_join(dat, levels(x)[[1]],
by = c("value" = "ID"))
}
dat
}
在ggplot中绘制多个栅格
您可以使用gplot_data将任何栅格转换为小标题.然后，您可以使用dplyr添加任何修改，并使用geom_tile在ggplot上进行绘制.有趣的一点是，只要fill选项具有可比性，您可以根据需要在任意时间使用geom_tile多次使用不同的栅格数据.否则，您可以使用以下技巧在背景栅格地图中删除NA值，并使用唯一的fill颜色.
# With gplot_data
library(ggplot2)
# Transform rasters as data frame
gplot_wrld_r <- gplot_data(wrld_r)
gplot_dist1 <- gplot_data(dist1)
# To define a unique fill colour for the world map,
# you need to remove NA values in gplot_wrld_r which
# can be done with dplyr::filter
ggplot() +
geom_tile(data = dplyr::filter(gplot_wrld_r, !is.na(value)),
aes(x = x, y = y), fill = "grey20") +
geom_tile(data = gplot_dist1,
aes(x = x, y = y, fill = value)) +
scale_fill_gradient("Distance",
low = 'yellow', high = 'blue',
na.value = NA) +
coord_quickmap()
在多边形上绘制栅格
当然，将背景图作为多边形对象，此技巧还可以让您在其上添加栅格:
wrld_simpl_sf <- sf::st_as_sf(wrld_simpl)
ggplot() +
geom_sf(data = wrld_simpl_sf, fill = "grey20",
colour = "white", size = 0.2) +
geom_tile(data = gplot_dist1,
aes(x = x, y = y, fill = value)) +
scale_fill_gradient("Distance",
low = 'yellow', high = 'blue',
na.value = NA)
gplot_data现在在此简单的R包中: https://github.com/statnmap/cartomisc
I am trying to overlay a raster layer onto a map in ggplot. The raster layer contains likelihood surfaces for each time point from a satellite tag. I also want to set cumulative probabilities(95%, 75%, 50%) on the raster layer.
I have figured out how to show the raster layer on the ggplot map, but the coordinates are not aligned with one another. I tried making each have the same projection but it does not seem to be working... I want them both to fit the boundaries of my model (xmin = 149, xmax = 154, ymin = -14, ymax = -8.75
Attached is my r code and the figure result:
#load data
ncname <- "152724-13-GPE3"
ncfname <- paste(ncname, ".nc", sep = "")
ncin <- nc_open(ncfname)
StackedObject<-stack("152724-13-GPE3.nc", varname = "monthly_residency_distributions")
MergedObject<-overlay(StackedObject,fun=mean )
MergedObject[is.na(MergedObject)]<-0
Boundaries<-extent(c(149, 154, -14, -8.75))
ExtendedObject<-extend(MergedObject, Boundaries)
Raster.big<-raster(ncol=1200,nrow=900,ext=Boundaries)
Raster.HR<-resample(x=ExtendedObject, y=Raster.big, method="bilinear")
Raster.HR@data@values<- Raster.HR@data@values/sum(Raster.HR@data@values)
RasterVals<-sort(Raster.HR@data@values)
Raster.breaks <- c(RasterVals[max(which(cumsum(RasterVals)<= 0.05 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.25 ))], RasterVals[max(which(cumsum(RasterVals)<= 0.50 ))], 1)
Raster.cols<-colorRampPalette(c("yellow","orange","red"))
RasterCols<- c(Raster.cols(3))
#Create Map
shape2 <- readOGR(dsn = "/Users/shannonmurphy/Desktop/PNG_adm/PNG_adm1.shp", layer = "PNG_adm1")
map<- crop(shape2, extent(149, 154, -14, -8.75))
projection(map)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
p <- ggplot() + geom_polygon(data = map, aes(x = long, y = lat, group = group), color = "black", size = 0.25) + coord_map()
projection(Raster.HR)<- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#plot raster and ggplot
par(mfrow=c(1,1))
plot(p)
par(mfrow=c(1,1), new = TRUE)
plot(Raster.HR, col=RasterCols, breaks=Raster.breaks, legend = NULL, bbox(map))
Please let me know if there is another package/line of code I should be using to do this! Appreciate any help
Ok I understand. You want to plot multiple raster layers on the ggplot or you want that the raster object is over your background polygon object. The problem with rasterVis::gplot is that it directly plot the raster and does not allow to add another one under or over. You remind me that I already had this need and modified function gplot to retrieve the data as a tibble so that you can then play with it as much as you want with dplyr and then ggplot2. Thanks for the reminder, I added it in my current github library for later use!
Let's use a reproducible example to show this function:
Create datasets
- Create a map of the world as a
Rasterto be use as background Raster map - Create a raster of data, here a distance from a point (limited to a maximum distance)
The code:
library(raster)
# Get world map
library(maptools)
data(wrld_simpl)
# Transform World as raster
r <- raster(wrld_simpl, res = 1)
wrld_r <- rasterize(wrld_simpl, r)
# Lets create a raster of data
pt1 <- matrix(c(100,0), ncol = 2)
dist1 <- distanceFromPoints(r, pt1)
values(dist1)[values(dist1) > 5e6] <- NA
plot(dist1)
# Plot both
plot(wrld_r, col = "grey")
plot(dist1, add = TRUE)
Function to extract (part of) raster values and transform as a tibble
#' Transform raster as data.frame to be later used with ggplot
#' Modified from rasterVis::gplot
#'
#' @param x A Raster* object
#' @param maxpixels Maximum number of pixels to use
#'
#' @details rasterVis::gplot is nice to plot a raster in a ggplot but
#' if you want to plot different rasters on the same plot, you are stuck.
#' If you want to add other information or transform your raster as a
#' category raster, you can not do it. With `SDMSelect::gplot_data`, you retrieve your
#' raster as a tibble that can be modified as wanted using `dplyr` and
#' then plot in `ggplot` using `geom_tile`.
#' If Raster has levels, they will be joined to the final tibble.
#'
#' @export
gplot_data <- function(x, maxpixels = 50000) {
x <- raster::sampleRegular(x, maxpixels, asRaster = TRUE)
coords <- raster::xyFromCell(x, seq_len(raster::ncell(x)))
## Extract values
dat <- utils::stack(as.data.frame(raster::getValues(x)))
names(dat) <- c('value', 'variable')
dat <- dplyr::as.tbl(data.frame(coords, dat))
if (!is.null(levels(x))) {
dat <- dplyr::left_join(dat, levels(x)[[1]],
by = c("value" = "ID"))
}
dat
}
Plot multiple rasters in ggplot
You can use gplot_data to transform any raster as a tibble. You are then able to add any modification using dplyr and plot on ggplot with geom_tile. The interesting point is that you can use geom_tile as many time as you want with different raster data, provided that fill option is comparable. Otherwise, you can use the trick below to remove NA values in the background raster map and use a unique fill colour.
# With gplot_data
library(ggplot2)
# Transform rasters as data frame
gplot_wrld_r <- gplot_data(wrld_r)
gplot_dist1 <- gplot_data(dist1)
# To define a unique fill colour for the world map,
# you need to remove NA values in gplot_wrld_r which
# can be done with dplyr::filter
ggplot() +
geom_tile(data = dplyr::filter(gplot_wrld_r, !is.na(value)),
aes(x = x, y = y), fill = "grey20") +
geom_tile(data = gplot_dist1,
aes(x = x, y = y, fill = value)) +
scale_fill_gradient("Distance",
low = 'yellow', high = 'blue',
na.value = NA) +
coord_quickmap()
Plot raster over polygons
Of course, with a background map as a polygon object, this trick also let you add your raster over it:
wrld_simpl_sf <- sf::st_as_sf(wrld_simpl)
ggplot() +
geom_sf(data = wrld_simpl_sf, fill = "grey20",
colour = "white", size = 0.2) +
geom_tile(data = gplot_dist1,
aes(x = x, y = y, fill = value)) +
scale_fill_gradient("Distance",
low = 'yellow', high = 'blue',
na.value = NA)
EDIT: gplot_data is now in this simple R package: https://github.com/statnmap/cartomisc
这篇关于在R中的ggplot2中的地图上覆盖栅格图层?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
