使用R ggplot绘制坦桑尼亚Mwanihana森林地图遇阻求助

问题背景 #

需要绘制坦桑尼亚Mwanihana森林区域的放大视图，目前仅能成功绘制数据点图层，但添加地形底图及坦桑尼亚地理范围时出现报错。

数据集示例 #

date     forestsite        primaryinterface secondaryinterface                   
<chr>    <chr>             <chr>            <chr>
8.27.13  UMNP-HQ           ecotourism       forest use/PD
8.28.13  UMNP-campsite3    ecotourism       forest use
8.28.13  UMNP-campsite3    ecotourism       forest use
8.28.13  UMNP-hondohondoa  ecotourism       forest use/PD
8.28.13  UMNP-hondohondob  ecotourism       forest use/PD/crop-raiding
8.29.13  UMNP-njokamoni    ecotourism       forest use
8.29.13  UMNP-mangabey     ecotourism       forest use   

anthrochange                     species       geometry
<chr>                            <chr>         <S3: sfc_POINT>      
small trails, snares             BWC           <S3: sfc_POINT>
roads, buildings                 RC            <S3: sfc_POINT>
roads, campsite                  RC            <S3: sfc_POINT>
roads, campsite                  BWC           <S3: sfc_POINT>
none                             SY            <S3: sfc_POINT>
clearcut road, houses, crops     YB            <S3: sfc_POINT>
trail to njokamoni               RC            <S3: sfc_POINT>
small trails, habituated troop   SM            <S3: sfc_POINT>

尝试代码 #

library(elevatr)
library(rgdal)

se.elevations <- get_elev_raster(locations = world, z = 7, clip = "locations") 
se.elevation.df <- raster::as.data.frame(se.elevations, xy = TRUE) %>%
  rename(elevation =  3) %>%
  na.omit() %>%
  filter(elevation >= 0)
xlimit <- c(36.85, 37)
ylimit <- c(-7.76, -7.88)

ggplot()+
  geom_tile(data = se.elevation.df, aes(x = x, y = y, fill = elevation))+
  scale_fill_gradientn(colours= c("gray", "black"))+
  geom_sf(data = world, fill = NA)+
  xlim(xlimit)+
  ylim(ylimit)+
  geom_sf(data = data, aes(color = forestsite, shape = primaryinterface), size = 3) +
    coord_sf(xlim = c(36.85, 37), ylim = c(-7.76, -7.88), expand = FALSE)

解决方法 #

核心问题分析 #

1. get_elev_raster使用world作为locations会获取全球高程数据，范围过大导致内存溢出或加载失败，且clip="locations"不适用于大区域裁剪
2. ggplot中同时使用xlim/ylimit和coord_sf(xlim/ylimit)会造成坐标系冲突
3. 未正确加载world对象，且目标是坦桑尼亚局部区域，无需加载全球边界

修正步骤 #

1. 安装并加载rnaturalearth包获取坦桑尼亚边界，缩小高程数据获取范围
2. 使用目标区域的边界框作为get_elev_raster的输入，改用clip="bbox"裁剪到目标范围
3. 移除xlim/ylimit，仅保留coord_sf的范围控制，避免坐标系冲突
4. 确保数据的坐标系与底图一致（均为WGS84，EPSG:4326）

修正后代码 #

library(elevatr)
library(ggplot2)
library(dplyr)
library(rnaturalearth)
library(sf)

# 获取坦桑尼亚边界
tanzania <- ne_countries(country = "tanzania", returnclass = "sf")

# 定义目标区域边界框
xlimit <- c(36.85, 37)
ylimit <- c(-7.88, -7.76) # 修正纬度范围顺序：从小到大
bbox <- st_bbox(c(xmin = xlimit[1], xmax = xlimit[2], ymin = ylimit[1], ymax = ylimit[2]), 
                crs = st_crs(4326))

# 获取目标区域高程数据
se.elevations <- get_elev_raster(locations = bbox, z = 10, clip = "bbox") # z值调高提升分辨率

# 转换为数据框用于ggplot
se.elevation.df <- raster::as.data.frame(se.elevations, xy = TRUE) %>%
  rename(elevation = 3) %>%
  na.omit() %>%
  filter(elevation >= 0)

# 绘制地图
ggplot() +
  # 地形底图
  geom_tile(data = se.elevation.df, aes(x = x, y = y, fill = elevation)) +
  scale_fill_gradientn(colours = c("#d9d9d9", "#8c8c8c", "#4a4a4a", "#1a1a1a"), 
                       name = "Elevation (m)") +
  # 坦桑尼亚边界
  geom_sf(data = tanzania, fill = NA, color = "black", size = 0.5) +
  # 数据点图层
  geom_sf(data = data, aes(color = forestsite, shape = primaryinterface), size = 3) +
  # 控制地图范围
  coord_sf(xlim = xlimit, ylim = ylimit, expand = FALSE, crs = 4326) +
  labs(title = "Mwanihana Forest, Tanzania", x = "Longitude", y = "Latitude") +
  theme_minimal()

额外注意事项 #

• 如果get_elev_raster下载缓慢，可尝试降低z值（分辨率），或检查网络连接
• 确保你的data对象的坐标系为EPSG:4326，若不是可使用st_transform(data, crs = 4326)转换
• 调整scale_fill_gradientn的颜色参数可优化地形可视化效果

内容的提问来源于stack exchange，提问作者Marnee Roundtree