如何在 R terra 中裁剪、投影和合并来自反子午线两侧的栅格数据？

Question

我正在从全球 latlon 栅格 (EPSG: 4326) 中裁剪数据，然后希望将该数据投影到合适的局部坐标参考系中。我正在裁剪的区域跨越了反子午线，因此我想裁剪反子午线两侧的全局数据，投影这两个“一半”，然后在投影它们后将它们合并在一起。

但是，当我投影两半时，它们最终的分辨率略有不同。我可以重新采样以获得相同的分辨率，然后合并它们，但生成的合并数据缺少跨反子午线的一些值。

我想知道为什么两半的分辨率不同，以及如何在不产生数据间隙的情况下合并它们。

下面的 Reprex 说明了这个问题

library(terra)
#> terra 1.7.55
library(geodata)

#equal area projection for Fiji from projectionwizard.org
fiji_crs <- "+proj=laea +lon_0=-181.8896484 +lat_0=-17.73775 +datum=WGS84 +units=m +no_defs"

#get temperature data for Fiji
fiji_temp <- bio_oracle(path=tempdir(), "Temperature", "Mean", benthic=FALSE, time="Present")

#split data into left and right hand side of the antimeridian (just taking the first raster from the stack of rasters for this example)
fiji_temp_lhs <- crop(fiji_temp[[1]], ext(176, 180, -21, -12))
fiji_temp_rhs <- crop(fiji_temp[[1]], ext(-180, -177, -21, -12))

#project the two halves - they now have slightly different resolution
fiji_temp_list <- lapply(list(fiji_temp_lhs, fiji_temp_rhs), function(x) project(x, fiji_crs))

fiji_temp_list
#> [[1]]
#> class       : SpatRaster 
#> dimensions  : 109, 48, 1  (nrow, ncol, nlyr)
#> resolution  : 9164.313, 9164.313  (x, y)
#> extent      : -230080.8, 209806.2, -364282.8, 634627.3  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-17.73775 +lon_0=-181.8896484 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : Present.Surface.Temperature.Mean 
#> min value   :                         25.97523 
#> max value   :                         29.11796 
#> 
#> [[2]]
#> class       : SpatRaster 
#> dimensions  : 109, 37, 1  (nrow, ncol, nlyr)
#> resolution  : 9190.435, 9190.435  (x, y)
#> extent      : 196542.6, 536588.7, -368078.4, 633679.1  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-17.73775 +lon_0=-181.8896484 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : Present.Surface.Temperature.Mean 
#> min value   :                         26.02962 
#> max value   :                         29.06818

# resample the RHS raster to have the same resolution as the LHS
fiji_temp_list[[2]] <- resample(fiji_temp_list[[2]], rast(extent = ext(fiji_temp_list[[2]]), res = res(fiji_temp_list[[1]]), crs = crs(fiji_temp_list[[2]])))

#can now merge the resulting rasters but there are missing values where the two halves of the data have been merged
fiji_temp_list |>
  sprc() |>
  merge() |>
  plot()

^{创建于 2023-11-09，使用 reprex v2.0.2}

注意：我将裁剪成两半，而不是使用单个多边形/边界框，因为涉及的计算机处理较少：单个多边形裁剪会产生跨越整个赤道的栅格。我还想避免投影未裁剪的全局数据，因为这也需要大量的内存/处理时间。

Answer 1

不要裁剪和变换栅格，而是先变换然后裁剪，如下所示：

fiji_crs <- "+proj=laea +lon_0=-181.8896484 +lat_0=-17.73775 +datum=WGS84 +units=m +no_defs"

fiji_temp <- geodata::bio_oracle(path=tempdir(), "Temperature", "Mean", benthic=FALSE, time="Present")

e <- terra::ext(c(-230080.8, 536588.7, -364282.8, 634627.3))

fiji_temp[[1]] |>
  terra::project(y = fiji_crs) |>
  terra::crop(e) |>
  terra::plot()
                                         

crs 投影后要裁剪的坐标是从您的

fiji_temp_list

栅格（xmin1 和 xmin2 的最小值等）给出的，但是您几乎可以自动计算两个扩展的它们，例如：

t <- terra::rast()
terra::ext(t) <- c(176, 180, -21, -12)
terra::project(t, y = fiji_crs)
#> class       : SpatRaster 
#> dimensions  : 370, 162, 1  (nrow, ncol, nlyr)
#> resolution  : 2690.974, 2690.974  (x, y)
#> extent      : -230080.8, 205857, -361033.1, 634627.3  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-17.73775 +lon_0=-181.8896484 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

^{创建于 2024-01-02，使用 reprex v2.0.2}