Closed Shirobakaidou closed 1 year ago
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.
Could you give it a try with recent development version via:
remotes::install_github("mapme-initiative/mapme.biodiversity")
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.>
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
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 singleNA
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))