Commit 499b5537 authored by Batjes, Niels's avatar Batjes, Niels
Browse files

Upload New File

parent 49f96d7b
---
title: "Accessing WoSIS from R -- 'Latest' Version"
author: "D G Rossiter"
date: "`r format(Sys.Date(), '%d-%B-%Y')`"
output:
word_document:
toc: yes
html_document:
fig_height: 4
fig_width: 6
number_section: yes
theme: spacelab
toc: yes
toc_float: yes
bibliography: wosis.bib
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warnings = FALSE, purl = FALSE)
```
# Introduction
This document shows how to access WoSIS "Latest" data from R. For access to WoSIS "Snapshot" data from R, see `WoSIS_Snapshot_with_R.Rmd` at https://git.wur.nl/Batje001/wosis/-/tree/master/R_scripts.
The "Latest" dynamic dataset contains the most recent version of standardised soil data served from WoSIS. Being dynamic, the dataset will grow once new point data are standardised, additional soil properties are considered, and/or when possible amendments are required.
For an overview of WoSIS, see https://www.isric.org/explore/wosis. This links to https://www.isric.org/explore/wosis/accessing-wosis-derived-datasets which explains the difference between snapshot and dynamic datasets, and how to access them.
The [Procedures Manual](http://dx.doi.org/10.17027/isric-wdc-2020-01) describes how the WoSIS database is built.
# Packages
If you do not have these on your system, install with `install.packages(..., dependencies=TRUE)` or via the R Studio package manager.
```{r load.package}
library(rgdal) # GDAL access from R
library(gdalUtils) # wrappers for GDAL utility programs that could be
# called from the command line
library(sp) # spatial data types, ASDAR structures
library(sf) # spatial data types -- Simple Features
library(ggplot2) # gpplot graphics
```
GDAL is used for spatial data import/export, coordinate systems etc.
Check for a valid GDAL installation. If you have installed a binary version of `rgdal`, GDAL should be valid.
```{r set.gdal, eval=FALSE, purl=FALSE}
gdal_setInstallation(rescan=TRUE, verbose=TRUE)
valid_install <- !is.null(getOption("gdalUtils_gdalPath"))
if (valid_install)
print("Valid GDAL found") else { stop("No valid GDAL") }
```
If you do not have a valid GDAL installation, see https://gdal.org and the package DESCRIPTION for `rgdal`.
Check the drivers supported by your installation, and their capabilities.
```{r}
drivers <- ogrDrivers()
# print(drivers)
ix <- grep("GPKG", drivers$name, fixed=TRUE)
drivers[ix,]
ix <- grep("ESRI", drivers$name, fixed=TRUE)
drivers[ix,]
```
We see that ESRI Shapefiles and [OGC GeoPackages](https://www.geopackage.org) are among the supported formats.
# WoSIS Web Feature Service
The "latest" (dynamic) version of WoSIS is provided via
[WFS (Web Feature Services)](http://www.opengeospatial.org/standards/wfs).
This is described in less technical terms in [Wikipedia](https://en.wikipedia.org/wiki/Web_Feature_Service).
WFS allows you to incorporate geographic data into your own GIS projects, unlike WMS (Web Map Service), which only displays geographic data within a GIS. Here we use R as the GIS, since it can handle both geographic and feature-space information.
## Specify WoSIS WFS and list available layers
Specify the web address of the "latest" version of WoSIS:
```{r specify.wfs.address}
wfs <- "WFS:http://maps.isric.org/mapserv?map=/map/wosis_latest.map"
```
List the layers in the WFS source with the `rgdal::ogrinfo` function. This can be somewhat slow, because it depends on the network and the remote server. This function lists information about an OGR-supported data source, including this WFS source. OGR stands for "OGR Simple Features Library", which is now part of GDAL, which stands for "Geospatial Data Abstraction Library", see https://gdal.org for details.
Note: The argument `q=TRUE` suppresses reporting of coordinate system, layer schema, extents, and feature count. All we want here is the layer list. The argument `so=TRUE` specifies "summary only", `ro=TRUE` means that the data is accessed in read-only mode.
```{r ogrinfo.wfs}
system.time(
layers.info <- ogrinfo(wfs, ro=TRUE, so=TRUE, q=TRUE)
)
cat(layers.info, sep = "\n")
```
There are `r length(layers.info)` layers. Most of the names refer to soil properties in a fairly obvious way: the database name `ms:wosis_latest:wosis_latest_` and then a property name, e.g., `bdfi33`. These names are explained
[in the on-line documentation](https://www.isric.org/explore/wosis/accessing-wosis-derived-datasets).
For example `bdfi33` is "Bulk density of the fine earth fraction, equilibrated at 33 kPa". units are $\textrm{mg} \; \textrm{kg}^{-1}$.
## Display site properties
The layer `"ms:wosis_latest:wosis_latest_profiles"` contains the site information.
The `rgdal::ogrinfo` function displays information about a dataset accessible via GDAL. Here we see the metadata for the site information. Note the `so` "summary only" option. This supresses listing of features, and shows only the summary information. We choose this because there are many profiles in the entire dataset.
```{r ogrinfo.site}
system.time(
profiles.info <-
ogrinfo(wfs, layer="ms:wosis_latest_profiles",
ro=TRUE, so=TRUE, q=TRUE)
)
cat(profiles.info, sep="\n")
```
This shows the metadata of the sites.
The `rgdal::ogrinfo` function also allows an optional limitation to a bounding box with the `spat` argument, as a four-element vector `(xmin, ymin, xmax, ymax)`. These are the coördinates in the layer's Coordinate Reference System (CRS), in this case geographic coördinates.
For example, all the profile (site) information from a $2^\circ \times 2^\circ$ tile in central Europe.
Note: the `q=FALSE` option lists the geometry and feature count.
```{r ogrinfo.eu.profiles.1}
system.time(
central.eu.profiles.info <-
ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
layer="ms:wosis_latest_profiles",
spat=c(6, 48, 8, 50))
)
head(central.eu.profiles.info, 8)
central.eu.profiles.info[41:60]
```
Show the number of features, the spatial extent, and the Coordinate Reference System (CRS):
```{r ogrinfo.eu.profiles.2}
ix.f <- grep("Feature Count", central.eu.profiles.info)
central.eu.profiles.info[ix.f]
ix.e <- grep("Extent", central.eu.profiles.info)
central.eu.profiles.info[ix.e]
ix.g <- grep("GEOGCRS", central.eu.profiles.info)
cat(paste(central.eu.profiles.info[ix.g:(ix.g+17)], collapse="\n"))
```
Here we see the number of profiles in this tile, the extent (a bit smaller than the bounding box we requested), and the CRS formatted as ["Well-Known Text"](http://docs.opengeospatial.org/is/12-063r5/12-063r5.html).
The CRS corresponds to [EPSG code](https://www.epsg.org) `4326`, i.e., geographic coördinates on the WGS84 datum.
Show the data field names and their data types. These are listed after the `Geometry Column` line of the metadata:
```{r show.fields}
ix.p <- grep("Geometry Column", central.eu.profiles.info)
n <- length(central.eu.profiles.info)
central.eu.profiles.info[ix.p+1:n]
```
These are also found in the [Procedures Manual](http://dx.doi.org/10.17027/isric-wdc-2020-01).
The `rgdal::ogrinfo` function also allows [SQL (Structured Query Language](https://en.wikipedia.org/wiki/SQL) [queries](https://en.wikipedia.org/wiki/SQL_syntax) to limit the extent of the search, by using the optional `where` argument. To use this we need to know the field names, which we saw in the previous output.
For example, the profiles from India:
```{r ogrinfo.india.profiles}
system.time(
india.profiles.info <-
ogrinfo(wfs, ro=TRUE, so=TRUE, q=FALSE,
layer="ms:wosis_latest_profiles",
where="country_name='India'")
)
```
Show the number of records and geographic extent:
```{r ogrinfo.india.profiles.2}
ix.f <- grep("Feature Count", india.profiles.info)
india.profiles.info[ix.f]
ix.e <- grep("Extent", india.profiles.info)
india.profiles.info[ix.e]
```
There are only `r india.profiles.info[ix.f]` profiles from India. Of course many more have been described, but despite the aim of WoSIS to be a complete world database, up till now it has proved impossible to develop a data-sharing arrangement with the [NBSSLUP](https://nbsslup.in/).
## Display layer properties
We need to know the layer's meaning, format, extent etc. before we read it into R.
Select the first property in the information list above. Again call `ogrinfo` but this time also specifying a layer within the WFS. Show only the summary, not the features, with the `so=TRUE` optional argument.
Note: the `q=FALSE` option lists the geometry and feature count.
```{r ogrinfo.properties, warnings}
system.time(
property.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
ro=TRUE, so=TRUE, q=FALSE)
)
cat(property.info, sep="\n")
ix.f <- grep("Feature Count", property.info)
property.info[ix.f]
```
This gives an explanation of ISRIC and WoSIS, and then the keywords showing this property.
We see this is "Bulk density of the fine earth fraction* equilibrated at 33 kPa". Refer to the procedures manual for how each property (here, bulk density) was determined. The footnote `*` gives the definition of this term for this dataset. See also the feature count -- this is how many layers have a bulk density value.
The `rgdal::ogrinfo` function allows [SQL (Structured Query Language](https://en.wikipedia.org/wiki/SQL) [queries](https://en.wikipedia.org/wiki/SQL_syntax) to limit the extent of the search, by using the optional `where` argument. To use this we need to know the field names, which we saw in the previous output.
For example, all the subsoil bulk densities from India:
```{r ogrinfo.bd.india}
system.time(
bd.india.info <- ogrinfo(wfs, layer="ms:wosis_latest_bdfi33",
where="country_name='India' AND upper_depth > 100",
ro=TRUE, so=TRUE, q=FALSE)
)
ix.f <- grep("Feature Count", bd.india.info)
bd.india.info[ix.f]
(n.records <- as.numeric(strsplit(bd.india.info[ix.f],
split=": ", fixed=TRUE)[[1]][2]))
```
Here there are only `r n.records` records
# Import WoSIS datasets to the client system
There seems to be no way to directly import from the WFS to an R workspace object, so there must first be an intermediate step: download the WFS layer in an appropriate format to a local directory, and then import as usual for GIS layers.
The `rgdal::ogr2ogr` function reads from one format on the server and writes to another in our client. The default output file format is an [ESRI Shapefile](https://en.wikipedia.org/wiki/Shapefile); other formats can be specified with the (optional) `f` argument.
The possible formats are listed [here](https://docs.geoserver.org/latest/en/user/services/wfs/outputformats.html). In this example we will use the default format,
i.e. the ESRI shapefile. See the next section for another format, CSV.
The `where` and `spat` arguments can also be used here. You will often want to limit the size of the object with one or both of these. The `ogr2ogr` function also allows an (optional) transformation to any EPSG-defined CRS with the `t_srs` argument. See `?ogr2ogr` for more options.
There are several choices of file format, here we explain three: [ESRI Shapefiles](#shape), [CSV files](#csv), and [OGC Geopackages](#gpkg).
We set up a directory on the local file system to receive the downloaded files:
```{r setup.local.dir}
wosis.dir.name <- "./wosis_latest"
if (!file.exists(wosis.dir.name)) dir.create(wosis.dir.name)
```
Note that **import via `ogr2ogr` can be quite slow**, because it depends on the network and the remote server, which may have a speed limitation to avoid overloads. Many of these downloads may take 15-20 minutes clock time, while only requiring less than a minute of local computer time.
## ESRI shapefiles {#shape}
### WoSIS profiles as ESRI Shapefiles
First, download the profile information for the whole world.
```{r download.profiles}
layer.name <- "wosis_latest_profiles"
system.time(
ogr2ogr(src=wfs,
dst=wosis.dir.name,
layer=layer.name,
f = "ESRI Shapefile",
overwrite=TRUE,
skipfailures=TRUE)
)
```
The number of profiles can be restricted with a `spat` spatial extent or a `where` SQL query.
For example, to download just the Indian profiles, into a subdirectory:
```{r download.profiles.india}
wosis.dir.name.india <- "./wosis_latest/india"
if (!file.exists(wosis.dir.name.india)) dir.create(wosis.dir.name.india)
layer.name <- "wosis_latest_profiles"
system.time(
ogr2ogr(src=wfs,
dst=wosis.dir.name.india,
layer=layer.name,
f = "ESRI Shapefile",
where="country_name='India'",
overwrite=TRUE,
skipfailures=TRUE)
)
```
### Reading imported profiles into R
Now read the downloaded shapefile of profiles into an R `sp` object.
For shapefiles, the directory and layer names must be specified separately as two arguments, `dsn` ("data set name") and `layer`. Notice how the server name of this layer which begins with `ms:` has been changed to `ms_` during import, due to restrictions on file names on the local file system.
```{r import.profiles}
layer.name <- "ms_wosis_latest_profiles"
profiles <- readOGR(dsn=wosis.dir.name, layer=layer.name,
stringsAsFactors = FALSE)
class(profiles)
dim(profiles)
names(profiles@data)
head(profiles@data, 1)
```
The current database has `r dim(profiles)[2]` profiles. Quite a resource!
The Indian profiles:
```{r import.profiles.india}
profiles.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
stringsAsFactors = FALSE)
dim(profiles.india)
```
Here there are only `r dim(profiles.india)[2]` profiles, as we saw with `ogrinfo`, above.
### WoSIS layers as ESRI Shapefiles
Second, download all records for the first property, i.e., bulk density, and save it a subdirectory of the current path, and with the default format. The directory must be first created if does not already exist.
We need to know the formats which are supplied by the provider; these are specified with the `f=` argument.
```{r download.bd}
layer.name <- "wosis_latest_bdfi33"
system.time(
ogr2ogr(src=wfs,
dst=wosis.dir.name,
layer=layer.name,
f = "ESRI Shapefile",
overwrite=TRUE,
skipfailures=TRUE)
)
```
Because the destination format is `"ESRI Shapefile"` many of the field names have to be shortened, as shown by warnings such as "Warning 6: Normalized/laundered field name: 'profile_layer_id' to 'profile_la'".
Note that `spat` is not applicable in this query because there is no spatial information in the attribute tables. However, it is possible to include an SQL query with `where` to limit the size of the download:
```{r download.bd.india}
system.time(
ogr2ogr(src=wfs,
dst=wosis.dir.name.india,
layer=layer.name,
f = "ESRI Shapefile",
where="country_name='India'",
overwrite=TRUE,
skipfailures=TRUE)
)
```
### Reading imported layer into R
Now read the downloaded shapefile into an R `sp` object.
For shapefiles, the directory and layer names must be specified separately as two arguments, `dsn` ("data set name") and `layer`. Notice how the server name of this layer which begins with `ms:` has been changed to `ms_` during import, due to restrictions on file names on the local file system.
```{r import.bd}
# here strings are just that, not to be interpreted as R factors
layer.name <- "ms_wosis_latest_bdfi33"
bd33 <- readOGR(dsn=wosis.dir.name, layer=layer.name,
stringsAsFactors = FALSE)
class(bd33)
names(bd33@data)
head(bd33@data, 1)
```
Each record has some identification:
* `gml_id` the attribute name + `profile_la` (see below)
* `profile_id` profile internal ID
* `profile_la` profile + layer internal ID
* `country_na` country name
* `upper_dept` upper limit of layer from soil surface (excluding litter layer). cm
* `lower_dept` lower limit of layer from soil surface (excluding litter layer), cm
* `layer_name` layer name as assigned during original profile description
* `litter` *** ???
Each attribute has several names, with the following extensions:
* `value` -- one or more values, in the format {1:value; 2:value...}, which are duplicate measurements
* `value_avg` -- the average of the values
* `method` -- text description of the analytical method
* `date` -- one or more values, in the format {1:yyyy-mm-dd; 2:yyyy-mm-dd...},
which are the dates each of th duplicate measurements was added to the database (not the original measurement date, nor the field sampling date)
* `dataset_id` -- text code of original database
* `profile_code` -- text code of profile from original database
* `licence` -- text string of the Creative Commons ^[https://creativecommons.org/licenses/] license for this value, e.g. `CC-BY-NC`
So for example the attribute `bdfi33` has the following fields, shortened when input to a shapefile to the first ten characters; if these would be duplicated the name is further manipulated:
* `bdfi33_value` $\to$ `bdfi33_val`
* `bdfi33_value_avg` $\to$ `bdfi33_v_1`
* `bdfi33_method` $\to$ `bdfi33_met`
* `bdfi33_date` $\to$ `bdfi33_dat`
* `bdfi33_dataset_id` $\to$ `bdfi33_d_1`
* `bdfi33_profile_code` $\to$ `bdfi33_pro`
* `bdfi33_licence` $\to$ `bdfi33_lic`
The shapefile has been imported as a `SpatialPointsDataFrame`, a class within the `sp` package, with the correct CRS, as we saw from the `ogrinfo`.
Examine the format of the attribute, this is in field `*_val`:
```{r head.bd}
head(bd33@data$bdfi33_val)
```
The format is `{seq:val[,seq:val]}` where the `seq` is an integer on `[1...]` indicating which measurement number -- note that there can be more than one measurement per property, e.g., repeated lab. measurements, and `val` is the numeric value.
But the average value has its own field, so if we only want the average, it is prepared for us. We see an example here, from six rows chosen to show several profiles with their layers:
```{r example.bd}
bd33@data[75:80, c("profile_id","upper_dept","lower_dept","bdfi33_val","bdfi33_v_1")]
```
Here is the India example. We make a histogram of the representative values.
```{r import.bd.india}
layer.name <- "ms_wosis_latest_bdfi33"
bd33.india <- readOGR(dsn=wosis.dir.name.india, layer=layer.name,
stringsAsFactors = FALSE)
class(bd33.india)
dim(bd33.india)
(profile.id.india <- unique(bd33.india@data$profile_id))
names(bd33.india@data)
hist(bd33.india@data$bdfi33_v_1, main="Bulk density, Layers in India",
xlab="g cm-3")
```
Here there are `r dim(bd33.india)[1]` records in `length(unique(bd33.india@data$profile_id))` profiles. These profiles can be found in the profiles database:
```{r}
```
## CSV files {#csv}
Another output format for `ogr2ogr` is the CSV 'comma-separated values' plain-text file. These typically have one header row giving the name of each column (field), and then one row (case, tuple) per observation
### WoSIS profiles as CSV files
For example, read the profile information for the $2^\circ \times 2^\circ$ tile in central Europe.
(Note: `overwrite=TRUE` does not seem to work for `.csv` files, so any previous must be manually deleted).
```{r download.profiles.csv}
wosis.dir.name.ceu <- "./wosis_latest/central_europe"
if (!file.exists(wosis.dir.name.ceu)) dir.create(wosis.dir.name.ceu)
src.layer.name <- "ms:wosis_latest_profiles"
dst.layer.name <- "wosis_latest_profiles_ceu"
(dst.target.name <- paste0(wosis.dir.name.ceu,"/",dst.layer.name,".csv"))
if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
ogr2ogr(src=wfs,
dst=dst.target.name,
layer=src.layer.name,
f="CSV",
spat=c(6, 48, 8, 50),
overwrite=TRUE)
round(file.info(dst.target.name)$size/1024,1)
```
This file is about `r round(file.info(dst.target.name)$size/1024)` Kb.
###Read imported profiles into R
The `read.csv` function reads from a CSV file into an R `data.frame`.
Read the profiles (sites) from central Europe:
```{r import.profiles.csv}
layer.name <- "wosis_latest_profiles_ceu"
system.time(
profiles.ceu <- read.csv(paste0(wosis.dir.name.ceu, "/",layer.name,".csv"),
stringsAsFactors = FALSE)
)
names(profiles.ceu)
```
### WoSIS layers as CSV files
Get the the bulk density for the layers in these profiles. Note that this query can not be limited by coordinates, since they are not included in this table. So we get all the layers and then limit to the profiles of interest. This is a very large file, about 97Mb.
```{r download.bd.csv}
src.layer.name <- "wosis_latest_bdfi33"
dst.layer.name <- "wosis_latest_bdfi33"
(dst.target.name <- paste0(wosis.dir.name,"/",dst.layer.name,".csv"))
if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
ogr2ogr(src=wfs,
dst=dst.target.name,
layer=src.layer.name,
f="CSV",
overwrite=TRUE)
file.info(dst.target.name)$size
```
### Read imported layers into R
The bulk density per-layer.
```{r import.bd.csv}
layer.name <- "wosis_latest_bdfi33"
system.time(
bd33.pts <- read.csv(file=paste0(wosis.dir.name,"/",layer.name,".csv"),
stringsAsFactors = FALSE)
)
dim(bd33.pts)
names(bd33.pts)
```
Notice that here we have the complete field names,
not truncated as they were in the shapefiles.
These can now be processed as above, i.e., the shapefiles.
## OGC Geopackage {#gpkg}
WoSIS WFS is also available as a OGC GeoPackage, using `f="GPKG` argument.
(Note: `overwrite=TRUE` does not seem to work for `.gpkg` files, so any previous must be manually deleted).
For example, the profiles as a GeoPackage:
```{r download.profile.gpkg}
layer.name <- "wosis_latest_profiles"
(dst.target.name <- paste0(wosis.dir.name,"/", layer.name, ".gpkg"))
if (file.exists(dst.target.name)) { file.remove(dst.target.name) }
system.time(
ogr2ogr(src=wfs,
dst=dst.target.name,
layer=layer.name,
f = "GPKG",
overwrite=TRUE,
skipfailures=TRUE)
)
```
This Geopackage is about 52 Mb.
### Reading imported profiles into R
The Geopackage can be read with the `st_read` function of the `sf` "Simple Features" package.
```{r read.profile.gpkg}
profiles.gpkg <- st_read(dst.target.name)
class(profiles.gpkg)
dim(profiles.gpkg)
names(profiles.gpkg)
```
Many of these names have no contents at the whole-profile (site) level.
Which countries have the most profiles?
```{r}
head(sort(table(profiles.gpkg$country_name), decreasing=TRUE))
```
These have georeference and so can be mapped. Here is an example of the profiles from Brazil, showing their source datasets:
```{r map.profiles.gpkg, fig.width=12, fig.height=12}
ggplot(data=profiles.gpkg[(profiles.gpkg$country_name=="Brazil"), ]) +
aes(col=dataset_id) +
geom_sf(shape=21, size=0.8, fill="black")
```
# Spatial objects
Once the WoSIS profiles have been imported to R, either from ESRI Shapefiles or CSV files, they can be converted to a spatial object in the `sp` package, by specifying the coördinates and the Coordinate Reference System (CRS).
```{r coordinates.profiles}
coordinates(profiles.ceu) <- c("longitude", "latitude")
proj4string(profiles.ceu) <- CRS("+init=epsg:4326")
class(profiles.ceu)
```
These can also be represented as Simple Features from the `sf` "Simple Features" package:
```{r coordinates.profiles.sf}
profiles.ceu.sf <- st_as_sf(profiles.ceu)
class(profiles.ceu.sf)
```
Review some site information, e.g., the WRB Reference Soil Groups:
```{r profiles.wrb}
table(profiles.ceu@data$cwrb_reference_soil_group)
table(profiles.ceu.sf$cwrb_reference_soil_group)
```
Note that most of these profiles do not have a WRB classification.
And a map:
```{r map.profiles.wrb}
ggplot(data=profiles.ceu.sf) +
aes(col=cwrb_reference_soil_group) +
geom_sf()
```
Note that most of these profiles do not have a WRB classification.
# Working with WoSIS as a `SoilProfileCollection`
The `aqp` "Algorithms for Quantitive Pedology" package [@Beaudette_2013] defines data structures and functions specific to soil profile data, i.e., with site and linked layer information.
Load the package:
```{r load.aqp}
library(aqp) # Algorithms for Quantitative Pedology
```
Convert the bulk density structure to a `SoilProfileCollection`, a data type defined in `aqp`. This data type has separate structures for the site (profile) and horizons.
Initialize with the horizons from the data frame created from the WoSIS layer.
The `aqp::depth` function initializes the SoilProfileCollection object. The formula has the field name of the profile on the left, and the the field names of the horizon boundaries on the right. These fields are in the WoSIS layer.
```{r bd.aqp}
ds.aqp <- bd33@data
depths(ds.aqp) <- profile_id ~ upper_dept + lower_dept
is(ds.aqp)
slotNames(ds.aqp)
str(ds.aqp@site)
str(ds.aqp@horizons)
head(ds.aqp@site)
head(ds.aqp@horizons[c(2,5,6,7,9)],12)
```
Note how the horizons have been grouped into sites, in the `@site` slot, and the per-horizon (by depth) values are in the `@horizons` slot. Here we have `r dim(ds.aqp@horizons)[1]` horizons in `r dim(ds.aqp@site)[1]` profiles.
Now this `SoilProfileCollection` can be used for many `aqp` functions. For example, here is the depth distribution of average bulk density of the components for the first 24 listed profiles, labelled by genetic horizon:
```{r plot.bd.spc,fig.width=12, fig.height=8}
plotSPC(ds.aqp[1:24,], name="layer_name", color='bdfi33_v_1')
```
A few layers in this set of profiles are missing bulk density.
# References
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment