WebDAV: R - Get XY info

Get SoilGrids layers values for a list of points

Load libraries

library(terra)
library(sf)

Set variables of interest

We define the variables for the soil property and layer of interest. See here for the naming conventions

voi = "soc" # variable of interest
depth = "0-5cm"
layer = "mean"

voi_layer = paste(voi,depth,layer, sep="_") # layer of interest 

We set other variables necessary for the WCS call for all kinds of requests

webdav_path = '/vsicurl?max_retry=3&retry_delay=1&list_dir=no&url=https://files.isric.org/soilgrids/latest/data/'

Read in point data and re-project to Homolosine

data=read.csv("my_points_data.csv")
data <- data.frame(
  id = 1:5,  # Unique IDs
  longitude = c(34.78, 35.12, 34.95, 35.05, 34.88),
  latitude = c(-1.29, -1.31, -1.27, -1.30, -1.28),
  value = runif(5, 10, 50)  # Some random values
)
spdata=st_as_sf(data,coords = c("longitude", "latitude"), crs = 4326)

igh='+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs'
spdata_igh=st_transform(spdata, igh)

data_igh=data.frame(st_coordinates(spdata_igh),id=spdata_igh$id)

Use the extract function of terra to get the values from the pixels

fun_pixel_values=function(data,VOI,VOI_LYR){
        extract(
            x=rast(paste0(webdav_path,"/",VOI,"/", VOI_LYR,".vrt")),
            y=data[, c('X', 'Y')],
            xy = FALSE,
            ID = FALSE)
}

value_pixels = fun_pixel_values(data_igh,voi,voi_layer)

The value_pixels vector can be binded with data_igh paying attention to order. The extract function of terra can also return more information, including the coordinates themselves for a safer workflow (in that case set xy = TRUE).

Last updated on 2025-12-12