我正在尝试在ggplot中将栅格层叠加到地图上。栅格层包含来自卫星标记的每个时间点的可能性表面。我还想在栅格层上设置累积概率(95%、75%、50%)。
我已经找出了如何在ggplot地图上显示栅格层,但是坐标不相互对齐。我尝试让它们都具有相同的投影,但似乎没有奏效... 我希望它们都适配我的模型边界(xmin = 149,xmax = 154,ymin = -14,ymax = -8.75)
附上我的r代码和结果图:
我已经找出了如何在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))
请告诉我是否有其他包/代码行可用于完成此操作!感谢任何帮助
rasterVis::gplot
直接在 ggplot 中添加栅格图。 - Sébastien Rochette