感谢
这个
帖子，我已经能够计算出一些 KUD，并将其屏蔽到形状文件的区域中，然后使用
terra
软件包中的
ggplot
和
SpatRaster
绘制这些图形
不过，我想将其映射到 ggmap 上。我当前的代码如下。
# 从 shapefile 创建 sf 对象，项目到 WGS84 UTM zone 31N/EPSG:32631
Abberton_utm <- st_read("WFD_Lake_Water_Bodies_Cycle_2.shp")%>%。
st_as_sf() %>%
st_transform(32631)
# 过滤数据以去除不准确的位置
AB_filt <- AB_DAT %>%
dplyr::filter(HPE <= 2)
# 绘制 KUD 数据的代码
AB_KUD <- AB_filt %>%
dplyr::arrange(Transmitter, DateTime) %>%
#dplyr::filter(week == "24") %>%
dplyr::select(Species, X_UTM, Y_UTM)
coordinates(AB_KUD) <- c("X_UTM","Y_UTM")
crs(AB_KUD) <- "EPSG:32631"；
# 创建 estUDm 对象
kud <- kernelUD(AB_KUD[,1], h = "href") # href = 引用带宽
# 为 rasterize() 创建 SpatRaster 模板
# 只需执行一次，无需在循环中复制
r <- rast(ext(Abberton_utm), ncol = 100, nrow = 100, crs = "EPSG:32631")
# 从 kud 中创建第 95 百分位数多边形 sf，用于初始 ggplot()
sf_poly <- getverticeshr(kud, percent = 95) %>%
st_as_sf() %>%
st_set_crs(32631)
# 将 sf_poly 转换为 SpatRaster，裁剪并屏蔽为 Abberton_utm
tmpr <- terra::rasterize(sf_poly, r) %>%
crop(., Abberton_utm, mask = TRUE)
# 为 tmpr 中的单元格分配百分位值
tmpr[] <- ifelse(is.na(tmpr[]), NA, 95)
# 创建初始图表
p <- ggplot() +
geom_spatraster(data = tmpr, show.legend = FALSE)
这只是计算 95% 的 KUD，并使用 ggplot 对其进行映射。不过，我希望用下面的代码在 ggmap 上绘制。
# 您获取地图的代码
Abberton_eastern <- c(0.875, 51.827)
AbRes <- get_map(location=Abberton_eastern, source="google", maptype='stamen_terrain', zoom=14)
随后，我使用了一个循环，该循环以 5% 的间隔计算 KUD，并将它们叠加在绘图中。但是，我似乎无法实现第一步。复制我的代码所需的所有数据都在 这里
明白了！你想在ggmap上叠加使用
adehabitatHR
生成的核密度估计（KDE）图，而不是ggplot。这是一个实现目标的细分方法，以及克服你面临障碍的代码示例：
1. 投影对齐
- ggmap使用WGS84（EPSG:4326）： ggmap从Google地图或Stamen设计等来源获取地图图块，这些地图图块使用WGS84坐标系。
-
将你的空间数据转换为WGS84：
使用
st_transform()将Abberton_utm、sf_poly和用于创建tmpr栅格的范围转换为EPSG:4326。
2. 从栅格到数据框进行ggmap兼容性
-
ggmap使用
geom_tile或geom_raster： 这些geoms需要一个具有经度、纬度和值的data.frame。 -
转换你的栅格：
使用
as.data.frame(rasterToPoints(tmpr))将tmpr栅格转换为适合ggmap的data.frame。
3. 在循环中绘制
- 针对每个百分比级别迭代： 在你的循环中，为每个KDE百分比级别创建单独的栅格，将其投影到WGS84，并转换为data.frame。
-
使用
geom_tile或geom_raster添加图层： 在ggmap()图上使用这些data.frame来叠加你的KDE。
代码示例
library(sf)
library(raster)
library(ggmap)
library(adehabitatHR)
library(terra)
# ...（你的数据加载和预处理代码）...
# 1. 投影对齐
Abberton_WGS84 <- st_transform(Abberton_utm, crs = 4326)
# 2. 创建一个函数来准备ggmap的栅格
raster_for_ggmap <- function(sf_polygon, raster_template, percentile) {
tmpr <- terra::rasterize(sf_polygon, raster_template) %>%
crop(., Abberton_WGS84, mask = TRUE)
tmpr[] <- ifelse(is.na(tmpr[]), NA, percentile)
# 转换为WGS84
tmpr_WGS84 <- project(tmpr, "+init=epsg:4326")
# 转换为data.frame
as.data.frame(rasterToPoints(tmpr_WGS84))
}
# 3. 使用ggmap绘制
# 获取地图
Abberton_eastern <- c(0.875, 51.827)
AbRes <- get_map(location = Abberton_eastern, source = "google",
maptype = 'stamen_terrain', zoom = 14)
# 创建基本ggmap
p <- ggmap(AbRes)
# 循环并添加KDE图层
percentages <- seq(5, 95, by = 5) # 设置要绘制的百分比级别
for (perc in percentages) {
# 为当前百分比创建sf多边形
sf_poly <- getverticeshr(kud, percent = perc) %>%
st_as_sf() %>%
st_set_crs(32631) %>%
st_transform(crs = 4326) # 转换为WGS84
# 为ggmap准备栅格
df_raster <- raster_for_ggmap(sf_poly, r, perc)
# 使用geom_tile添加KDE图层
p <- p +
geom_tile(data = df_raster, aes(x = x, y = y, fill = layer), alpha = 0.5)
}
# 自定义和显示最终图
p +
scale_fill_gradient(name = "KDE Percentile", low = "transparent", high = "blue") +
labs(title = "Kernel Density Estimation on ggmap")
这将创建一个ggmap，其中叠加了表示不同KDE百分比级别的多个半透明图层。调整颜色、透明度和整体美学以适合你的可视化偏好。
其他提示：
- 性能： 对于大量百分比级别或高分辨率栅格，此过程可能需要一段时间才能绘制。考虑在循环外预处理栅格以提高效率。
-
tmap包： 如果你想探索绘制空间数据的其他选项，tmap包提供了一种强大且灵活的方法来创建具有多个图层的静态和交互式地图。
祝你的ggmap KDE可视化好运！
标签：r,geospatial,ggmap,terra From： 78539515