如何在R语言底图上叠加核密度估计（KDE）展示物种家域

核心问题分析 #

1. ggmap API限制：Google地图服务现在强制要求API密钥，导致get_map报错，需改用免费无API限制的底图方案。
2. 投影不匹配：你的KDE多边形（kde50）使用UTM投影（epsg:32748），而底图多采用WGS84经纬度投影，必须统一投影才能叠加。
3. 栅格转换错误：raster(kde50)无效，因为kde50是多边形对象，栅格化需基于kernelUD的原始结果。

方案一：静态底图（ggplot2 + 矢量地图） #

无需API，使用开源矢量底图绘制静态叠加图：

# 加载必要包
library(sp)
library(ggplot2)
library(mapdata)

# 读取并处理GPS数据（保留你的原始代码）
gps_data_121130 <- read.csv('gps_data_121130.csv')
gps_points <- SpatialPoints(gps_data_121130[, c("Longitude", "Latitude")], 
                            proj4string = CRS("+proj=longlat +datum=WGS84"))
cord.UTM <- spTransform(gps_points, CRS("+init=epsg:32748"))
kde <- kernelUD(cord.UTM, h = "href", grid = 100, kern = "bivnorm")
kde50 <- getverticeshr(kde, percent = 50)

# 将KDE多边形转换回WGS84投影（与底图匹配）
kde50_wgs84 <- spTransform(kde50, CRS("+proj=longlat +datum=WGS84"))
# 转换为ggplot支持的data.frame格式
kde_df <- fortify(kde50_wgs84)

# 获取西雅图所在的金县矢量底图
king_county <- subset(map_data("county", region = "washington"), subregion == "king")

# 绘制叠加图
ggplot() +
  geom_polygon(data = king_county, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  geom_point(data = gps_data_121130, aes(x = Longitude, y = Latitude), 
             color = "blue", size = 1) +
  geom_polygon(data = kde_df, aes(x = long, y = lat, group = group), 
               fill = "red", alpha = 0.3, color = "darkred", lwd = 1) +
  coord_sf(xlim = range(gps_data_121130$Longitude) + c(-0.01, 0.01),
           ylim = range(gps_data_121130$Latitude) + c(-0.01, 0.01)) +
  labs(title = "50% KDE物种家域", x = "经度", y = "纬度") +
  theme_bw()

方案二：交互式底图（leaflet） #

生成可交互的网页地图，无需API，支持缩放、点击查看信息：

library(sp)
library(leaflet)

# 基于你的原始代码生成kde50后，转换投影
kde50_wgs84 <- spTransform(kde50, CRS("+proj=longlat +datum=WGS84"))

# 创建交互式地图
leaflet() %>%
  addTiles() %>%  # 添加OpenStreetMap免费底图
  addPolygons(data = kde50_wgs84, fillColor = "red", fillOpacity = 0.3, 
              color = "darkred", weight = 2) %>%
  addMarkers(data = gps_data_121130, lng = ~Longitude, lat = ~Latitude, 
             popup = ~paste("部署ID:", DeployID, "<br>时间:", Date)) %>%
  setView(lng = mean(gps_data_121130$Longitude), 
          lat = mean(gps_data_121130$Latitude), zoom = 14)

栅格化KDE的正确方法 #

若需要生成KDE密度栅格并叠加，需基于kernelUD的结果转换：

library(raster)

# 从kernelUD生成栅格
kde_raster <- raster(kde)
# 转换为WGS84投影
kde_raster_wgs84 <- projectRaster(kde_raster, crs = CRS("+proj=longlat +datum=WGS84"))

# 可在ggplot中叠加栅格
ggplot() +
  geom_raster(data = as.data.frame(kde_raster_wgs84, xy = TRUE), 
              aes(x = x, y = y, fill = layer)) +
  geom_point(data = gps_data_121130, aes(x = Longitude, y = Latitude), color = "white") +
  scale_fill_viridis_c(option = "plasma") +
  coord_sf()

内容的提问来源于stack exchange，提问作者Alison Meeth