使用terra包转换栅格数据坐标系后出现意外结果求助

我有一套采用EPSG:4326（单位为度）的栅格数据，想要转换为单位为米的投影坐标系，选用了投影参数+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs，也可接受其他米制投影坐标系，但转换后得到意外结果——原本预期栅格单元值与EPSG:4326下完全一致，实际却有明显差异。

图示说明：转换前后的栅格可视化结果呈现明显不同，投影后的栅格出现值的模糊/变化

测试数据集：包含经纬度坐标与对应栅格值的CSV文件

使用代码：

data_df <- read.csv("C:/Users/Sophie/Downloads/Test.csv")
data_r <- tidyterra::as_spatraster(tibble::as_tibble(data_df, xy = TRUE), xycols = 1:2, crs = "EPSG:4326")
plot(data_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

test_r <- terra::project(data_r, "+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs")
plot(test_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

问题原因 #

terra::project()函数默认采用双线性插值作为重采样方法。地理坐标系转投影坐标系时，栅格的单元格位置、分辨率都会发生变化，插值计算会对原栅格值进行重新计算，最终导致转换后的单元值与原数据不一致。

解决方法 #

投影时指定重采样方法为最近邻插值（"near"），该方法会直接选取距离新单元格中心最近的原单元格值，能100%保留原栅格的单元值。修改后的代码如下：

test_r <- terra::project(
  data_r, 
  "+proj=lcc +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs",
  method = "near"
)
plot(test_r, col = rev(colorspace::divergingx_hcl(100, "Spectral")))

可选米制投影坐标系推荐 #

可根据数据所在区域选择合适的EPSG米制投影，例如：

• 北美区域：EPSG:26910（UTM Zone 10N）、EPSG:3857（Web墨卡托）
• 全球通用：EPSG:3857

使用EPSG代码替换投影参数即可，示例：

test_r <- terra::project(data_r, "EPSG:3857", method = "near")

内容的提问来源于stack exchange，提问作者Sophie Père