如何在一个Crs中合并一个Spatraster，延伸到另一个CR，然后合并延伸

Question

我有多个Sentinel-2颗粒，它们遍布多个UTM区域。我提出了一个简单的情况，即2个区域（例如，EPSG：32610和EPSG：32609）。我正在寻找一种有效的方法来使马赛克跨越更大的地理区域。

提供的脚本有几个问题。首先，这似乎是一个过于复杂的过程，必须有一种更简单的方法来解决此问题。对于想要横跨多个UTM区域的卫星图像或合并卫星图像的人们，这种情况可能经常遇到这种情况。我在使用 resample（）和 project（）的问题上遇到的问题是，分辨率从10 x 10 m变为带小数的稍微调整的分辨率（例如，10.01354 mx 10.02389 m ）；这是有道理的。

因此，循环首先将Spatraster的CRS放在目录中的CRS中 crss_r 。接下来，我使用一系列的IF语句列出了（a）在东部的较小CRS和（b）西部的较大CRS。在东方和西方共享相同CRS的Spatraster合并在一起（East Merge，West Merge）。东部设置为多边形的范围，该多边形足够大，可以包含东方和西部的所有尖峰。这应该允许东方和西部之间的合并，并与多边形进行最后的作物。

east1_2021 <- rast(ncols = 10, nrows = 10, xmin = 687000, xmax = 688000, ymin = 6253400, ymax = 6254400, crs="EPSG:32609")
evals1 <- sample(c(-1:1), 1000, replace=TRUE)

east2_2021 <- rast(ncols = 10, nrows = 10, xmin = 687500, xmax = 688500, ymin = 6253400, ymax = 6254400, crs="EPSG:32609")
evals2 <- sample(c(-1:1), 1000, replace=TRUE)

east3_2021 <- rast(ncols = 10, nrows = 10, xmin = 689300, xmax = 690300, ymin = 6253400, ymax = 6254400, crs="EPSG:32609")
evals3 <- sample(c(-1:1), 1000, replace=TRUE)

west1_2021 <- rast(ncols = 10, nrows = 10, xmin = 758600, xmax = 759600, ymin = 6290020, ymax = 6291020, crs="EPSG:32610")
wvals1 <- sample(c(-1:1), 1000, replace=TRUE)

west2_2021<- rast(ncols = 10, nrows = 10, xmin = 758800, xmax = 759800, ymin = 6290220, ymax = 6291220, crs="EPSG:32610")
wvals2 <- sample(c(-1:1), 1000, replace=TRUE)

west3_2021<- rast(ncols = 10, nrows = 10, xmin = 755800, xmax = 756800, ymin = 6288220, ymax = 6289220, crs="EPSG:32610")
wvals3 <- sample(c(-1:1), 1000, replace=TRUE)

values(east1_2021) <- evals1
values(east2_2021) <- evals2
values(east3_2021) <- evals3
values(west1_2021) <- wvals1
values(west2_2021) <- wvals2
values(west3_2021) <- wvals3

samp_rasters <- list()
samp_rasters[[1]] <- east1_2021
samp_rasters[[2]] <- east2_2021
samp_rasters[[3]] <- west1_2021
samp_rasters[[4]] <- west2_2021
samp_rasters[[5]] <- east3_2021
samp_rasters[[6]] <- west3_2021


E2021.list <- list()
E2021e.list <- list()
E2021w.list <- list()

crss_r <- data.frame()

for(i in 1:length(samp_rasters)){
  E2021.list[[i]] <- samp_rasters[[i]]
crss_r <- rbind(crss_r, as.numeric(crs(E2021.list[[i]], describe=TRUE)$code))
}

unique(crss_r)
#X32609
#1  32609
#4  32610

## Counters
  e = 1 # East counter
  w = 1 # West counter

  y = 1 # First in East
  z = 1 # First in West

  for(i in 1:length(samp_rasters)){
    
    if(crs(E2021.list[[i]], describe=TRUE)$code == 32609 & y == 1){
      e = 1
      E2021e.list[[e]] <- E2021.list[[i]]
      y <- y+1
    }
    
    else if(crs(E2021.list[[i]], describe=TRUE)$code == 32610 & z == 1){
      w = 1
      E2021w.list[[w]] <- E2021.list[[i]]
      z <- z+1
    }
    
    else if(crs(E2021.list[[i]], describe=TRUE)$code == 32609 & y != 1){
      e <- e+1
      E2021e.list[[e]] <- E2021.list[[i]]
    }

    else if(crs(E2021.list[[i]], describe=TRUE)$code == 32610 & z != 1){
      w <- w+1
      E2021w.list[[w]] <- E2021.list[[i]]
    }
  
    else{next}
  }

  ##E2021e <- sprc(E2021e.list) ## create a SpatRasterCollection east?
  ##E2021w <- sprc(E20212.list) ## create a SpatRasterCollection west? 
        
  E2021_east <- do.call(mosaic, E2021e.list)
  W2021_west <- do.call(mosaic, E2021w.list)
      

示例polygon ext（polygon.shp），我将最接近的偶数汇总到最接近的数字：

#SpatExtent : 606452.305334048, 869193.508446992, 6145080.67822892, 6378308.32505462 (xmin, xmax, ymin, ymax)

  e <- ext(606500, 869200, 6145100, 6378400)


  E2020_eastext <- terra::extend(E2020_east, ext(e)) ##
  W2020_westext <- terra::extend(W2020_west, ext(e)) ##

  E2021_eastext <- terra::extend(E2021_east, ext(polygon.shp))
  W2021_westext <- terra::extend(W2021_west, ext(polygon.shp))
  
  Full_2021 <- merge(E2021_eastext, W2021_westext)
  
  Full_2021 <- crop(Full_2021, polygon.shp, mask=TRUE)

使用我的实际数据在最终合并或马赛克中存在问题 - 来自不同区域的图像在错误的位置（湖泊）在西方，东方都在镜像）。

此外，一旦合并或马赛克对齐，我很想知道是否有更简单的方法可以跨区域完成这一目标？我计划运行这个越来越多的地理区域，这使这个过程比我想要的要复杂得多。

Answer 1

这是您的示例数据。但是我要问的第一个问题是您是如何陷入困境的？如果您使用的是Sentinel数据，它们并非都具有不同的范围？

library(terra)
east1_2021 <- rast(ncols = 10, nrows = 10, xmin = 687000, xmax = 688000, ymin = 6253400, ymax = 6254400, crs="EPSG:32609", vals=sample(c(-1:1), 100, replace=TRUE))
east2_2021 <- rast(ncols = 10, nrows = 10, xmin = 687500, xmax = 688500, ymin = 6253400, ymax = 6254400, crs="EPSG:32609", vals=sample(c(-1:1), 100, replace=TRUE))
east3_2021 <- rast(ncols = 10, nrows = 10, xmin = 689300, xmax = 690300, ymin = 6253400, ymax = 6254400, crs="EPSG:32609", vals=sample(c(-1:1), 100, replace=TRUE))
west1_2021 <- rast(ncols = 10, nrows = 10, xmin = 758600, xmax = 759600, ymin = 6290020, ymax = 6291020, crs="EPSG:32610", vals=sample(c(-1:1), 100, replace=TRUE))
west2_2021 <- rast(ncols = 10, nrows = 10, xmin = 758800, xmax = 759800, ymin = 6290220, ymax = 6291220, crs="EPSG:32610", vals=sample(c(-1:1), 100, replace=TRUE))
west3_2021 <- rast(ncols = 10, nrows = 10, xmin = 755800, xmax = 756800, ymin = 6288220, ymax = 6289220, crs="EPSG:32610", vals=sample(c(-1:1), 100, replace=TRUE))
samp_rasters <- list(east1_2021, east2_2021, east3_2021, west1_2021, west2_2021, west3_2021)

