Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WKT string as srs in cube_view #101

Open
jvandoninck opened this issue Jul 3, 2024 · 0 comments
Open

WKT string as srs in cube_view #101

jvandoninck opened this issue Jul 3, 2024 · 0 comments

Comments

@jvandoninck
Copy link

There seems to be a strange behaviour when defining the srs in gdalcubes::cube_view as a WKT. I created a data cube, but used terra:crs() to set the srs from an existing spatRaster object. In the reproducible example below a similar way to set srs using the example from the package.

library(gdalcubes)
library(terra)

#Data cube (from L8 cube gdalcubes package)
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"), quiet = TRUE)

L8.col = image_collection(file.path(tempdir(), "L8.db"))
v = cube_view(extent=list(left=388941.2, right=766552.4, 
                          bottom=4345299, top=4744931, t0="2018-01", t1="2018-06"),
              srs=terra::crs("EPSG:32618"), nx = 497, ny=526, dt="P1M",
              aggregation="first", resampling="bilinear")
L8.cube = raster_cube(L8.col, v) 
L8.cube = select_bands(L8.cube, c("B02", "B03", "B04", "B05", "B06"))

No problem so far. An error only occurs when using this data cube in a subsequent reduce_time() function, but only when using the FUN argument, and not the expr argument. In the example below two ways to create a median value composite. I have no clue what is causing this or whether it is machine dependent. When using the proj4 or EPSG notation both functions execute correctly.

#Median compositing function 1 - with expr
cube_composite.median <- function(cube){
  names <- gdalcubes::bands(cube)$name
  comp_args <- list(x=cube)
  comp_args <- c(comp_args, lapply(names, function(x){paste0("median","(",x,")")}))
  composite <- do.call(gdalcubes::reduce_time, comp_args)

  return(composite)
}

#Median compositing function 2 - with FUN
cube_composite.median2 <- function(cube){
  names <- gdalcubes::bands(cube)$name
  composite <- gdalcubes::reduce_time(cube, names=names, FUN=function(x){
    x_median <- apply(x, 1, median, na.rm=TRUE)
    return(x_median)
  })
  return(composite)
}

#Apply the two median composite functions 
L8.median <- cube_composite.median(L8.cube)
plot(L8.median, rgb=3:1)

L8.median2 <- cube_composite.median2(L8.cube)
plot(L8.median2, rgb=3:1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant