Getting spatial data on soil profiles in R

Soil structure and health is critical to water availability, nutrient cycling, and plant productivity. By extension, soils will have strong associations with invertebrate community assemblage, diversity, and abundance and is thus a invertebrate science.

Modern spatial datasets on soil are available to help us intergrate variation in soil characteristics in the prediction of invertebrate processes. In Austalia, one of the newer soil data sets is Soil and Landscape grid of Australia. The data set is available as gridded rasters with a cell resolution of 3" (roughly 90m) and contains a range of variables including Bulk Density (Whole Earth), Organic Carbon, Clay, Silt, Sand, pH, and more. In addition, the data product includes information of landscape level variables such as Aspect, Relief 1000m Radius, Relief 300m Radius, Topographic Wetness Index, Topographic Position Index, and more. See the home page for more details.

Here we see how the data set can be quickly accessed and visualised inside R through the slga package available on github.

First we install the package.


Then load the package.


We can see what soil products are available.

Type Product Short_Name Code
Soil National_3D NAT_3D ACLEP_AU_TRN_N
Soil National NAT ACLEP_AU_NAT_C
Soil South_Australia SA ACLEP_AU_SAT_D
Soil Tasmania TAS ACLEP_AU_TAS_N
Soil Western_Australia WA ACLEP_AU_WAT_D
Landscape Prescott Index PSIND SRTM_attributes_3s_ACLEP_AU
Landscape Net Radiation [January] NRJAN SRTM_attributes_3s_ACLEP_AU
Landscape Net Radiation [July] NRJUL SRTM_attributes_3s_ACLEP_AU
Landscape Total Shortwave Sloping Surf [January] TSJAN SRTM_attributes_3s_ACLEP_AU
Landscape Total Shortwave Sloping Surf [July] TSJUL SRTM_attributes_3s_ACLEP_AU
Landscape Slope [percent] SLPPC SRTM_attributes_3s_ACLEP_AU
Landscape Slope [percent] Median 300m Radius SLMPC SRTM_attributes_3s_ACLEP_AU
Landscape Slope Relief Class RELCL SRTM_attributes_3s_ACLEP_AU
Landscape Aspect ASPCT SRTM_attributes_3s_ACLEP_AU
Landscape Relief [1000m radius] REL1K SRTM_attributes_3s_ACLEP_AU
Landscape Relief [300m radius] REL3C SRTM_attributes_3s_ACLEP_AU
Landscape Topographic Wetness Index TWIND SRTM_attributes_3s_ACLEP_AU
Landscape TPI Mask TPMSK SRTM_attributes_3s_ACLEP_AU
Landscape SRTM_TopographicPositionIndex TPIND SRTM_attributes_3s_ACLEP_AU
Landscape Contributing Area [partial] CAPRT SRTM_attributes_3s_ACLEP_AU
Landscape MrVBF MRVBF SRTM_attributes_3s_ACLEP_AU
Landscape Plan Curvature PLNCV SRTM_attributes_3s_ACLEP_AU
Landscape Profile Curvature PRFCV SRTM_attributes_3s_ACLEP_AU
Name Code Units Transformation
Available Water Capacity AWC % None
Bulk Density (Fine Earth) BDF g/cm None
Bulk Density (Whole Earth) BDW g/cm None
Cation Exchange Capacity CEC meq/100g Log
Cation Exchange Capacity (Effective) ECE meq/100g Log
Clay CLY % None
Coarse Fragments CFG % None
Depth of Regolith DER Meters None
Depth of Soil DES Meters None
Electrical Conductivity ECD dS/m Log
Organic Carbon SOC % Log
pH CaCl2 PHC pH Units None
pH Water PHW pH Units None
Sand SND % None
Silt SLT % None
Total Nitrogen NTO % Log
Total Phosphorus PTO % Log

Now let’s get some data on my favourite place in Victoria - Gariwerd (Grampians National Park). Define the bounding box of the area of interest aoi using longitudes and latitudes in the form c(xmin, ymin, xmax, ymax). Using the national dataset product NAT we get the estimated value VAL of the percentage of sand SND at a depth of 0-5 cm depth = 1.

aoi <- c(142.138939, -37.645513, 142.694520, -36.843094)
sand_percent <- get_soils_data(product = 'NAT', attribute = 'SND',
                        component = 'VAL', depth = 1,
                        aoi = aoi, write_out = FALSE)

We can also extract information for a single point of interest (poi).

poi <- c(142.798591, -37.214808)
organic_carbon_percent <- get_soils_point(product = 'NAT', attribute = 'SOC',
                        component = 'ALL', depth = 1,
                        poi = poi)
##   NAT_SOC_VAL_000_005 NAT_SOC_CLO_000_005 NAT_SOC_CHI_000_005
## 1            5.296265            1.480427            22.69777

With one line of code we can see the data interactively usign the mapview package, but for more customisability I highly recommend leaflet for interactive maps (mapview is a wrapper for this) and tmap for static maps. Note that the number of pixels is slightly reduced during plotting with mapview.