Soil Grids access

[kristoffer wild]

I was wondering if anyone has any easy workarounds for querying the soilgrids.org database for soil hydraulic properties when running the argument 'soilgrids =1' in any of the microclimate functions (e.g., micro_NCEP, micro_era5). At the moment, I query SoilGrids via the WCS service (maps.isric.org) for bdod, clay, silt, and sand at standard depths, run NicheMapR::pedotransfer to derive PE/BB/BD/KS, then run the microclimate model with soilgrids = 0 and pass those soil parameters directly. This takes some time. Do you know of a better workaround or a more direct link? It's been about 3 weeks since the argument/link has been down.

Best,

Kris 

[NicheMapR]

Hi Kris,

Perhaps not an easy workaround, but you could do this in Julia (currently implemented a fork from the main branch https://github.com/EcoJulia/RasterDataSources.jl soon):

(@v1.12) pkg> add https://github.com/mrke/RasterDataSources.jl#soil-grids

using RasterDataSources
using RasterDataSources: SOILGRIDS_QUANTILES

ENV["RASTERDATASOURCES_PATH"] = "where you store your rasters" 

for layer in (:sand, :silt, :clay, :bdod)
    for depth in depths(SoilGrids, layer)
        for quantile in SOILGRIDS_QUANTILES
            getraster(SoilGrids, layer; depth=depth, quantile=quantile)
        end
    end
end

To access the data after it downloads:

using Rasters, ArchGDAL, RasterDataSources
using Extents, Proj

all_depths = RasterDataSources.depths(SoilGrids)  # ("0-5cm","5-15cm",...)

# Download/patch VRTs for all depths (returns Vector of NamedTuples)
results = getraster(SoilGrids, (:sand, :silt, :clay); depth=collect(all_depths), quantile="mean")

# Set up coordinate transform once
stack0 = RasterStack(results[1], lazy=true)
trans = Proj.Transformation("EPSG:4326", convert(String, crs(stack0)); always_xy=true)
x, y = trans(150.2093, -33.8688)  # location from lon/lat to in IGH meters

# Extract point value at each depth
rows = map(zip(all_depths, results)) do (depth, paths)
    stack = RasterStack(paths, lazy=true)
    pt = stack[Near(x), Near(y)]
    (depth=depth, clay=pt[:clay], sand=pt[:sand], silt=pt[:silt])
end

df = DataFrame(rows)

and if you want to extract a plot a subset do something like:

result_clay = getraster(SoilGrids, :clay; depth="0-5cm", quantile="mean")
clay_raster = Raster(result_clay, lazy=true)

trans = Proj.Transformation("EPSG:4326", convert(String, crs(clay_raster)); always_xy=true)

lon, lat = 144.9631, -33.8136  # Melbourne
x_c, y_c = trans(lon, lat)
x_c, y_c = lon, lat
half = 25000.0  # 25 km → 5×5 km box
half = 0.25
# Crop in native IGH CRS

clay_region = read(clay_raster[
    X(Rasters.Between(x_c - half, x_c + half)),
    Y(Rasters.Between(y_c - half, y_c + half))
])

# Raw values are in g/kg; replace missing with NaN, divide by 10 for % clay (g/100g)
clay_pct = replace_missing(clay_region, NaN32) ./ 10f0

fig, ax, plt = heatmap(clay_pct;
    colormap=:YlOrBr,
    axis=(; title="Clay 0–5 cm, 5×5 km", aspect=DataAspect()))
Colorbar(fig[1, 2], plt; label="Clay (%)")
fig

There's also functionality to download the Soil and Landscape Grid of Australia in the same way. The download code for that is:

using RasterDataSources
using RasterDataSources: SLGA_COMPONENTS

for layer in layers(SLGA)
    for depth in depths(SLGA, layer)
        for component in SLGA_COMPONENTS
            getraster(SLGA, layer; depth=depth, component=component)
        end
    end
end

Then run the same code above but replace SoilGrids with SLGA, just use decimal degrees longitude/latitude in the query for that and delete:  quantile = "mean". Note that SoilGrids just downloads VRT (virtual raster) so when you request data it takes a little while but for the SLGA it downloads the whole files which are big (~20 gigs for all depths for sand, silt, clay and bulk density).

All the best,
Mike