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 gdallocationinfo to get the values from the pixels

fun_pixel_values=function(rowPX,data,VOI,VOI_LYR){
    as.numeric(
        gdallocationinfo(
            srcfile=paste0(webdav_path,"/",VOI,"/", VOI_LYR,".vrt"),
            x=data[rowPX,"X"],
            y=data[rowPX,"Y"],
            geoloc=TRUE,
            valonly=TRUE))
}

value_pixels=unlist(lapply(1:3,function(x){fun_pixel_values(x,data_igh,voi,voi_layer)}))

The value_pixels vector can be binded with data_igh paying attention to order. gdallocationinfo can also return more information, including the coordinates themselves for a safer workflow.

Last updated on 2025-10-16