使用 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 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

它们都具有相同的投影。我不知道出了什么问题？

编辑

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

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 月 月份的栅格给了我这个：

有谁知道我是否可以检查我的栅格是否反转为多边形边界？

Answer 1

这似乎与 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" 时会导致一些不一致，请在 rgee 1.0.6 发布之前始终使用 via = "drive"。

Answer 2

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

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