使用 ggmap 在地图顶部绘制密度二维图

编程问答 2022-04-22

问题描述

我想根据空间点数据生成二维密度图。在背景中，我想显示一个开放的地图（例如雄蕊地形）。此外我想绘制奥地利的边界。两个数据集（数据点和边界）都是 epsg 4326 中的 shapefile。

我设法制作了这样的图（请参阅下面的屏幕截图和代码 V1），但问题是一侧背景中的地图与另一侧的绘制点和奥地利边界之间存在偏移侧面，如下所示。

2D-Density Plot V1 - full
2D-Density Plot V1 - detail

这是代码（V1）：

library(sf)
library(rgdal)
library(ggplot2)
library(ggmap)

# Read point data (epsg: 4326)
sk <- st_read("points.shp") 

# Read country border polygon (epsg: 4326)
blogr4326 <- readOGR(<path>,<layer name>)
bl4326_df <- fortify(blogr4326)

# Austria Box (extent)
# Longitude: 9.6 to 16.94504
# Latitude: 46.52694 to 48.81667

map <- get_map(c(left = +9.6,bottom = 46.52694,right = +16.90,top = 48.99),color = "color",crop = FALSE)

hm_sk <- ggmap(map,extent = "panel",maprange=FALSE,darken=0.0) + 
  geom_point(data = sk,aes(x=X_wgs84,y=Y_wgs84)) +
  stat_density2d(data = sk,y=Y_wgs84,fill = ..density..,alpha=cut(..density..,breaks=c(-Inf,0.08,Inf))),contour = FALSE,bins=16,geom = 'raster',n=500) + 
  ggtitle("Schwarzkiefer 2016/2020") + xlab("X_wgs84") + ylab("X_wgs84") +
  scale_fill_distiller(palette= "Spectral",direction=-1,limits = c(0.08,8.50)) + 
  scale_alpha_manual(values=c(0,0.7),guide="none") +
  geom_polygon(data=bl4326_df,aes(long,lat,group=group),color='black',fill='NA',inherit.aes = TRUE) +
  coord_fixed(1.5)

hm_sk

我发现这种转变是由于背景中的地图在投影 epsg:3857 中而我的 shapefile 在投影 epsg:4326 中引起的，如本 post 中所述。所以我将我的 shapefile 投影到 epsg 3857 并将提供的代码插入到我的代码中，如您所见（V2）：

library(sf)
library(rgdal)
library(ggplot2)
library(ggmap)

# Read point data (epsg: 3857)
sk <- st_read("points.shp")

# Read country border polygon (epsg: 3857)
blogr3857 <- readOGR(<path>,<layer name>)
bl3857_df <- fortify(blogr3857)

# Austria Box (extent)
# Longitude: 9.6 to 16.94504
# Latitude: 46.52694 to 48.81667

map <- get_map(c(left = +9.6,crop = FALSE)

#-------------------------------------------------------------------------------------------------
# Following code according to this link to avoid the shift between map and country border polygon:
# https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-raster

# Define a function to fix the bBox to be in epsg:3857
ggmap_bBox <- function(map) {
  if (!inherits(map,"ggmap")) stop("map must be a ggmap object")
  # Extract the bounding Box (in lat/lon) from the ggmap to a numeric vector,# and set the names to what sf::st_bBox expects:
  map_bBox <- setNames(unlist(attr(map,"bb")),c("ymin","xmin","ymax","xmax"))
  
  # Convert the bBox to an sf polygon,transform it to 3857,# and convert back to a bBox (convoluted,but it works)
  bBox_3857 <- st_bBox(st_transform(st_as_sfc(st_bBox(map_bBox,crs = 4326)),3857))
  
  # Overwrite the bBox of the ggmap object with the transformed coordinates 
  attr(map,"bb")$ll.lat <- bBox_3857["ymin"]
  attr(map,"bb")$ll.lon <- bBox_3857["xmin"]
  attr(map,"bb")$ur.lat <- bBox_3857["ymax"]
  attr(map,"bb")$ur.lon <- bBox_3857["xmax"]
  map
}

# Use the function:
map <- ggmap_bBox(map)
#-------------------------------------------------------------------------------------------------


hm_sk <- ggmap(map,extent = "device",maprange=FALSE) +#,darken=0.0) + 
  coord_sf(crs = st_crs(3857)) + # f867orce the ggplot2 map to be in 3857
  geom_point(data = sk,aes(x=X_PM,y=Y_PM)) +
  stat_density2d(data = sk,y=X_PM,n=500) + 
  scale_fill_distiller(palette= "Spectral",8.50)) +
  scale_alpha_manual(values=c(0,guide="none") +
  geom_polygon(data=bl3857_df,inherit.aes = FALSE) +
  ggtitle("Schwarzkiefer 2016/2020") + xlab("X_3857") + ylab("X_3857")
  
hm_sk

现在，偏移问题解决了，但密度图不再可见（仅绘制了地图、点和边界），如下所示：

2D-Density Plot V2 - full
2D-Density Plot V2 - detail

任何建议，我如何生成正确对齐并包含密度图的图？非常感谢提前！

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

ggmapggplot2rrstat-density2d