R Under development (unstable) (2024-03-24 r86185 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)) stars object with 2 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA's geomatrix.tif 74 107 123 126.6857 132 255 727 dimension(s): from to offset delta refsys x/y x 1 44 1840902 3 WGS 84 / UTM zone 11N [x] y 1 44 1144003 -3 WGS 84 / UTM zone 11N [y] > # 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)) > mean(as.vector(x2[[1]]==y[[1]]), na.rm = TRUE) [1] 0.9958644 > unique(as.vector(x2[[1]]-y[[1]])) [1] NA 0 -33 -25 25 49 17 > > 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.43 0.43 2.81