CODEWITHSUNDEEP
rrasterprojectiongoogle-earth-enginergee
89_Simple
89_Simple

Reputation: 3805

raster does not align with shapefile after processing with rgee

I defined a polygon:

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)

enter image description here

The polygon covers countries around south-east Asia

For each pixel in the polygon, I want to calculate monthly sum of a given band for a given year as follows: 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_imagecollection_to_local(
      ic = Imagecollection_pr,
      region = polygon,
      dsn = paste0('pr_')
)

Reading a single month's file

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

Since this raster covers southeast asian countries, I downloaded the shapefile

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

Plotting them on top of each other:

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

enter image description here

There is a misalignment and I am not sure if it is the right raster that has been processed for the given polygon. I also checked if their projection is same

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

Both of them have the same projection as well. I cannot figure out what went wrong?

EDIT

As suggested in comments, I defined a new polygon covering Australia as follows:

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)

enter image description here

and repeated the above code. Plotting the raster again for the month of March on polygon gives me this:

enter image description here

Does anyone know if I can check if my raster is reversed w.r.t to polygon boundaries?

Upvotes: 0

Views: 435

Answers (2)

csaybar
csaybar

Reputation: 179

This seems to be related to rgdal rather than to the raster package. Some raster downloaded from GEE have data flipped with respect to y. I solved this problem, as follow:

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')
}

Recent versions of the earthengine Python API causes some inconsistencies when via = "getInfo" is used, please always use via = "drive" until the release of rgee 1.0.6.

Upvotes: 1

Robert Hijmans
Robert Hijmans

Reputation: 47191

There does not seem to be a misalignment. To plot all these countries in one step, you could do

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

Upvotes: 1

Related Questions