R Under development (unstable) (2024-11-06 r87293 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> suppressPackageStartupMessages(library(stars))
> geomatrix = system.file("tif/geomatrix.tif", package = "stars")
> x = read_stars(geomatrix)
> # can stars reproduce what gdal does, by default?
> x2 = st_warp(x, use_gdal = TRUE, no_data_value = -9999)
> y = st_warp(x, x2)
> if (interactive()) { plot(x2, breaks = "equal", axes=TRUE) }
> if (interactive()) { plot(y, breaks = "equal", axes=TRUE) }
> names(x2) = names(y)
> all.equal(x2, y) # yes?
[1] TRUE
> 
> # does gdal reproduce with stars template object?
> (x2 = setNames(st_warp(x, y, use_gdal = TRUE, no_data_value=-1), "file.tif"))
stars object with 2 dimensions and 1 attribute
attribute(s):
          Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
file.tif    74     115    123 126.4289     132  255  224
dimension(s):
  from to  offset delta                refsys point x/y
x    1 25 1840902  5.22 WGS 84 / UTM zone 11N FALSE [x]
y    1 25 1144003 -5.22 WGS 84 / UTM zone 11N FALSE [y]
> 
> # does gdal reproduce what stars does, default cell size?
> (x2 = st_warp(x, crs = st_crs(x), use_gdal = FALSE))
stars object with 2 dimensions and 1 attribute
attribute(s):
               Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
geomatrix.tif    74     115    123 126.4289     132  255  224
dimension(s):
  from to  offset delta                refsys x/y
x    1 25 1840902  5.22 WGS 84 / UTM zone 11N [x]
y    1 25 1144003 -5.22 WGS 84 / UTM zone 11N [y]
> (y = setNames(st_warp(x, x2, use_gdal = TRUE, debug = FALSE, no_data_value=-1), "file.tif"))
stars object with 2 dimensions and 1 attribute
attribute(s):
          Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
file.tif    74     115    123 126.4289     132  255  224
dimension(s):
  from to  offset delta                refsys point x/y
x    1 25 1840902  5.22 WGS 84 / UTM zone 11N FALSE [x]
y    1 25 1144003 -5.22 WGS 84 / UTM zone 11N FALSE [y]
> 
> # try with multiple bands:
> tif = system.file("tif/L7_ETMs.tif", package = "stars")
> (x1 = read_stars(tif))
stars object with 3 dimensions and 1 attribute
attribute(s):
             Min. 1st Qu. Median     Mean 3rd Qu. Max.
L7_ETMs.tif     1      54     69 68.91242      86  255
dimension(s):
     from  to  offset delta                     refsys point x/y
x       1 349  288776  28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y       1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band    1   6      NA    NA                         NA    NA    
> tifcp = tempfile(fileext = ".tif")
> file.copy(tif, tifcp)
[1] TRUE
> x1p = read_stars(tifcp, proxy = TRUE)
> st_dimensions(x1p)
     from  to  offset delta                     refsys point x/y
x       1 349  288776  28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y       1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band    1   6      NA    NA                         NA    NA    
> (x1a = st_warp(x1, crs = st_crs(4326)))
stars object with 3 dimensions and 1 attribute
attribute(s):
             Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
L7_ETMs.tif     1      54     69 68.93072      86  255 8778
dimension(s):
     from  to offset      delta refsys x/y
x       1 350 -34.92  0.0002592 WGS 84 [x]
y       1 352  -7.95 -0.0002592 WGS 84 [y]
band    1   6     NA         NA     NA    
> (x1b = setNames(st_warp(x1, x1p, use_gdal = TRUE, no_data_value=-1), "file.tif"))
stars object with 3 dimensions and 1 attribute
attribute(s):
          Min. 1st Qu. Median     Mean 3rd Qu. Max.
file.tif     1      54     69 68.91242      86  255
dimension(s):
     from  to  offset delta                     refsys point x/y
x       1 349  288776  28.5 SIRGAS 2000 / UTM zone 25S FALSE [x]
y       1 352 9120761 -28.5 SIRGAS 2000 / UTM zone 25S FALSE [y]
band    1   6      NA    NA                         NA    NA    
> 
> # does gdal reproduce what stars does? Smaller grid:
> x2 = st_warp(x, crs = st_crs(x), use_gdal = FALSE, cellsize = 3)
> # x2 = x2[,2:43,2:43]
> if (interactive()) { plot(x2, breaks = "equal", axes=TRUE, reset = FALSE) }
> if (interactive()) { plot(st_as_sfc(st_bbox(x2)), add = TRUE, col = NA, border = 'red') }
> ### doesn't work: FIXME: check with more recent GDAL:
> (y = setNames(st_warp(x, x2, use_gdal = TRUE, debug = FALSE, no_data_value=-1), "file.tif"))
stars object with 2 dimensions and 1 attribute
attribute(s):
          Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
file.tif    74     107    123 126.6584     132  255  727
dimension(s):
  from to  offset delta                refsys point x/y
x    1 44 1840902     3 WGS 84 / UTM zone 11N FALSE [x]
y    1 44 1144003    -3 WGS 84 / UTM zone 11N FALSE [y]
> if (interactive()) { plot(y, breaks = "equal") }
> names(x2) = names(y)
> # isTRUE(all.equal(x2, y, check.attributes=FALSE))
> m = mean(as.vector(x2[[1]]==y[[1]]), na.rm = TRUE)
> u = unique(as.vector(x2[[1]]-y[[1]]))
> 
> g = st_as_stars()
> attr(g, "dimensions")$x$offset = 0
> g$values = as.vector(t(matrix(seq_len(prod(dim(g))), 180, 360)))
> if (interactive()) { plot(g, axes=TRUE) }
> a = st_warp(g, st_as_stars(), use_gdal=FALSE)
> b = st_warp(g, st_as_stars(), use_gdal=TRUE, no_data_value = -1)
> all.equal(a, b, check.attributes = FALSE)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   2.39    0.43    2.82