Reclassify 2D stars array

With the raster package we can reclassify two dimensional spatial arrays as follows:

set.seed(1)

library(raster)
#> Loading required package: sp

r <- raster::raster(matrix(sample(1:9, 100, replace = TRUE), 10, 10))

raster::plot(r)

rcl_matrix <- as.matrix(data.frame(from = c(1,4,7), to = c(3,6,9), becomes = 1:3))

raster::plot(raster::reclassify(r, rcl = rcl_matrix, right = NA))

^{Created on 2022-01-07 by the reprex package (v2.0.1)}

I struggle to do the same operation using the stars package. It is crucial that I do not have to specify the slice name as the array I am doing the operation for has always two spatial dimensions but the name can vary. In my actual applications there can be tens of different from-to classes. Specifying each one manually would not work. A stars vignette instructs to do the operation using dplyr::case_when() or forcats::fct_recode(). I only manage to do this a clumsy way:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

rcl_matrix <- as.matrix(data.frame(from = c(1,4,7), to = c(3,6,9), becomes = 1:3))
sr <- stars::st_as_stars(matrix(sample(1:9, 100, replace = TRUE), 10, 10))

plot(dplyr::mutate(sr, A1 = dplyr::case_when(A1 < 4 ~ 1, A1 < 7 ~ 2, A1 < 10 ~ 3)))

# OR:

sr$A1[sr$A1 <= rcl_matrix[1,2]] <- rcl_matrix[1,3]
sr$A1[sr$A1 <= rcl_matrix[2,2] & sr$A1 > rcl_matrix[1,2]] <- rcl_matrix[2,3]
sr$A1[sr$A1 <= rcl_matrix[3,2] & sr$A1 > rcl_matrix[1,2]] <- rcl_matrix[3,3]

plot(sr)

^{Created on 2022-01-07 by the reprex package (v2.0.1)}

How to do the reclassification for stars objects as elegantly/easily as for raster objects?

Solution 1:

Here is one way of doing it for numeric rasters. Note that you'll have to set the first break lower than in the raster::reclassify example.

set.seed(1)
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

rcl <- data.frame(from = c(1,4,7), to = c(3,6,9), becomes = 1:3)
sr <- stars::st_as_stars(matrix(sample(1:9, 100, replace = TRUE), 10, 10))

plot(cut(sr, c(0, rcl$to), labels = rcl$becomes))

^{Created on 2022-01-13 by the reprex package (v2.0.1)}

Different patterns in the outcomes stem from random sampling. The outcomes, in this case, are identical except that raster::reclassify returns numeric, while stars::cut returns factor:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(raster)
#> Loading required package: sp

dt <- matrix(rep(1:10, 10), 10, 10)
rcl <- data.frame(from = c(1,4,7), to = c(3,6,10), becomes = 1:3)

r <- raster::raster(dt)
sr <- stars::st_as_stars(dt)

r_rec <- raster::reclassify(r, rcl = rcl, right = NA)
sr_cut <- cut(sr, c(0, rcl$to), labels = rcl$becomes)

all.equal(as.matrix(r_rec), 
          matrix(as.numeric(as.character(unlist(sr_cut[[1]]))),nrow=nrow(sr_cut[[1]]))
)
#> [1] TRUE

^{Created on 2022-01-13 by the reprex package (v2.0.1)}