Uninterpretable Output of the calc_indicators() function for soil properties

[Shirobakaidou]

I was trying to calculate soil indicators using the calc_indicators() function. The expected output should be three values (mean, median, sd) for each feature, but the eventual output was a single NA for all features.

My script is as follows:

library(wdpar) library(geodata) library(sf) library(ggplot2) library(mapme.biodiversity) library(tidyverse)

Download gadm of MDG

gadm <- gadm(country = "Madagascar", resolution = 1, level = 0, path = paste0(getwd(),'/data')) %>% st_as_sf()

Download wdpa

wdpa <- wdpa_fetch("Madagascar", wait = TRUE, download_dir = paste0(getwd(), '/data'))

Make Bounding Box

bbox <- st_as_sf(st_as_sfc(st_bbox(wdpa)))

Make Hexagon Grid

for (i in c(5)) {
cellarea <- i * (1e+6)
cellsize <- 2 * sqrt(cellarea / ((3 * sqrt(3) / 2))) * sqrt(3) / 2
hexa <- st_make_grid(x = bbox, cellsize = cellsize, square = FALSE)
st_write(hexa, paste0(getwd(), "/data/grids/hexa_", i, "km2.gpkg"), append= FALSE) }

Load Hexagons

hexa <- st_read(paste0(getwd(), "/data/grids/hexa_5km2.gpkg"))

Crop hexagons to Country

gadm <- st_transform(gadm, crs = st_crs(hexa))
intersection <- unlist(st_intersects(gadm, hexa))
hexa <- hexa[sort(intersection),]

Export cropped Hexagon Grids

write_sf(hexa, paste0(getwd(), "/data/grids/hexa_5km2_sub.gpkg"), append = FALSE)

Load cropped Hexagons

aoi <- st_read(paste0(getwd(), "/data/grids/hexa_5km2_sub.gpkg"))
aoi <- st_transform(aoi, "EPSG:4326")

ncores = 12
sample = FALSE
if(sample){ set.seed(123) index = sample(1:nrow(aoi), 1000, replace = FALSE) aoi = aoi[index, ] }

aoi = init_portfolio(aoi, years = 2000:2020, outdir = "./data", tmpdir = "./tmp", cores = ncores, verbose = TRUE)

Download Raster Data: Soilgrid

get_resources(aoi, resources = c("soilgrids"), layers = c("clay", "nitrogen", "phh2o", "soc"), depths = c("0-5cm", "5-15cm", "30-60cm"), stats = c("mean"))

Calculate Indicators

grid_soil <- calc_indicators(x = aoi, indicators = "soilproperties", stats_soil = c("mean", "median", "sd"), engine = "extract")
grid_soil_unnested <- grid_soil %>% unnest(soilproperties)
print(unique(grid_soil_unnested$value))

[goergen95]

Thanks, can confirm this bug. It was caused by an incorrect return class of the internal functions. I am trying to merge the changes into main today.

[goergen95]

Could you give it a try with recent development version via:

remotes::install_github("mapme-initiative/mapme.biodiversity")

[Shirobakaidou]

Could you give it a try with recent development version via:

remotes::install_github("mapme-initiative/mapme.biodiversity")

After removing the old version and installing the newest version, another error message occurred when I was trying to run get_resources(aoi, resources = c("soilgrids"), layers = c("clay", "nitrogen", "phh2o", "soc"), depths = c("0-5cm", "5-15cm", "30-60cm"), stats = c("mean"))

The error message is called

<error/vctrs_error_scalar_type>
Error in stop_vctrs():
! Input must be a vector, not an environment.

Backtrace:
1. mapme.biodiversity::get_resources(...)
16. vctrs:::stop_scalar_type(<fn>(<env>), "")
17. vctrs:::stop_vctrs(msg, "vctrs_error_scalar_type", actual = x)
<simpleWarning: Download for resource soilgrids failed. Returning unmodified portfolio object.>

[goergen95]

The cause of this might be related to issues with your data (similar to #117). Using the same fix for the CRS I see the following:

library(geodata)
#> Loading required package: terra
#> terra 1.6.47
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE
remotes::install_github("mapme-initiative/mapme.biodiversity")
#> Skipping install of 'mapme.biodiversity' from a github remote, the SHA1 (a897bef4) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(mapme.biodiversity)
library(tidyverse)

data_dir <- "./issues"
dir.create(data_dir, showWarnings = FALSE)
dir.create(file.path(data_dir, "data"), recursive = T, showWarnings = FALSE)
dir.create(file.path(data_dir, "tmp"), recursive = T, showWarnings = FALSE)

gadm <- gadm(country = "Madagascar",
             resolution = 1,
             level = 0,
             path = file.path(data_dir, "data")) %>% st_as_sf()

aoi <- st_read("hexa_5km2_sub.gpkg")
#> Reading layer `hexa_5km2_sub' from data source 
#>   `/home/darius/Desktop/mapme/hexa_5km2_sub.gpkg' using driver `GPKG'
#> Simple feature collection with 120076 features and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 4165767 ymin: -3162891 xmax: 4873395 ymax: -1537018
#> CRS:           unknown
st_crs(aoi) <- st_crs(6933) # fix CRS
aoi <- st_transform(aoi, "EPSG:4326")
aoi$id <- 1:nrow(aoi)

# restrict analysis to polygons off of Madagascar mainland
mdg <- st_cast(gadm, "POLYGON")[1,]
#> Warning in st_cast.sf(gadm, "POLYGON"): repeating attributes for all
#> sub-geometries for which they may not be constant
intersection <- st_filter(aoi, mdg)
aoi <- aoi[!aoi$id %in% intersection$id, ]

ncores = 6
aoi = init_portfolio(aoi,
                     years = 2000:2020,
                     outdir = file.path(data_dir, "data"),
                     tmpdir = file.path(data_dir, "tmp"),
                     cores = ncores,
                     verbose = TRUE)

aoi <- get_resources(aoi,
                     resources = c("soilgrids"),
                     layers = c("clay"),
                     depths = c("0-5cm"),
                     stats = c("mean"))
#> Starting process to download resource 'soilgrids'........
#> Starting to download data for layer 'clay', depth '0-5cm', and stat 'mean'. This may take a while...

grid_soil <- calc_indicators(aoi,
                             indicators = "soilproperties",
                             stats_soil = c("mean", "median", "sd"),
                             engine = "exactextract")

grid_soil %>%
  unnest(soilproperties) |>
  pull(mean) |>
  unique() |>
  head()
#> [1]      NaN 27.23902 26.55611 30.30000 25.25556 25.57976

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