alternative to depths<- for constructing a SoilProfileCollection

[brownag]

Based on the perennial challenges with teaching the depths<- SPC construction syntax, I have been thinking about providing an alternative.

We have had the SoilProfileCollection() function for several years. It is a slightly-safe (but not safe enough to use in general) low-level wrapper around the new("SoilProfileCollection") constructor that is used in one place internally.

If we define an S3 generic "SoilProfileCollection" function we could support methods for both character (current method) and data.frame (new method), and possibly other object types as input to create a S4 SoilProfileCollection object

I have prototyped this here: https://github.com/brownag/aqp/tree/spc-constructor1

With this setup we can keep the character method (which has never been recommended for public use but is exported), and add a data.frame method that doesnt rely on being on LHS of assignment operator or have formula/RHS argument syntax. Internally the prototype implementation uses depths<-. I am not recommending we deprecate that, but perhaps supersede it (suggest that it not be used in new code)

Also, I think we could easily make the character method internal, it is only used in a couple places, and as far as I am aware it has not been legitimately used in the wild. The new prototype data.frame method should be safer and also fully backward compatible provided that the arguments are explicitly specified (not relying on positional)

[dylanbeaudette]

Sure, have an example of the new syntax handy?

[brownag]

Here is what it currently looks like:

library(aqp)
#> This is aqp 2.3.1

data(sp4)

spc <- SoilProfileCollection(sp4, idcol = "id", depthcols = c("top", "bottom"))

spc
#> SoilProfileCollection with 10 profiles and 30 horizons
#> profile ID: id  |  horizon ID: hzID 
#> Depth range: 16 - 49 cm
#> 
#> ----- Horizons (6 / 30 rows  |  10 / 14 columns) -----
#>      id hzID top bottom name   K   Mg  Ca CEC_7 ex_Ca_to_Mg
#>  colusa    1   0      3    A 0.3 25.7 9.0  23.0        0.35
#>  colusa    2   3      8  ABt 0.2 23.7 5.6  21.4        0.23
#>  colusa    3   8     30  Bt1 0.1 23.2 1.9  23.7        0.08
#>  colusa    4  30     42  Bt2 0.1 44.3 0.3  43.0        0.01
#>   glenn    5   0      9    A 0.2 21.9 4.4  18.8        0.20
#>   glenn    6   9     34   Bt 0.3 18.9 4.5  27.5        0.20
#> [... more horizons ...]
#> 
#> ----- Sites (6 / 10 rows  |  1 / 1 columns) -----
#>         id
#>     colusa
#>      glenn
#>      kings
#>   mariposa
#>  mendocino
#>       napa
#> [... more sites ...]
#> 
#> Spatial Data:
#> [EMPTY]

# infers id if if there is one and only one shared column between horizons and site
# also can initialize spatial data or other data.frame slots in same call
spc2 <- SoilProfileCollection(
  sp4,
  site = data.frame(
    id = unique(sp4$id),
    x = 1:10,
    y = 1:10
  ),
  coordinates = c("x", "y"),
  crs = "EPSG:4326"
)

spc2
#> SoilProfileCollection with 10 profiles and 30 horizons
#> profile ID: id  |  horizon ID: hzID 
#> Depth range: 16 - 49 cm
#> 
#> ----- Horizons (6 / 30 rows  |  10 / 14 columns) -----
#>      id hzID top bottom name   K   Mg  Ca CEC_7 ex_Ca_to_Mg
#>  colusa    1   0      3    A 0.3 25.7 9.0  23.0        0.35
#>  colusa    2   3      8  ABt 0.2 23.7 5.6  21.4        0.23
#>  colusa    3   8     30  Bt1 0.1 23.2 1.9  23.7        0.08
#>  colusa    4  30     42  Bt2 0.1 44.3 0.3  43.0        0.01
#>   glenn    5   0      9    A 0.2 21.9 4.4  18.8        0.20
#>   glenn    6   9     34   Bt 0.3 18.9 4.5  27.5        0.20
#> [... more horizons ...]
#> 
#> ----- Sites (6 / 10 rows  |  3 / 3 columns) -----
#>         id x y
#>     colusa 1 1
#>      glenn 2 2
#>      kings 3 3
#>   mariposa 8 8
#>  mendocino 9 9
#>       napa 4 4
#> [... more sites ...]
#> 
#> Spatial Data:
#>   CRS: EPSG:4326; x: 1 to 10; y: 1 to 10


as(spc2, 'sf')
#> only site data are extracted
#> Simple feature collection with 10 features and 3 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 10 ymax: 10
#> Geodetic CRS:  WGS 84
#>                id  x  y      geometry
#> 1          colusa  1  1   POINT (1 1)
#> 2           glenn  2  2   POINT (2 2)
#> 3           kings  3  3   POINT (3 3)
#> 4        mariposa  8  8   POINT (8 8)
#> 5       mendocino  9  9   POINT (9 9)
#> 6            napa  4  4   POINT (4 4)
#> 7      san benito  5  5   POINT (5 5)
#> 8          shasta  6  6   POINT (6 6)
#> 9  shasta-trinity 10 10 POINT (10 10)
#> 10         tehama  7  7   POINT (7 7)

[dylanbeaudette]

Nice. I like this. The depths() <- constructor was cute, but not great for programmatic construction. Thanks for the addition. Merge when you are happy with it.

[brownag]

OK, now that I know you are amenable to it I will spend some time polishing it up and making sure I've handled the edge cases.