用 rgee 处理后栅格不与 shapefile 对齐

Question

我定义了一个多边形：

library(rgee)
ee_Initialize()
  
polygon <- ee$Geometry$Polygon(
list(
    c(91.17, -13.42), 
    c(154.10, -13.42), 
    c(154.10, 21.27), 
    c(91.17, 21.27),
    c(91.17, -13.42)
))

Map$addLayer(polygon)

多边形覆盖东南亚国家

对于多边形中的每个像素，我想计算给定年份给定波段的月总和，如下所示： month_vec <- 1:12
pr_ls <- 列表()

for(m in seq_along(month_vec)){

   month_ref <- month_vec[m]

   pr_ls[[m]] <-  
        ee$ImageCollection("NASA/NEX-GDDP")$ 
        filterBounds(polygon)$ # filter it by polygon
        select('pr')$ # select rainfall
        filter(ee$Filter$calendarRange(2000, 2000, "year"))$ # filter the year
        filter(ee$Filter$calendarRange(month_ref, month_ref, "month"))$ # filter the month 
        filter(ee$Filter$eq("model","ACCESS1-0"))$ # filter the model
        sum() # sum the rainfall
}

Imagecollection_pr <- ee$ImageCollection(pr_ls) 

ee_imagecollection_to_local(
      ic = Imagecollection_pr,
      region = polygon,
      dsn = paste0('pr_')
)

读取一个月的文件

my_rast <- raster(list.files(pattern = '.tif', full.names = TRUE)[1])

由于此栅格涵盖东南亚国家，因此我下载了 shapefile

sea_shp <- getData('GADM', country = c('IDN','MYS','SGP','BRN','PHL'), level = 0)

将它们绘制在彼此之上：

plot(my_rast)
plot(sea_shp, add = T)

有一个错位，我不确定它是否是正确的光栅为给定的多边形处理。我还检查了他们的投影是否相同

crs(my_rast)
CRS arguments: +proj=longlat +datum=WGS84 +no_defs
crs(sea_shp)      
CRS arguments: +proj=longlat +datum=WGS84 +no_defs

两者的投影也是一样的。我不知道哪里出了问题？

编辑

按照评论中的建议，我定义了一个覆盖澳大利亚的新多边形如下：

polygon <- ee$Geometry$Polygon(
   list(
     c(88.75,-45.26), 
     c(162.58,-45.26), 
     c(162.58,8.67), 
     c(88.75,8.67),
     c(88.75,-45.26)
   )
  )

 Map$addLayer(polygon)

并重复上面的代码。在多边形上再次绘制 3 月 月份的栅格给了我这个：

有谁知道我是否可以检查我的栅格是否被反转 w.r.t 到多边形边界？

Answer 1

好像没有错位。要一步绘制所有这些国家/地区，您可以这样做

x <- lapply(c('IDN','MYS','SGP','BRN','PHL'), function(i) getData('GADM', country = i, level = 0))
sea_shp <- bind(x)

Answer 2

这似乎与 rgdal 有关，而不是与 raster 包有关。从 GEE 下载的一些栅格数据相对于 y 翻转。我解决了这个问题，如下：

library(rgee)
library(raster)
ee_Initialize()

polygon <- ee$Geometry$Polygon(
  list(
    c(91.17, -13.42), 
    c(154.10, -13.42), 
    c(154.10, 21.27), 
    c(91.17, 21.27),
    c(91.17, -13.42)
  ))


month_vec <- 1:12
pr_ls <- list()


for(m in seq_along(month_vec)){
  month_ref <- month_vec[m]
  pr_ls[[m]] <-  
    ee$ImageCollection("NASA/NEX-GDDP")$ 
    filterBounds(polygon)$ # filter it by polygon
    select('pr')$ # select rainfall
    filter(ee$Filter$calendarRange(2000, 2000, "year"))$ # filter the year
    filter(ee$Filter$calendarRange(month_ref, month_ref, "month"))$ # filter the month 
    filter(ee$Filter$eq("model","ACCESS1-0"))$ # filter the model
    sum() # sum the rainfall
}

Imagecollection_pr <- ee$ImageCollection(pr_ls) %>% ee_get(0)

exp1 <- ee_imagecollection_to_local(
  ic = Imagecollection_pr,
  region = polygon,
  dsn = "pp_via_drive",
  via = "drive" # please always use "drive" or "gcs" until rgee 1.0.6 release
)

# One option
gdalinfo <- try (rgdal::GDALinfo(exp1))
if (isTRUE(attr(gdalinfo, "ysign") == 1)) {
  exp1_r <- flip(raster(exp1), direction='y')
}

最新版本的 earthengine Python API 在使用 via = "getInfo" 时会导致一些不一致，请始终使用 via = "drive" 直到 rgee 1.0.6 发布。

用 rgee 处理后栅格不与 shapefile 对齐

raster does not align with shapefile after processing with rgee

r

raster

projection

google-earth-engine

rgee