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

编程问答 2022-08-15

问题描述

我定义了一个多边形：

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

Map$addLayer(polygon)

enter image description here

多边形覆盖了东南亚地区的国家

对于多边形中的每个像素，我要计算给定年份的给定波段的每月总和，如下所示： 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)

enter image description here

没有对齐，我不确定它是否是正确的栅格 处理给定的多边形。我还检查了他们的投影是否相同

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,8.67),c(88.75,-45.26)
   )
  )

 Map$addLayer(polygon)

enter image description here

并重复以上代码。再次在多边形上的 3月月中绘制栅格，可以得到以下信息：

enter image description here

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

解决方法

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

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

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

library(rgee)
library(raster)
ee_Initialize()

polygon <- ee$Geometry$Polygon(
  list(
    c(91.17,-13.42),c(154.10,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')
}

当使用via =“ getInfo”时，最新版本的Earthengine Python API引起一些不一致，请始终使用via =“ drive”直到发布rgee 1.0.6。

google-earth-engineprojectionrrrasterrgee