appelmar / gdalcubes

Creating and analyzing Earth observation data cubes in R
https://gdalcubes.github.io
Other
120 stars 28 forks source link

Unexpected result with temporal aggregation #37

Closed goergen95 closed 4 years ago

goergen95 commented 4 years ago

I am getting a very strange result concerning temporal aggregation when initiating a cube view for yearly aggregated rasters. Here is a reproducible example based on the packages data. In fact, all of the data is from the year 2018, however, for the sake of my point, I am assuming that the GTiffs written to disk actually represent yearly aggregation values (this is how I am using the package in my use case). When a new image collection is initiated and the 1st of January is specified as the respective timestamp of the raster files the temporal aggregation delivers empty rasters. However, when a different date within the year is specified the temporal aggregation seems to work (at least for the several dates I tried so far). Can you point me in the right direction, whether this is a bug or am I not getting the complete picture of how the temporal aggregation works?

library(gdalcubes)
library(magrittr)

if (!file.exists(file.path(tempdir(), "L8.db"))) {
  L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"),
                         ".TIF", recursive = TRUE, full.names = TRUE)
  create_image_collection(L8_files, "L8_L1TP", file.path(tempdir(), "L8.db")) 
}

L8.col = image_collection(file.path(tempdir(), "L8.db"))
v = cube_view(extent = L8.col, dx = 500, dy = 500, srs = "EPSG:32618", dt = "P6M")
raster_cube(L8.col, v) %>%
  select_bands(c("B05", "B04")) %>%
  apply_pixel("(B05-B04)/(B05+B04)", "NDVI") %>%
  write_tif(dir =".", prefix = "gdalcubes_")

files = list.files(".", pattern = "gdalcubes_")
y.col = create_image_collection(files, date_time = c("2017-01-01", "2018-01-01"))
v2 = cube_view(extent = y.col, dx = 500, dy = 500, srs = "EPSG:32618", dt = "P1Y")

raster_cube(y.col, v2) %>%
  plot(zlim = c(-1,1))

y2.col = create_image_collection(files, date_time = c("2017-12-31", "2018-12-31"))
raster_cube(y2.col, v2) %>%
  plot(zlim = c(-1,1))

Output of the first plot call (1st of January):

1st_january

Output of the second plot call (31st of December):

31st_december

appelmar commented 4 years ago

Thanks a lot for your help! This was indeed a bug when images in a collection have year and month only (wrong string datetime comparison in the C++ part).