forked from r-spatial/stars
-
Notifications
You must be signed in to change notification settings - Fork 0
/
warp.R
37 lines (33 loc) · 1.31 KB
/
warp.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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)
y = st_warp(x, x2)
plot(x2, breaks = "equal", axes=TRUE)
plot(y, breaks = "equal", axes=TRUE)
names(x2) = names(y)
all.equal(x2, y) # yes?
if (FALSE) { # FIXME: to be removed when sf >= 0.7-4 is on CRAN
# does gdal reproduce with stars template object?
(x2 = st_warp(x, y, use_gdal = TRUE))
# does gdal reproduce what stars does, default cell size?
(x2 = st_warp(x, crs = st_crs(x), use_gdal = FALSE))
(y = st_warp(x, x2, use_gdal = TRUE, debug = FALSE))
# try with multiple bands:
tif = system.file("tif/L7_ETMs.tif", package = "stars")
(x1 = read_stars(tif))
(x1p = read_stars(tif, proxy = TRUE))
(x1a = st_warp(x1, crs = st_crs(4326)))
(x1b = st_warp(x1, x1p, use_gdal = TRUE))
# 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]
plot(x2, breaks = "equal", axes=TRUE, reset = FALSE)
plot(st_as_sfc(st_bbox(x2)), add = TRUE, col = NA, border = 'red')
### doesn't work: FIXME: check with more recent GDAL:
#(y = st_warp(x, x2, use_gdal = TRUE, debug = FALSE))
#plot(y, breaks = "equal")
#names(x2) = names(y)
#all.equal(x2, y)
}