better coercion from SpatRaster #484
Conversation
|
Thanks, Robert! Yes, I can see how For now, I guess if |
|
I think this nails it pretty much. Do you have an example with time series data, preferably with different amount of time instances per source? The attribute name of the stars object is still odd (same as first band name); what should it be? library(terra)
# terra version 1.5.2
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
s <- rast(system.file("ex/logo.tif", package="terra"))
x <- c(s[[2]], s[[3:1]])
names(x) <- letters[1:4]
### x has two sources (which happen to be the same file) but never the bands as in the file
x
# class : SpatRaster
# dimensions : 77, 101, 4 (nrow, ncol, nlyr)
# resolution : 1, 1 (x, y)
# extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
# sources : logo.tif
# logo.tif (3 layers)
# names : a, b, c, d
# min values : 0, 0, 0, 0
# max values : 255, 255, 255, 255
#class : SpatRaster
#dimensions : 77, 101, 4 (nrow, ncol, nlyr)
#resolution : 1, 1 (x, y)
#extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#sources : logo.tif
# logo.tif (3 layers)
#names : a, b, c, d
#min values : 0, 0, 0, 0
#max values : 255, 255, 255, 255
### note the two values for sid and the band numbers
sources(x, bands=TRUE)
# sid source bands
# 1 1 /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/terra/ex/logo.tif 2
# 2 2 /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/terra/ex/logo.tif 3
# 3 2 /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/terra/ex/logo.tif 2
# 4 2 /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/terra/ex/logo.tif 1
# sid source bands
#1 1 C:/soft/R/R-4.1.2/library/terra/ex/logo.tif 2
#2 2 C:/soft/R/R-4.1.2/library/terra/ex/logo.tif 3
#3 2 C:/soft/R/R-4.1.2/library/terra/ex/logo.tif 2
#4 2 C:/soft/R/R-4.1.2/library/terra/ex/logo.tif 1
### the result of this is not correct (and has a warning)
st_as_stars(x)
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# a 0 138 204 186.448 254 255
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 101 0 1 unnamed NA NULL [x]
# y 1 77 77 -1 unnamed NA NULL [y]
# band 1 4 NA NA NA NA a,...,d
st_as_stars(s)
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# red 0 138 206 186.8136 254 255
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 101 0 1 unnamed NA NULL [x]
# y 1 77 77 -1 unnamed NA NULL [y]
# band 1 3 NA NA NA NA red , green, blue |
|
Thanks, perhaps something like this then: There is no particular use case for allowing combining (subsets of) different sources, it is a general usability feature. |
|
What you suggested was actually more like this, and probably better |
|
Although stored by source, I consider "time" a property of the SpatRaster, with each layer having a single timestamp. That can be a POSIX or a Date or "raw" (just a number). Either all layers have a timestamp, or they are all NA. But they may not be consecutive. |
|
great, we now have library(terra)
# terra version 1.5.2
s <- rast(system.file("ex/logo.tif", package="terra"))
time(s) <- as.Date("2001-05-04") + 0:2
ss <- c(s, s)
time(ss) <- as.Date("2001-05-04") + 0:5
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
(x <- c(ss[[1:3]], ss[[1:2]]))
# class : SpatRaster
# dimensions : 77, 101, 5 (nrow, ncol, nlyr)
# resolution : 1, 1 (x, y)
# extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
# sources : logo.tif (3 layers)
# logo.tif (2 layers)
# names : red, green, blue, red, green
# min values : 0, 0, 0, 0, 0
# max values : 255, 255, 255, 255, 255
# time : 2001-05-04 to 2001-05-06
st_as_stars(x)
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# red 0 136 203 185.6155 254 255
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 101 0 1 unnamed NA NULL [x]
# y 1 77 77 -1 unnamed NA NULL [y]
# time 1 5 NA NA Date NA 2001-05-04,...,2001-05-05 |
|
(still a funny attribute name: |
|
I am not sure, but I take it that an attribute is a (sub)-dataset? If so, there is no direct match, as a SpatRaster always is a single dataset. So perhaps just use a generic name like "data"? |
|
Actually, there are also Also, any chance of getting this to CRAN soon? It would help as I could then update terra for UCRT. |
submitted! |
|
On CRAN now... |
|
There still seems to be some friction when converting between
|
|
|
|
@edzer, thanks. |
|
That might work now. |
|
@edzer, I updated |
|
OK, can you please prepare a bit more of a reprex for this? |
|
This is what I get: library(raster) #‘3.5.19’
# Loading required package: sp
library(terra)
# terra 1.5.21
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
packageVersion("sf")
# [1] ‘1.0.7’
r<- rast(nrow=100, ncol=100, nlyr=2)
set.seed(5)
values(r[[1]])<- rnorm(n = ncell(r), mean = 20, sd = 3)
set.seed(5)
values(r[[2]])<- sample(1:7, size = ncell(r), replace = TRUE)
r[[2]]<- as.factor(r[[2]])
levels(r[[2]])<- data.frame(ID=1:7, features = c("c1","c2", "c3", "c4", "c5", "c6", "c7"))
names(r[[2]])<- "features"
droplevels(st_as_stars(r))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# values 1 4 7.275819 12.01187 19.98289 30.81064
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 100 -180 3.6 WGS 84 NA NULL [x]
# y 1 100 90 -1.8 WGS 84 NA NULL [y]
# band 1 2 NA NA NA NA lyr.1 , features
droplevels(st_as_stars(r[[2]]))
# stars object with 2 dimensions and 1 attribute
# attribute(s):
# values
# c1:1429
# c2:1381
# c3:1441
# c4:1443
# c5:1408
# c6:1429
# c7:1469
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 100 -180 3.6 WGS 84 NA NULL [x]
# y 1 100 90 -1.8 WGS 84 NA NULL [y] |
|
It does seem
|
|
A simpler way to look at my post directly above is that the number of levels is set to 1-255 when reading in the file with |
|
It looks like |
|
That should work better now. |
|
@edzer, it seems to be mostly there now. This edge case however is still not accounted for: if classes are defined starting at 0 instead of 1, but that class 0 class is not present in the raster. When converting from a |
|
Thanks for pushing the edge case! This gives library(terra)
# terra 1.5.40
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is TRUE
library(raster)
# Loading required package: sp
# Warning message:
# no function found corresponding to methods exports from ‘raster’ for: ‘area’
r<- rast(nrow=10, ncol=10)
set.seed(5)
values(r)<- sample(1:2, size = ncell(r), replace = TRUE) #3 classes but only 2 are present (class 0 missing)
values(r)[1]<- NA
levels(r)<- data.frame(ID=0:2, category = c("c1", "c2", "c3")) #Define 3 classes starting from 0
freq(r) #raster only has classes c2 and c3
# layer value count
# 1 1 c2 48
# 2 1 c3 51
r_stars<- st_as_stars(r)
unique(r_stars["values"][[1]]) #stars says that raster instead is composed of class c1 and c2
# [1] <NA> c2 c3
# Levels: c2 c3
writeRaster(r, "rnew.tif", overwrite=TRUE)
levels(terra::rast("rnew.tif"))[[1]]
# value category
# 1 0 c1
# 2 1 c2
# 3 2 c3
unique(stars::read_stars("rnew.tif")[["rnew.tif"]]) # when read in from file stars loses a factor level
# [1] <NA> c2 c3
# Levels: c2 c3This loses factor level |
|
Thanks @edzer! I still think it'd be best if all factor levels are retained since On a related note, if the final class (e.g. c3 in this example) is missing, conversion from |
|
Good catch; that should work now. |
|
@edzer thanks. The classes are now mapped correctly and there is no error. There is some slight inconsistency on how things are handled though as in some cases the first/last class will be dropped if it is not present and in others they will be kept. I think keeping all classes would be the optimal behavior as then the first and last factor levels are not treated as "special" (i.e. any class in the middle will be present in the levels whether or not it is present). |
* handle 0 and gaps in attribute and color tables
|
I think this pretty much nails it: zeroes, and gaps in attribute / color tables, and roundtripping. Note that this (currently) needs the respective branches in library(terra)
# terra 1.5.43
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
# Class 1 missing
r0 <- rast(nrow=2, ncol=2)
values(r0) <- 0:3
values(r0)[1] <- NA
r1 <- r0 + 1
# 1: all classes
levels(r0) <- data.frame(ID=0:3, category = c("c0", "c1", "c2", "c3")) #Define 3 classes starting from 0
levels(r1) <- data.frame(ID=1:4, category = c("c1", "c2", "c3", "c4")) #Define 3 classes starting from 1
writeRaster(r0, "r0.tif", overwrite=TRUE)
writeRaster(r1, "r1.tif", overwrite=TRUE)
all.equal(levels(stars::read_stars("r0.tif")[["r0.tif"]]),
levels(st_as_stars(r0)[["values"]]))
# [1] TRUE
all.equal(levels(stars::read_stars("r1.tif")[["r1.tif"]]),
levels(st_as_stars(r1)[["values"]]))
# [1] TRUE
lc = read_stars(system.file("tif/lc.tif", package = "stars"))
write_stars(lc, "lc.tif")
lc2 = read_stars("lc.tif")
lc
# stars object with 2 dimensions and 1 attribute
# attribute(s):
# lc.tif
# Evergreen Forest : 456
# Herbaceuous : 270
# Open Water : 252
# Developed, Low Intensity : 81
# Developed, Medium Intensity: 48
# (Other) : 142
# NA's :2615
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 84 3092415 3000 Albers Conical Equal Area FALSE NULL [x]
# y 1 46 59415 -3000 Albers Conical Equal Area FALSE NULL [y]
lc2
# stars object with 2 dimensions and 1 attribute
# attribute(s):
# lc.tif
# Evergreen Forest : 456
# Herbaceuous : 270
# Open Water : 252
# Developed, Low Intensity : 81
# Developed, Medium Intensity: 48
# (Other) : 142
# NA's :2615
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 84 3092415 3000 Albers Conical Equal Area FALSE NULL [x]
# y 1 46 59415 -3000 Albers Conical Equal Area FALSE NULL [y]
all.equal(lc, lc2)
# [1] TRUE |
|
Hmmn, @edzer I have the latest development version of remotes::install_github("r-spatial/sf")
#> Skipping install of 'sf' from a github remote, the SHA1 (4f550e55) has not changed since last install.
#> Use `force = TRUE` to force installation
remotes::install_github("r-spatial/stars")
#> Skipping install of 'stars' from a github remote, the SHA1 (78b96b92) has not changed since last install.
#> Use `force = TRUE` to force installation
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
lc = read_stars(system.file("tif/lc.tif", package = "stars"))
#> Warning in read_stars(system.file("tif/lc.tif", package = "stars")): categorical
#> data values starting at 0 are shifted with one to start at 1
write_stars(lc, "lc.tif")
lc2 = read_stars("lc.tif")
lc
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#> lc.tif
#> :2615
#> Evergreen Forest : 456
#> Herbaceuous : 270
#> Open Water : 252
#> Developed, Low Intensity : 81
#> Developed, Medium Intensity: 48
#> (Other) : 142
#> dimension(s):
#> from to offset delta refsys point values x/y
#> x 1 84 3092415 3000 Albers Conical Equal Area FALSE NULL [x]
#> y 1 46 59415 -3000 Albers Conical Equal Area FALSE NULL [y]
lc2
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#> lc.tif
#> :2615
#> Evergreen Forest : 456
#> Herbaceuous : 270
#> Open Water : 252
#> Developed, Low Intensity : 81
#> Developed, Medium Intensity: 48
#> (Other) : 142
#> dimension(s):
#> from to offset delta refsys point values x/y
#> x 1 84 3092415 3000 Albers Conical Equal Area FALSE NULL [x]
#> y 1 46 59415 -3000 Albers Conical Equal Area FALSE NULL [y]
all.equal(lc, lc2)
#> [1] "Component \"lc.tif\": Attributes: < Component \"colors\": Lengths (256, 255) differ (string compare on first 255) >"Created on 2022-06-17 by the reprex package (v2.0.1) |
|
You need the branches: |
|
Thanks @edzer! Seems to be working great! |
|
Categories are ignored for later layers though if there's more than one categorical raster. library(terra)
#> terra 1.5.43
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
r<- rast(nrow=10, ncol=10)
set.seed(5)
values(r)<- sample(1:3, size = ncell(r), replace = TRUE)
levels(r) <- data.frame(ID=1:3, category = c("c1", "c2", "c3"))
r2<- rast(nrow=10, ncol=10)
set.seed(6)
values(r2)<- sample(1:3, size = ncell(r2), replace = TRUE)
levels(r2) <- data.frame(ID=1:3, category = c("x1", "x2", "x3"))
r3<- c(r,r2)
names(r3)<- c("lyr.1", "lyr.2")
stars_r3<- st_as_stars(r3)
#> Warning in st_as_stars.SpatRaster(r3): ignoring categories/levels for all but
#> first layerCreated on 2022-06-17 by the reprex package (v2.0.1) |
|
That is by design, as currently bands in a > c(st_as_stars(r), st_as_stars(r2))
stars object with 2 dimensions and 2 attributes
attribute(s):
values values.1
c1:34 x1:28
c2:32 x2:37
c3:34 x3:35
dimension(s):
from to offset delta refsys point values x/y
x 1 10 -180 36 WGS 84 NA NULL [x]
y 1 10 90 -18 WGS 84 NA NULL [y]and we could, in principle, make it an option for |
|
Not sure where the library(terra)
# terra 1.5.43
#> terra 1.5.43
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
r<- rast(nrow=10, ncol=10)
set.seed(5)
values(r)<- sample(1:3, size = ncell(r), replace = TRUE)
levels(r) <- data.frame(ID=1:3, category = c("c1", "c2", "c3"))
r2<- rast(nrow=10, ncol=10)
set.seed(6)
values(r2)<- sample(1:3, size = ncell(r2), replace = TRUE)
levels(r2) <- data.frame(ID=1:3, category = c("x1", "x2", "x3"))
r3<- c(r,r2)
names(r3)<- c("lyr.1", "lyr.2")
st_as_stars(r3)
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in (function (x) : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in x$.self$finalize() : attempt to apply non-function
# Error in (function (x) : attempt to apply non-function
# stars object with 2 dimensions and 2 attributes
# attribute(s):
# lyr.1 lyr.2
# c1:34 x1:28
# c2:32 x2:37
# c3:34 x3:35
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 10 -180 36 WGS 84 NA NULL [x]
# y 1 10 90 -18 WGS 84 NA NULL [y]
(a = c(st_as_stars(r3[[1]]), st_as_stars(r3[[2]])))
# stars object with 2 dimensions and 2 attributes
# attribute(s):
# lyr.1 lyr.2
# c1:34 x1:28
# c2:32 x2:37
# c3:34 x3:35
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 10 -180 36 WGS 84 NA NULL [x]
# y 1 10 90 -18 WGS 84 NA NULL [y]
(b = st_redimension(a))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# lyr.1.lyr.2
# c1:34
# c2:32
# c3:34
# x1:28
# x2:37
# x3:35
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 10 -180 36 WGS 84 NA NULL [x]
# y 1 10 90 -18 WGS 84 NA NULL [y]
# new_dim 1 2 NA NA NA NA lyr.1, lyr.2
(c = split(b))
# stars object with 2 dimensions and 2 attributes
# attribute(s):
# lyr.1 lyr.2
# c1:34 c1: 0
# c2:32 c2: 0
# c3:34 c3: 0
# x1: 0 x1:28
# x2: 0 x2:37
# x3: 0 x3:35
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 10 -180 36 WGS 84 NA NULL [x]
# y 1 10 90 -18 WGS 84 NA NULL [y]
(d = droplevels(c))
# stars object with 2 dimensions and 2 attributes
# attribute(s):
# lyr.1 lyr.2
# c1:34 x1:28
# c2:32 x2:37
# c3:34 x3:35
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 10 -180 36 WGS 84 NA NULL [x]
# y 1 10 90 -18 WGS 84 NA NULL [y] |
|
@edzer awesome! The info seems to get transferred properly. library(terra)
#> terra 1.5.44
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(tmap)
r<- rast(nrow=10, ncol=10)
set.seed(5)
values(r)<- sample(1:3, size = ncell(r), replace = TRUE)
levels(r) <- data.frame(ID=1:3, category = c("c1", "c2", "c3"))
r2<- rast(nrow=10, ncol=10)
set.seed(6)
values(r2)<- sample(1:3, size = ncell(r2), replace = TRUE)
levels(r2) <- data.frame(ID=1:3, category = c("x1", "x2", "x3"))
r3<- c(r,r2)
names(r3)<- c("lyr.1", "lyr.2")
r3_stars<- st_as_stars(r3)
plot(r3_stars)plot(r3_stars["lyr.2"])tm_shape(r3_stars)+
tm_raster()Created on 2022-06-19 by the reprex package (v2.0.1) |
It does, but it doesn't plot secondary attributes; the multi-band version assumes the array constitutes single factor: plot(st_redimension(r3_stars))
|
|
Thanks! I still need to get a bit more familiar with |
The motivation for this pull request is two fold. First, an upcoming change to
terra::sourcesbreaks the currentstars::st_as_stars; with the proposed changes it works with the CRAN and the development version. Second, the coercion from a file source ignored the bands; it always assumed that all bands were included. It also ignores the fact that there can be multiple source files (and perhaps the same source file multiple times). With the current version, you getThis pull request changes the behavior to taking the first file and the bands used from that file:
That is still suboptimal, but not wrong especially if a warning is given about dropped layers.
Better is possible, and I commented a bit more in the code, but I wanted to propose this small step first