2

In the example below, I use st_apply() to scan over the space dimension of x and return a vector corresponding to the band dimension, resulting in stars object y. As x and y are in principle identical, subtracting the two should give a stars object filled with zeros. This doesn't work, as you can see in the second picture.

The problem is that the dimensions have been rearranged. How do I reshuffle them so that x-y works properly?

Mark

PS This is a rather contrived example, but it stems from a realworld example smoothing over the non-spatial dimension.

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
#> Registered S3 methods overwritten by 'stars':
#>   method             from
#>   st_bbox.SpatRaster sf  
#>   st_crs.SpatRaster  sf
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
x
#> 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 values x/y
#> x       1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
#> y       1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
#> band    1   6      NA    NA                           NA    NA   NULL
plot(x)
#> downsample set to c(1,1,1)


y <- st_apply(x,c("x","y"),FUN=function(x) x)
y
#> 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 values
#> function(x) x    1   6      NA    NA                           NA    NA   NULL
#> x                1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL
#> y                1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL
#>               x/y
#> function(x) x    
#> x             [x]
#> y             [y]
d <- x-y
plot(d)
#> downsample set to c(1,1,1)

Created on 2022-02-02 by the reprex package (v2.0.1)

Session info 
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.6.3 (2020-02-29)
#>  os       Ubuntu 20.04.3 LTS          
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       Etc/UTC                     
#>  date     2022-02-02                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date       lib source        
#>  abind       * 1.4-5   2016-07-21 [1] CRAN (R 3.6.3)
#>  assertthat    0.2.1   2019-03-21 [1] CRAN (R 3.6.3)
#>  class         7.3-15  2019-01-01 [4] CRAN (R 3.5.2)
#>  classInt      0.4-3   2020-04-07 [1] CRAN (R 3.6.3)
#>  cli           3.1.0   2021-10-27 [1] CRAN (R 3.6.3)
#>  crayon        1.4.2   2021-10-29 [1] CRAN (R 3.6.3)
#>  curl          4.3.2   2021-06-23 [1] CRAN (R 3.6.3)
#>  DBI           1.1.2   2021-12-20 [1] CRAN (R 3.6.3)
#>  digest        0.6.29  2021-12-01 [1] CRAN (R 3.6.3)
#>  dplyr         1.0.7   2021-06-18 [1] CRAN (R 3.6.3)
#>  e1071         1.7-9   2021-09-16 [1] CRAN (R 3.6.3)
#>  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 3.6.3)
#>  evaluate      0.14    2019-05-28 [1] CRAN (R 3.6.3)
#>  fansi         1.0.0   2022-01-10 [1] CRAN (R 3.6.3)
#>  fs            1.5.0   2020-07-31 [1] CRAN (R 3.6.3)
#>  generics      0.1.0   2020-10-31 [1] CRAN (R 3.6.3)
#>  glue          1.6.0   2021-12-17 [1] CRAN (R 3.6.3)
#>  highr         0.9     2021-04-16 [1] CRAN (R 3.6.3)
#>  htmltools     0.5.1.1 2021-01-22 [1] CRAN (R 3.6.3)
#>  httr          1.4.2   2020-07-20 [1] CRAN (R 3.6.3)
#>  KernSmooth    2.23-16 2019-10-15 [4] CRAN (R 3.6.1)
#>  knitr         1.33    2021-04-24 [1] CRAN (R 3.6.3)
#>  lifecycle     1.0.1   2021-09-24 [1] CRAN (R 3.6.3)
#>  lwgeom        0.2-8   2021-10-06 [1] CRAN (R 3.6.3)
#>  magrittr      2.0.1   2020-11-17 [1] CRAN (R 3.6.3)
#>  mime          0.11    2021-06-23 [1] CRAN (R 3.6.3)
#>  pillar        1.6.4   2021-10-18 [1] CRAN (R 3.6.3)
#>  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 3.6.3)
#>  proxy         0.4-26  2021-06-07 [1] CRAN (R 3.6.3)
#>  purrr         0.3.4   2020-04-17 [1] CRAN (R 3.6.3)
#>  R6            2.5.1   2021-08-19 [1] CRAN (R 3.6.3)
#>  Rcpp          1.0.7   2021-07-07 [1] CRAN (R 3.6.3)
#>  reprex        2.0.1   2021-08-05 [1] CRAN (R 3.6.3)
#>  rlang         0.4.12  2021-10-18 [1] CRAN (R 3.6.3)
#>  rmarkdown     2.10    2021-08-06 [1] CRAN (R 3.6.3)
#>  rstudioapi    0.13    2020-11-12 [1] CRAN (R 3.6.3)
#>  sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 3.6.3)
#>  sf          * 1.0-5   2021-12-17 [1] CRAN (R 3.6.3)
#>  stars       * 0.5-3   2021-06-08 [1] CRAN (R 3.6.3)
#>  stringi       1.7.3   2021-07-16 [1] CRAN (R 3.6.3)
#>  stringr       1.4.0   2019-02-10 [1] CRAN (R 3.6.3)
#>  tibble        3.1.6   2021-11-07 [1] CRAN (R 3.6.3)
#>  tidyselect    1.1.1   2021-04-30 [1] CRAN (R 3.6.3)
#>  units         0.7-2   2021-06-08 [1] CRAN (R 3.6.3)
#>  utf8          1.2.2   2021-07-24 [1] CRAN (R 3.6.3)
#>  vctrs         0.3.8   2021-04-29 [1] CRAN (R 3.6.3)
#>  withr         2.4.3   2021-11-30 [1] CRAN (R 3.6.3)
#>  xfun          0.25    2021-08-06 [1] CRAN (R 3.6.3)
#>  xml2          1.3.2   2020-04-23 [1] CRAN (R 3.6.3)
#>  yaml          2.2.1   2020-02-01 [1] CRAN (R 3.6.3)
#> 
#> [1] /home/mapa/R/x86_64-pc-linux-gnu-library/3.6
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library
Mark Payne
  • 557
  • 5
  • 12

1 Answers1

1

One solution that I have found is to use aperm() to rearrange the order of the dimensions in y before doing the arithmetic operations:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
#> Registered S3 methods overwritten by 'stars':
#>   method             from
#>   st_bbox.SpatRaster sf  
#>   st_crs.SpatRaster  sf
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
y <- st_apply(x,c("x","y"),FUN=function(x) x)
y <- aperm(y,c(2,3,1))
d <- x-y
d
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>              Min. 1st Qu. Median Mean 3rd Qu. Max.
#> L7_ETMs.tif     0       0      0    0       0    0
#> dimension(s):
#>      from  to  offset delta                       refsys point values x/y
#> x       1 349  288776  28.5 UTM Zone 25, Southern Hem... FALSE   NULL [x]
#> y       1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE   NULL [y]
#> band    1   6      NA    NA                           NA    NA   NULL
plot(d)
#> Warning in plot.stars(d): breaks="quantile" leads to a single class; maybe try
#> breaks="equal" instead?
#> downsample set to c(1,1,1)

Created on 2022-02-03 by the reprex package (v2.0.1)

Mark Payne
  • 557
  • 5
  • 12
  • Yes, we had misunderstood each other. So I deleted my answer. Anyway, the main thing is that you were able to find a solution. And, sure, `aperm()` is the right way to do what you want. Cheers – lovalery Feb 03 '22 at 11:16