1

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

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

2D 密度图 V1 - 完整
的 2D 密度图 V1 - 详细信息

这是代码（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, aes(x=X_WGS84, 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 中，如本文所述。所以我将我的 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, bottom = 46.52694, right = +16.90, top = 48.99), color = "color", 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) +#, extent = "panel", 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, aes(x=X_PM, y=X_PM, fill = ..density.., alpha=cut(..density..,breaks=c(-Inf,0.08,Inf))), contour = FALSE, bins=16, geom = 'raster', n=500) + 
  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=bl3857_df, aes(long, lat, group=group), color='black', fill='NA', inherit.aes = FALSE) +
  ggtitle("Schwarzkiefer 2016/2020") + xlab("X_3857") + ylab("X_3857")
  
hm_sk

现在，移位问题解决了，但密度图不再可见（仅绘制地图、点和边界），如您在此处看到的：

2D-Density Plot V2 - 完整
的 2D-Density Plot V2 - 详细信息

有什么建议，我如何制作一个正确对齐并包含密度图的图？提前非常感谢！

4

0 回答 0