投影后丢失的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_散斑器是

> 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

这太神秘了。也许你可以提出问题。

但对于你想要达到的目标，最好还是这样做。也就是说，避免将光栅数据转换为你想要的光栅。

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)