更改栅格的crs以匹配简单特征点对象的crs

Question

假设我有一个 MODIS LAI 栅格和一个具有两个不同 crs 的简单要素 (

sf

) 对象（点）。我需要转换栅格的 crs 以匹配

sf

数据的 crs（反之亦然也可以，但我更愿意摆脱 MODIS 正弦系统）。各种尝试并没有让我找到解决方案，即将栅格和

sf

放在同一个 CRS 中以进行进一步的可视化、处理等。请帮助我如何做到这一点？以下是我尝试过的一些选项：

library(terra)
library(sf)

lai <- terra::rast("lai.tif")
terra::crs(lai) #check crs of lai

plots <- readRDS("plots.rds")
sf::st_crs(plots) #check crs of sf

#Trial 1
plots_tr <- sf::as_Spatial(plots) #transform the sf and get the crs
lai_pr1 <- terra::project(lai, "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") #project the raster using the crs of the transformed sf
terra::plot(lai_pr1) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

#Trial 2
lai_pr2 <- lai
terra::crs(lai_pr2) <- "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" # replace lai crs
terra::plot(lai_pr2) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

#Trial 3
lai_pr3 <- terra::project(lai, "epsg:25831") #project the raster using the EPSG code of the sf data (i.e., ETRS89 / UTM zone 31N)
terra::plot(lai_pr3) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

链接到运行此示例的数据：

lai<- https://github.com/frandadamo/frandadamo/blob/main/lai.tif

绘图<- https://github.com/frandadamo/frandadamo/blob/main/plots.rds

也许问题是

sf

的crs是WKT，而栅格的crs不是？此外，使用

sf

而不使用

sf::as_Spatial()

不会改变结果。我还尝试过用

st_transform()

更改 sf 的 crs。我阅读了已经发布的许多问题，但无法得到解决方案。

Maybe spTransform()

？

我希望这不是太愚蠢。非常感谢。

Answer 1

固定 MODIS 文件的坐标参考系后，使用

sf::st_transform()

或

terra::project()

非常简单，具体取决于您要保留的 crs。为了进行分析，我将继续重新投影矢量数据并保持栅格不变，但如果您愿意，可以随意删除 MODIS 投影（只是不要忘记重新采样）。另外，绝对不需要使用

{sp}

。

library(terra)
#> terra 1.7.71
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

lai <- terra::rast("lai.tif")
plots <- readRDS("plots.rds")

# fix MODIS crs
modis_crs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
terra::crs(lai) <- modis_crs

# option 1: reproject vector data 
plots_reproj <- sf::st_transform(plots, modis_crs)

terra::plot(lai, main = "MODIS Sinusoidal")
plot(plots_reproj, add = TRUE, pch = 19, col="black", cex = 0.1)

# option 2: reproject raster data
lai_reproj <- terra::project(lai, "epsg:25831")

terra::plot(lai_reproj, main = "EPSG:25831")
plot(plots, add = TRUE, pch = 19, col="black", cex = 0.1)

^{创建于 2024-04-27，使用 reprex v2.1.0}