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, 2022 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, 2022 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, 2022 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, 2022 at 7:13
  • Random: this could be an off-label use of Google Earth Engine Aug 30, 2022 at 12:16

2 Answers 2


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, 2022 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, 2022 at 10:24

As your comment mentions, you can use geodata::soil_world which permits download of SoilGrids data without going through the buggy SoilGrids R tutorial. I pair geodata and terra and avoid use of raster and other to-be outdated packages.

UPDATE Note that you will have to download the data, but you can be specific about your downloads based on variables

#sample random coordinates using terra
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
s <- spatSample(r, 10, "random", xy=TRUE, values=FALSE)

#set coordinate system to homosline crs as SoilGrids uses this
igh <- '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
sites <- vect(s, crs = igh) 
#You may include geom = (c"x", "y") in vect() if using a data.frame of coordinates

#download data using soil_world; giving an example of 1 variable
rpH <- soil_world(var="phh2o", depth=15, stat="mean", path="Rdata/") #I created an Rdata folder in my Rproj so I can organize the data better

#extract the SoilGrids data based on the coordinate
values_pH <- terra::extract(rpH, sites)
values_pH #check if any NA values? If so, may need to add na.rm = TRUE in extract()
sites <- cbind.data.frame(crds(sites), values_pH) #binding together data and the original coordinates.

This took me a while to figure out, but I hope it's helpful!

  • Thank you for your time. Really helpful! Mar 7, 2023 at 12:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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