如何解决用rgee处理后栅格未与shapefile对齐
我定义了一个多边形:
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)
多边形覆盖了东南亚地区的国家
对于多边形中的每个像素,我要计算给定年份的给定波段的每月总和,如下所示: 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
它们两者的投影也相同。我不知道出了什么问题?
编辑
如评论中所建议,我定义了一个覆盖澳大利亚的新多边形,如下所示:
polygon <- ee$Geometry$Polygon(
list(
c(88.75,-45.26),c(162.58,8.67),c(88.75,-45.26)
)
)
Map$addLayer(polygon)
并重复以上代码。再次在多边形上的 3月月中绘制栅格,可以得到以下信息:
有人知道我是否可以检查栅格是否反转为多边形边界吗?
解决方法
似乎没有错位。要一步一步绘制所有这些国家/地区,您可以
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。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。