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

wondering why terra::roll() will not match window_time() #104

Open
zackarno opened this issue Oct 3, 2024 · 0 comments
Open

wondering why terra::roll() will not match window_time() #104

zackarno opened this issue Oct 3, 2024 · 0 comments

Comments

@zackarno
Copy link

zackarno commented Oct 3, 2024

I'm probably misunderstanding something. But I cannot make terra::rolll() match the results of window_time() which i thought would be possible.

Below is a reprex where I try 2 methods for window_time() and then play with the parameters in terra::roll() in several different ways. Nonetheless I can't get them to match up.

library(terra)
#> terra 1.7.78
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Warning: package 'sf' was built under R version 4.4.1
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(gdalcubes)
#> 
#> Attaching package: 'gdalcubes'
#> The following objects are masked from 'package:terra':
#> 
#>     animate, crop, size

# create image collection from example Landsat data only 
# if not already done in other examples
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"), 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-07"),
              srs="EPSG:32618", nx = 400, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.nir = select_bands(L8.cube, c("B05"))

# first w/ mean expr
L8.nir.rollmean_expr = window_time(L8.nir, window = c(1,1), expr = "mean(B05)")  

# then with kernel
L8.nir.rollmean_kernel = window_time(L8.nir, window = c(1,1),kernel = rep(1/3,3))

# convert to terra raster so i can see the values better
window_result_expr <- terra::rast(stars::st_as_stars(L8.nir.rollmean_expr))
#> Warning in CPL_read_gdal(as.character(x), as.character(options),
#> as.character(driver), : GDAL Message 1: The dataset has several variables that
#> could be identified as vector fields, but not all share the same primary
#> dimension. Consequently they will be ignored.
window_result_kernel <- terra::rast(stars::st_as_stars(L8.nir.rollmean_kernel))
#> Warning in CPL_read_gdal(as.character(x), as.character(options),
#> as.character(driver), : GDAL Message 1: The dataset has several variables that
#> could be identified as vector fields, but not all share the same primary
#> dimension. Consequently they will be ignored.


# take original cube and convert to terra
r_from_cube <- terra::rast(stars::st_as_stars(L8.cube))
#> Warning in CPL_read_gdal(as.character(x), as.character(options),
#> as.character(driver), : GDAL Message 1: The dataset has several variables that
#> could be identified as vector fields, but not all share the same primary
#> dimension. Consequently they will be ignored.
r_from_cube
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :        7889,         NaN,        8616,        9259,        9782,        9375, ... 
#> max values  :       21454,         NaN,       31587,       50520,       40355,       47991, ... 
#> time (raw)  : 1 to 7

terra_result_na_T <- terra::roll(
  x=r_from_cube,
  n=3, 
  fun = "mean",
  na.rm=T
)
terra_result_na_T_circular <- terra::roll(
  x=r_from_cube,
  n=3, 
  fun = "mean",
  na.rm=T, 
  circular =T
)
terra_result_na_F <- terra::roll(
  x=r_from_cube,
  n=3, 
  fun = "mean",na.rm=F
)
terra_result_na_F_circular <- terra::roll(
  x=r_from_cube,
  n=3, 
  fun = "mean",na.rm=F, circular=T
)

window_result_expr
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :        5278,        5278,        5383,        5383,        5564,        5510, ... 
#> max values  :       27253,       33477,       57497,       45797,       45797,       48476, ... 
#> time (raw)  : 1 to 7
window_result_kernel
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :         NaN,         NaN,         NaN,        6395,    6580.333,    6528.333, ... 
#> max values  :         NaN,         NaN,         NaN,       38198,   38943.333,   35322.000, ... 
#> time (raw)  : 1 to 7
terra_result_na_T
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :        7889,        7889,        8829,      8829.0,      9259.0,        9394, ... 
#> max values  :       21454,       31587,       50520,     36977.5,     36977.5,       39473, ... 
#> time (raw)  : 1 to 7
terra_result_na_T_circular
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :        7889,        7889,        8829,      8829.0,      9259.0,        9394, ... 
#> max values  :       20388,       31587,       50520,     36977.5,     36977.5,       39473, ... 
#> time (raw)  : 1 to 7
terra_result_na_F
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :         NaN,         NaN,         NaN,    9978.333,    9861.333,    10357.33, ... 
#> max values  :         NaN,         NaN,         NaN,   32611.000,   32742.333,    31030.00, ... 
#> time (raw)  : 1 to 7
terra_result_na_F_circular
#> class       : SpatRaster 
#> dimensions  : 423, 400, 7  (nrow, ncol, nlyr)
#> resolution  : 944.028, 944.7565  (x, y)
#> extent      : 388941.2, 766552.4, 4345299, 4744931  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
#> source(s)   : memory
#> names       : time2~01-01, time2~02-01, time2~03-01, time2~04-01, time2~05-01, time2~06-01, ... 
#> min values  :         NaN,         NaN,         NaN,    9978.333,    9861.333,    10357.33, ... 
#> max values  :         NaN,         NaN,         NaN,   32611.000,   32742.333,    31030.00, ... 
#> time (raw)  : 1 to 7

Created on 2024-10-03 with reprex v2.1.0

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