您在UTM 9和10中具有输入数据。

 sapply(samp_rasters, \(i)crs(i,proj=T)) |> unique()
 #[1] "+proj=utm +zone=9 +datum=WGS84 +units=m +no_defs" 
 #[2] "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"

坐标参考您希望合并的数据所在 - 因为通常不应该是两者中的任何一个（但请参见下面的情况下，有关您选择两者之一的情况）。我会假设经度/纬度。

outcrs <- "+proj=longlat"

一般方法可能是首先制作模板输出栅格。要制作一个，我首先计算输入数据的程度。

fun <- function(r) {
    p <- as.polygons(r, extent=TRUE) 
    crds(project(p, outcrs))
}

xy <- do.call(rbind, lapply(samp_rasters, fun) ) |> apply(2, range)
e <- round(ext(as.vector(xy)) + 0.05 , 1)
e
#SpatExtent : -126, -118.7, 56.3, 56.7 (xmin, xmax, ymin, ymax)

使用以下范围，然后选择与输入数据（100m 〜.01度）匹配的空间分辨率

temp <- rast(e, res=.01, crs=outcrs) 
temp
#class       : SpatRaster 
#dimensions  : 40, 730, 1  (nrow, ncol, nlyr)
#resolution  : 0.01, 0.01  (x, y)
#extent      : -126, -118.7, 56.3, 56.7  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 

现在将所有射手投射到输出模板上

y <- lapply(samp_rasters, \(x) project(x, temp))

结果合并

z <- merge(rast(y[1:3]), rast(y[4:6]))
names(z) <- paste0("step", 1:3)
z
#class       : SpatRaster 
#dimensions  : 40, 730, 3  (nrow, ncol, nlyr)
#resolution  : 0.01, 0.01  (x, y)
#extent      : -126, -118.7, 56.3, 56.7  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       :       step1,       step2,       step3 
#min values  : -0.68459326, -0.08613852, -0.99006349 
#max values  :   0.0181742,   0.6328486,  -0.1931030 

，并将您也可以与 tapp 这样的 （这实际上与 Mosaic 相同，但是与 Merge 相同，如果没有重叠值），

zz <- tapp(rast(y), rep(1:3, 2), mean, na.rm=TRUE)
names(zz) <- paste0("step", 1:3)

但请注意，一旦您拥有模板，您就可以进行快捷方式使用强加

e <- impose(sprc(samp_rasters[1:3]), temp)
w <- impose(sprc(samp_rasters[4:6]), temp)
m <- merge(e, w)

（为了提高效率，您可以通过裁剪整个区域模板的相关部分来使用不同的模板），

您不应需要扩展 at 代码> MERGE （或者 Mosaic ）可以考虑到这一点。

如果您想查看结果，则首先需要汇总值，因为在此示例中，有两个小区域，值介于两者之间。

a <- aggregate(m, 100, mean, na.rm=T)
plot(a,1)

---

替代方案。您可能会有一种情况，要保留两个来源之一的CR。对于UTM来说，这通常不是一个好主意，因为您可能会造成很多失真。但是，如果一个区域涵盖的数据不远，那是合理的。在这种情况下，您可以做

outcrs <- crs(samp_rasters[[1]])
xy <- do.call(rbind, lapply(samp_rasters, fun) ) |> apply(2, range)
e <- round(ext(as.vector(xy)))
temp <- rast(e, res=100, crs=outcrs) 
e <- impose(sprc(samp_rasters[1:3]), temp)
w <- impose(sprc(samp_rasters[4:6]), temp)
m <- merge(e, w)
m
#class       : SpatRaster 
#dimensions  : 767, 4391, 3  (nrow, ncol, nlyr)
#resolution  : 100, 100  (x, y)
#extent      : 687000, 1126100, 6253400, 6330100  (xmin, xmax, ymin, ymax)
#coord. ref. : WGS 84 / UTM zone 9N (EPSG:32609) 
#source      : memory 
#names       : lyr.1, lyr.1, lyr.1 
#min values  :    -1,    -1,    -1 
#max values  :     1,     1,     1