【问题标题】：Overlay raster layer on map in ggplot2 in R?在R中的ggplot2中的地图上覆盖栅格图层？
【发布时间】：2017-11-04 22:25:33
【问题描述】：

我正在尝试将栅格图层叠加到 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))

如果我应该使用另一个包/代码行来执行此操作，请告诉我！感谢任何帮助

【问题讨论】：

嗨，我刚刚注意到你拒绝了我的编辑。唯一的变化是内联该图像。你现在的状态还好吗？
嗨@yacc 对不起，我看错了，我以为你说要删除它。请随时再次编辑，对此感到抱歉。
没关系。编辑这样的问题很常见，因为 1) 新用户没有嵌入图像的权限（很有趣，但我认为这是一种防止垃圾邮件的方法），以及 2) 改进编辑的帖子会有所回报。例如，我要争取白色奖牌。这意味着要进行 80 次编辑才能获得银牌。 :) 这枚奖牌很有趣，我不知道，但这很有趣。那好吧。
可以使用rasterVis::gplot直接在ggplot中添加栅格。
@SébastienRochette 我试过了，但很难添加我想要的其他栅格参数，例如 50、75、95% 的可能性和与之配套的配色方案。我也有点不清楚究竟使用 rasterVis 包中的什么函数来将此栅格图层添加到 ggplot。

标签： r ggplot2 raster

【解决方案1】：

好的，我明白了。您想在 ggplot 上绘制多个栅格图层，或者您希望栅格对象位于背景多边形对象之上。 rasterVis::gplot 的问题在于它直接绘制栅格并且不允许在下方或上方添加另一个栅格。你提醒我，我已经有这个需求并修改了函数 gplot 以将数据作为 tibble 检索，这样你就可以使用 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

【讨论】：