I have a dataset of ~200 observations, from scientific papers describing agroforestry, with associated coordinates (lat., lon.) dispersed all over the temperate world (US, Europe, China, etc.). I would like to extract SoilGrids' location-specific soil information (e.g. bulk density, sand, clay, SOC, nitrogen, etc.) from each of the observations – and save it in a dataframe using R (RStudio).

I thought there must be an approach where I can download extracted soil information from ONLY the observations/locations needed. Hence, avoid the need of first; downloading the big tiff files, and then; extracting the data from individual observations/locations.

Problem summarised

SoilGrids (res. 250 m) provide great data on predicted global soil properties at a high spatial resolution. The issue is that the data is stored in big (Geo)tif files that take up lots of space on the local disc when downloading and processing – an issue that becomes even more pressing when needing several soil variables.

Since, I am dealing with point locations and I don't have an area of interest (aoi) specified in a polygon it is not appropriate for me to follow the examples of Ivan Lizarazo (link: https://rpubs.com/ials2un/soilgrids_webdav) or Luis de Sousa (link: https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/wcs_from_R.md)

Attempted approaches

I have been trying the package soilDB from Beaudette et al. (2022) (link: http://ncss-tech.github.io/soilDB/index.html). The function I used to get SoilGrids data called fetchSoilGrids() is quite slow because it automatically extract ALL SoilGrids information from every location. Furthermore, the output of this function is not a data.frame nor tibble but in a class called 'SoilProfileCollection', making it hard to perform additional data processing.

Looking for an approach

Hence, I am looking for an approach where I can extract SoilGrids data directly from the ISRIC Webdav archive (link: https://files.isric.org/soilgrids/latest/data/) for ALL my observations but ONLY for a specified set of variables – returned as a data.table or similar.

The approach could for instance be a custom made function that downloads SoilGrids data by iterating through a .csv list of files (soil variables / soil layers) and a data.frame list of coordinates. The output of this function should then be stored on a local folder in a named file.

Code and data

Subset of observations

database.lat.long <-
database.lat.long %>% 
slice_sample(n = 10)


       "1",   "1",  50.6184127,          5.1146644,      "22",          "20",          "50",       "NA",           "NA",           "NA",      "CROP",       "CON",     "TILL",   "TRUE",  "FALSE",   "TRUE",  "FALSE",
      "99",   "2", 50.78118355,   4.24521861730643,      "NA",           "3",          "48",       "NA",            "4",           "47", "TREE.CROP",       "CON",     "TILL",   "TRUE",   "TRUE",   "TRUE",  "FALSE",
      "72",   "2",   50.666245,   4.55004408365023,      "NA",          "NA",          "NA",       "NA",           "NA",           "NA",        "NA",        "NA",       "NA",  "FALSE",  "FALSE",  "FALSE",  "FALSE",
     "175",  "12",  51.1638175,         10.4478313,      "28",           "4",          "83",       "NA",           "NA",           "NA", "SOIL.CROP",        "NA",       "NA",  "FALSE",   "TRUE",   "TRUE",  "FALSE",
      "45",  "15",   49.040403, -0.866531100775681,      "NA",          "NA",          "NA",       "NA",           "NA",           "NA",        "NA",       "CON",       "NA",   "TRUE",  "FALSE",  "FALSE",  "FALSE",
     "164",   "1",  55.6522136,          12.292049,      "27",          "10",          "27",       "27",           "NA",           "NA", "TREE.CROP",       "ORG",     "TILL",   "TRUE",   "TRUE",   "TRUE",   "TRUE",
     "156",   "1",  43.1147515,          1.6082177,      "NA",          "NA",          "NA",       "NA",            "5",           "11", "TREE.CROP",       "CON", "RED.TILL",  "FALSE",  "FALSE",  "FALSE",  "FALSE",
     "116",  "11",   47.560991,          7.8431299,      "60",           "1",          "60",       "NA",            "1",           "60", "TREE.CROP",       "CON",     "TILL",   "TRUE",  "FALSE",   "TRUE",  "FALSE",
     "162",   "1", 51.40866085,   9.97784694635362,       "6",          "NA",          "NA",        "6",            "1",            "6", "TREE.CROP",       "CON", "RED.TILL",  "FALSE",   "TRUE",   "TRUE",  "FALSE",
      "24",   "5",     43.7238,            4.03747,      "18",           "6",          "41",       "NA",            "6",           "41", "TREE.CROP",       "CON",     "TILL",  "FALSE",  "FALSE",   "TRUE",  "FALSE"



database.locations.df <- 
  database.lat.long %>% 
  select(ID.P, ID.S,
         latitude, longitude) %>%
  dplyr::mutate_at(vars(ID.P, ID.S), as.numeric) %>%
  rowwise() %>% 
  mutate(ID.PAIR = paste0(sort(c(ID.P, ID.S)), collapse = " ")) %>% 
  separate(ID.PAIR, " ", into = c("ID.P", "ID.S")) %>%
  mutate(ID = row_number()) %>%
  dplyr::rename("lat" = latitude,
                "lon" = longitude,
                "id" = ID) %>%
  select(id, lat, lon) %>%

x <- fetchSoilGrids(x = database.locations.df,
                    progress = TRUE)



[1] "SoilProfileCollection"
[1] "aqp"

Any ideas or suggestions to solve this problem? Also, how to convert the 'SoilProfileCollection' data type into a data.frame or similar?

  • If the files are uniformly formatted and their server supports HTTP RANGE requests or the equivalent, you could get just the portions of the files you need, but it's a long shot. Aug 25 at 14:09
  • I now found this here on the soilDB package where they explain how to convert to data.frame and Accessing, Setting, and Replacing Data ncss-tech.github.io/AQP/aqp/aqp-intro.html#1_Introduction it is quite helpful Aug 25 at 14:19
  • @KamauLindhardt Welcome to GIS.StackExchange. The question on the title is already answered below, however you make many more in your post. That is not how this site functions, please read the Tour carefully and break down your post in various questions if necessary. Aug 30 at 6:52
  • Thanks for the welcoming introduction @Luís, I will mark the answer of@Jorge as an answer to my question. Is that OK? Aug 30 at 7:13
  • Random: this could be an off-label use of Google Earth Engine Aug 30 at 12:16

1 Answer 1


There is an example for a single point queries using direct webdav access:


SoilDB is a nice package but (for what I could see), it make a full data request to the REST API, overloading the service.

Aside from the webdav you can also look into direct REST API request for point data only, if you indicate in the exact layer, depth and property the response will be relatively fast. You can find the REST API Docs here:


And then program your R code to use a generic REST API


  • Thank you @jorge I have seen this example for a single point queries using direct webdav access.. The only problem is that it needs me to specify individual soil variables (voi), depths and layers and I am interested in at least 8 of them. As I see it now it does not seem possible to specify them in for instance a vector to iterate over. Aug 26 at 8:27
  • I have also recently stumbled upon the geodata package by Robert Hijmans ( github.com/rspatial/geodata ) that seems to have the option of downloading soil data from SoilGrids via the function soil_world() Though I haven't tried it yet.. Sep 7 at 10:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.