同UTM投影下SpatialPolygons绘图错位的成因与解决方法问询

问题成因 #

• CRS字符串格式错误：你定义的aeqd.crs和UTMstring嵌套了多余的双引号（如"\"+proj=aeqd...\""），传递给CRS()函数时会导致投影参数解析错误。CRS()需要纯proj4格式字符串，额外的双引号会让R错误识别投影规则，最终引发投影转换后的图层位置偏移。
• 投影转换未统一基准：两次调用projectRaster()时未指定统一的输出栅格模板，可能导致两个栅格的范围、分辨率不匹配，进一步放大错位问题。

解决方法 #

步骤1：修正CRS字符串格式 #

去掉嵌套双引号，正确定义投影参数：

library(magrittr)
library(raster)
library(rgeos)

# 修正CRS字符串，移除多余嵌套双引号
aeqd.crs <- "+proj=aeqd +lat_0=25.6871 +lon_0=-79.29336 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
UTMstring <- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"

步骤2：统一投影转换的目标栅格 #

先完成第一个栅格的投影转换，以此为基准转换第二个栅格，确保两者范围、分辨率完全对齐：

# 处理r1并转换到UTM
r1 <- raster("r1")
proj4string(r1) <- CRS(aeqd.crs)
r1 <- raster::setMinMax(r1)
r1_utm <- projectRaster(r1, crs = UTMstring) 

# 处理r2，以r1_utm为模板转换，强制对齐
r2 <- raster("r2")
proj4string(r2) <- CRS(aeqd.crs)
r2 <- raster::setMinMax(r2)
r2_utm <- projectRaster(r2, crs = UTMstring, res = res(r1_utm), ext = extent(r1_utm))

步骤3：重分类与重叠计算（栅格版，更高效） #

直接用栅格工具计算重叠，避免转换多边形的额外误差：

# 定义重分类规则
recl <- matrix(data = c(0, 0.05, 0,
                        0.05, 1, 1), ncol = 3, byrow = TRUE)

# 标准化r1并分类
r1_utm[is.na(r1_utm)] <- 0
r1_norm <- r1_utm / max(r1_utm[], na.rm = TRUE)
r1_recl <- reclassify(r1_norm, recl)
r1_recl[r1_recl == 0] <- NA

# 标准化r2并分类
r2_utm[is.na(r2_utm)] <- 0
r2_norm <- r2_utm / max(r2_utm[], na.rm = TRUE)
r2_recl <- reclassify(r2_norm, recl)
r2_recl[r2_recl == 0] <- NA

# 计算重叠区域栅格
overlap_raster <- overlay(r1_recl, r2_recl, fun = function(x, y) {
  return(ifelse(!is.na(x) & !is.na(y), 1, NA))
})

# 计算重叠面积（单位：平方公里）
common_area <- cellStats(area(overlap_raster), sum) / 1000000
print(common_area)

# 绘图验证
plot(r1_recl, col = "green", legend = FALSE)
plot(r2_recl, col = "red", add = TRUE, legend = FALSE)
plot(overlap_raster, col = "orange", add = TRUE, legend = FALSE)

若需转换为多边形计算（兼容原需求） #

确保栅格对齐后再转换，避免错位：

p1 <- as(r1_recl, 'SpatialPolygons')
p2 <- as(r2_recl, 'SpatialPolygons')
# 强制统一CRS（双重保险）
proj4string(p2) <- proj4string(p1)

overlap <- gIntersection(p1, p2)
common_area <- gArea(overlap) / 1000000
print(common_area)

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