使用terra::extract()提取栅格值时地图边界/日界线出现伪影

使用terra::extract()将全球栅格数据提取到投影坐标系下的多边形网格时，在地图右侧边界（日界线附近）出现虚假伪影，设置na.rm=FALSE时伪影消失但会丢失有效数据。示例代码如下：

示例数据与操作 #

# 加载所需库
library(tidyverse)
library(terra)
library(sf)

# 获取全球高程数据
elevation <- geodata::elevation_global(res = 10, download = TRUE, path = tempdir()) 

# 创建200x200km的方形网格（EPSG:6933米制投影）
grid_200km <- terra::rast(
  extent = ext(-17367530.4451614, 17367530.4451614, -6045743.25594052, 7166680.28573498), 
  crs = "EPSG:6933",
  resolution = 200000
) %>% as.polygons() %>% sf::st_as_sf()

问题重现代码 #

# 提取栅格值到网格
testExtract <- terra::extract(elevation, 
                              grid_200km,                        
                              fun=mean,                             
                              na.rm=TRUE,                            
                              bind=TRUE                             
) 

testExtract_sf <- st_as_sf(testExtract) %>% drop_na(wc2.1_10m_elev)

# 绘图可见右侧伪影
ggplot(testExtract_sf) +
  ggtitle("na.rm = TRUE（存在伪影）") +
  geom_sf(aes(fill = wc2.1_10m_elev)) +
  scale_fill_viridis_c() +  
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), legend.position="none")

成因分析 #

• 核心原因是EPSG:6933投影的全球连续性：该投影将全球映射为一个连续的矩形，左右边界对应地理坐标的±180°经线（日界线），在投影坐标系中这两个边界是“相连”的。
• terra::extract()处理最右侧的网格时，会错误地把日界线另一侧（最左侧）的栅格像素纳入计算——这些像素在真实地理上属于完全不同的区域，但投影坐标系下被判定为与网格重叠。
• 当na.rm=TRUE时，跨日界线的海洋NA像素被忽略，只计算了另一侧有陆地值的像素，从而产生虚假的异常值伪影；na.rm=FALSE时，混合NA后结果直接为NA，伪影被隐藏，但同时丢失了有效网格的数据。

解决方法 #

方法1：投影并裁剪栅格到网格范围（推荐） #

先将原始地理坐标栅格转换为与网格一致的投影，同时裁剪到网格的精确范围，彻底避免跨日界线的像素被误提取：

# 将高程栅格转换为EPSG:6933，同时裁剪到网格的范围
elevation_proj <- terra::project(elevation, "EPSG:6933", ext = terra::ext(grid_200km))

# 执行提取操作
testExtract <- terra::extract(elevation_proj, 
                              grid_200km,                        
                              fun=mean,                             
                              na.rm=TRUE,                            
                              bind=TRUE                             
) 

testExtract_sf <- st_as_sf(testExtract) %>% drop_na(wc2.1_10m_elev)

# 绘图验证伪影消失
ggplot(testExtract_sf) +
  ggtitle("修正后结果（投影+裁剪）") +
  geom_sf(aes(fill = wc2.1_10m_elev)) +
  scale_fill_viridis_c() +  
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), legend.position="none")

方法2：过滤跨日界线的网格单元 #

如果不需要转换投影，可以提前识别并排除跨越±180°经线的网格：

# 将网格转换为地理坐标（EPSG:4326），计算每个多边形的经度范围
grid_geo <- st_transform(grid_200km, "EPSG:4326")
grid_geo$lon_span <- st_bbox(grid_geo)[["xmax"]] - st_bbox(grid_geo)[["xmin"]]

# 过滤掉经度跨度超过170°的网格（正常网格不会跨日界线，跨度远小于180°）
grid_filtered <- grid_geo %>% 
  filter(lon_span < 170) %>% 
  st_transform("EPSG:6933")

# 使用过滤后的网格提取数据
testExtract <- terra::extract(elevation, 
                              grid_filtered,                        
                              fun=mean,                             
                              na.rm=TRUE,                            
                              bind=TRUE                             
) 

testExtract_sf <- st_as_sf(testExtract) %>% drop_na(wc2.1_10m_elev)

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