Stacking rasters for small AOIs throws error

[carostolle]

Hi,

we are trying to load a stars object that has a DEM - retrieved through the elevatr package - as well as some terrain metrics calculated via whitebox as attributes. However, for small areas the combining of the two separate stars objects fails with the error: Error in c.stars(structure(list(elevation2edc4a34fdfd.tif = structure(c(328.039428710938, : don't know how to merge arrays: please specify parameter along The same code for a bigger area runs without problem, see in the example below. When using terra to stack the rasters it works for both AOIs. After downgrading to stars v.0.5-5 the code also works for the smaller AOI.

Thanks for any hints on what the issue could be!

library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE
library(stars)
#> Loading required package: abind
library(elevatr)
library(whitebox)

temp_files = c("elevation", "slope") |>
  vapply(
    \(f) tempfile(pattern = f, fileext = ".tif")
    , FUN.VALUE = character(1L)
  )

on.exit(unlink(temp_files), add = TRUE, after = FALSE)

## AOI small
df = data.frame(
  lon = c(11.15525, 11.15622)
  , lat =  c(49.76438, 49.76510)
)

## AOI big
# df = data.frame(
#     lon = c(11.1127, 11.20474)
#     , lat =  c(49.74578, 49.7775)
# )

aoi = st_as_sf(
  df
  , coords = c("lon", "lat")
  , crs = 4326
) |>
  st_bbox() |>
  st_as_sfc() |>
  st_as_sf()

dem = get_elev_raster(
  st_transform(aoi, crs = 3035)
  , z = 14
  , clip = "locations"
)
#> Mosaicing & Projecting
#> Clipping DEM to locations
#> Note: Elevation units are in meters.

dem = stars::st_as_stars(dem)
dem = setNames(dem, nm = "elevation")

## slope
stars::write_stars(dem, dsn = temp_files[["elevation"]])

# calculate slope stars
whitebox::wbt_slope(
  temp_files[["elevation"]]
  , output = temp_files[["slope"]]
)

slope = stars::read_stars(temp_files[["slope"]])
dem = stars::read_stars(temp_files[["elevation"]])

# dimensions
all.equal(
  st_dimensions(slope)
  , current = st_dimensions(dem)
)
#> [1] "Attributes: < Component \"raster\": Component \"blocksizes\": Mean relative difference: 25 >"

stars = stars::read_stars(temp_files)
#> Error in c.stars(structure(list(elevation2edc4a34fdfd.tif = structure(c(328.039428710938, : don't know how to merge arrays: please specify parameter along

# stack using terra
library(terra)
#> terra 1.6.7

slope_t = terra::rast(temp_files[["slope"]])
dem_t = terra::rast(temp_files[["elevation"]])

# dimensions
all.equal(
  dim(slope_t)
  , current = dim(dem_t)
)
#> [1] TRUE

stack = terra::rast(temp_files)

^{Created on 2022-09-09 with reprex v2.0.2}

[edzer]

Could you provide a reprex that doesn't use whitebox? Installing that from CRAN leads to an error,

sh: 1: /home/edzer/R/x86_64-pc-linux-gnu-library/4.0/whitebox/WBT/whitebox_tools: not found

on my side, suggesting that the CRAN package is not enough.

[carostolle]

Hi thanks for the quick reply. Usually whitebox should work after running whitebox::install_whitebox(), but I will try to figure out a reprex without whitebox.

[edzer]

Thanks, that works!

[edzer]

Great report - should work now!

[carostolle]

Yes, that seems to have done the trick! Thank you!

[alexyshr]

Hi @carostolle

I have installed elevatr, sf and stars from GitHub, and also whitebox::install_whitebox(). I get next error when running get_elev_layer to create variable dem

Accessing raster elevation [=========================] 100%
Mosaicing & Projecting
Clipping DEM to locations
Error: OGRCreateCoordinateTransformation(): transformation not available
In addition: Warning message:
In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot find coordinate operations from `EPSG:3035' to `ENGCRS["ETRS89-extended / LAEA Europe",EDATUM["Unknown engineering datum"],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

Any suggestions to test your code?

Thanks!

[alexyshr]

It works in crs epsg 4326!

[edzer]

@alexyshr can you provide a reprex of the thing that does not work?

[alexyshr]

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(stars)
#> Loading required package: abind

library(elevatr)
#> elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
#> of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
#> deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
#> will be dropped in future versions, so please plan accordingly.
library(whitebox)

temp_files <- c("elevation", "slope") |>
  vapply(
    \(f) tempfile(pattern = f, fileext = ".tif"),
    FUN.VALUE = character(1L)
  )
on.exit(unlink(temp_files), add = TRUE, after = FALSE)

# AOI small
df <- data.frame(
  lon = c(11.15525, 11.15622),
  lat = c(49.76438, 49.76510)
)

## AOI big
# df = data.frame(
#     lon = c(11.1127, 11.20474)
#     , lat =  c(49.74578, 49.7775)
# )

aoi <- st_as_sf(
  df,
  coords = c("lon", "lat"),
  crs = 4326
) |>
  st_bbox() |>
  st_as_sfc() |>
  st_as_sf()

#without trasformin to 3035 it works
aoi <- st_transform(aoi, crs = 3035)

#the issue is in elevatr::get_elev_raster
dem <- get_elev_raster(
  aoi,
  z = 14,
  clip = "locations"
)
#> Mosaicing & Projecting
#> Clipping DEM to locations
#> Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
#> GDAL Error 6: Cannot find coordinate operations from `EPSG:3035' to
#> `ENGCRS["ETRS89-extended / LAEA Europe",EDATUM["Unknown engineering
#> datum"],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
#> Error in eval(expr, envir, enclos): OGRCreateCoordinateTransformation(): transformation not available
# #Mosaicing & Projecting
# #Clipping DEM to locations
# #Note: Elevation units are in meters.
# # it returns a raster layer

# dem <- stars::st_as_stars(dem)
# dem <- setNames(dem, nm = "elevation")

# ## slope
# stars::write_stars(dem, dsn = temp_files[["elevation"]])

# # calculate slope stars
# whitebox::wbt_slope(
#   temp_files[["elevation"]],
#   output = temp_files[["slope"]]
# )

# slope <- stars::read_stars(temp_files[["slope"]])
# dem <- stars::read_stars(temp_files[["elevation"]])

# # dimensions
# all.equal(
#   st_dimensions(slope),
#   current = st_dimensions(dem)
# )

# stars <- stars::read_stars(temp_files)

# # stack using terra
# library(terra)

# slope_t <- terra::rast(temp_files[["slope"]])
# dem_t <- terra::rast(temp_files[["elevation"]])

# # dimensions
# all.equal(
#   dim(slope_t),
#   current = dim(dem_t)
# )

# stack <- terra::rast(temp_files)

^{Created on 2024-04-07 with reprex v2.0.2}

Session info

``` r sessioninfo::session_info() #> - Session info --------------------------------------------------------------- #> setting value #> version R version 4.3.3 (2024-02-29 ucrt) #> os Windows 10 x64 (build 19045) #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate English_United States.utf8 #> ctype English_United States.utf8 #> tz America/New_York #> date 2024-04-07 #> pandoc 2.17.1.1 @ C:/Users/500596~1/AppData/Local/Pandoc/ (via rmarkdown) #> #> - Packages ------------------------------------------------------------------- #> package * version date (UTC) lib source #> abind * 1.4-5 2016-07-21 [1] CRAN (R 4.3.0) #> class 7.3-22 2023-05-03 [2] CRAN (R 4.3.3) #> classInt 0.4-10 2023-09-05 [1] CRAN (R 4.3.1) #> cli 3.6.2 2023-12-11 [2] CRAN (R 4.3.3) #> codetools 0.2-19 2023-02-01 [2] CRAN (R 4.3.3) #> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.1) #> curl 5.2.1 2024-03-01 [2] CRAN (R 4.3.3) #> DBI 1.2.2 2024-02-16 [1] CRAN (R 4.3.3) #> digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.1) #> e1071 1.7-14 2023-12-06 [2] CRAN (R 4.3.3) #> elevatr * 1.0.0.9999 2024-04-07 [1] Github (jhollist/elevatr@041d88e) #> evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.1) #> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.1) #> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.1) #> glue 1.7.0 2024-01-09 [2] CRAN (R 4.3.3) #> hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.1) #> htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.3.1) #> httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.1) #> KernSmooth 2.23-22 2023-07-10 [1] CRAN (R 4.3.1) #> knitr 1.44 2023-09-11 [1] CRAN (R 4.3.1) #> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.3) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.1) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.1) #> prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.3.1) #> progress 1.2.3 2023-12-06 [1] CRAN (R 4.3.3) #> progressr 0.14.0 2023-08-10 [1] CRAN (R 4.3.1) #> proxy 0.4-27 2022-06-09 [1] CRAN (R 4.3.1) #> purrr 1.0.2 2023-08-10 [1] CRAN (R 4.3.1) #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.3.2) #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.3.1) #> R.oo 1.26.0 2024-01-24 [1] CRAN (R 4.3.2) #> R.utils 2.12.3 2023-11-18 [1] CRAN (R 4.3.2) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.1) #> Rcpp 1.0.12 2024-01-09 [1] CRAN (R 4.3.3) #> reprex 2.0.2 2022-08-17 [1] CRAN (R 4.3.1) #> rlang 1.1.3 2024-01-10 [2] CRAN (R 4.3.3) #> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.1) #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.1) #> sf * 1.0-17 2024-04-06 [1] Github (r-spatial/sf@2867dd2) #> slippymath 0.3.1 2019-06-28 [1] CRAN (R 4.3.3) #> stars * 0.6-5 2024-04-07 [1] Github (r-spatial/stars@c645e77) #> styler 1.10.2 2023-08-29 [1] CRAN (R 4.3.2) #> terra 1.7-74 2024-04-07 [1] Github (rspatial/terra@29fe968) #> units 0.8-5 2023-11-28 [2] CRAN (R 4.3.3) #> vctrs 0.6.5 2023-12-01 [2] CRAN (R 4.3.3) #> whitebox * 2.3.4 2023-11-18 [1] CRAN (R 4.3.3) #> withr 3.0.0 2024-01-16 [1] CRAN (R 4.3.3) #> xfun 0.40 2023-08-09 [1] CRAN (R 4.3.1) #> yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0) #> #> [1] C:/Users/500596972/AppData/Local/R/win-library/4.3 #> [2] C:/Program Files/R/R-4.3.3/library #> #> ------------------------------------------------------------------------------ ```

[edzer]

That works for me:

[alexyshr]

Thank you!

I will install all those packages from the source again. Any other advice will be appreciated.

Best!