terra 栅格值在投影后丢失

Question

SpatRaster 的值在更改坐标参考系统后丢失。我看不出任何原因。

library(terra)

ext <-
  terra::ext(
    9757195,
    9853641,
    734695,
    799794 
  )

r <-
  terra::rast(ext,
              resolution = 2000,
              crs = "EPSG:6933")

我创建了一个 SpatVector 点几何图形，然后与我的栅格叠加并确定这些点落在栅格的哪些像元中。此操作在预计的 CRS 中完成。

coord_vec <- data.frame( x = c(9849641, 9761195), y = c(795794.8, 738695.7))

coord_vec <- terra::vect(coord_vec, 
                         crs =  "EPSG:6933", geom=c("x", "y"))
r2_ <-
  terra::rasterize(x = coord_vec, y = r)

我想返回大地坐标系，但值丢失了。

r2_proj <- terra::project(x = r2_,
               y = "epsg:4326")

投影前的r2_spatraster

> r2_
class       : SpatRaster 
dimensions  : 33, 48, 1  (nrow, ncol, nlyr)
resolution  : 2000, 2000  (x, y)
extent      : 9757195, 9853195, 734695, 800695  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / NSIDC EASE-Grid 2.0 Global (EPSG:6933) 
source      : memory 
name        : lyr.1 
min value   :     1 
max value   :     1

投影后，值丢失。

> r2_proj 
class       : SpatRaster 
dimensions  : 27, 52, 1  (nrow, ncol, nlyr)
resolution  : 0.01927436, 0.01927436  (x, y)
extent      : 101.1252, 102.1275, 5.768228, 6.288636  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : memory 
name        : lyr.1 
min value   :   NaN 
max value   :   NaN

此工作流已经针对许多点和范围的数据集进行了测试，因此这种意外输出似乎是由这些点和范围的值生成的。

当我将 gdal 设置为 FALSE 时，我会获得非空值，因此这似乎是 GDAL-warp 算法的结果。

terra::project(x = r2_,
               y = "epsg:4326", gdal = F)

> terra::project(x = r2_,
+                y = "epsg:4326", gdal = F)
class       : SpatRaster 
dimensions  : 27, 52, 1  (nrow, ncol, nlyr)
resolution  : 0.01927436, 0.01927436  (x, y)
extent      : 101.1252, 102.1275, 5.768228, 6.288636  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : memory 
name        : lyr.1 
min value   :   0.5 
max value   :   0.5

Answer 1

这很神秘。也许你可以 raise an issue.

不过你想达到的效果，还是这样比较好。即避免对栅格数据进行变换，栅格化为你想要的栅格。

library(terra)
ext <-  ext( 9757195, 9853641, 734695, 799794  )
r <- rast(ext, resolution = 2000, crs = "EPSG:6933")
 
coord_vec <- data.frame( x = c(9849641, 9761195), y = c(795794.8, 738695.7))
coord_vec <- vect(coord_vec, crs="EPSG:6933", geom=c("x", "y"))

v <- project(coord_vec, "epsg:4326")
# project the raster template by using rast(r)
r2 <- project(rast(r), "epsg:4326")
r2 <- rasterize(v, r2)

terra 栅格值在投影后丢失

terra raster values lost after projection

r

raster

terra