calc_AddRequestedSoilLayers not always working

[Zachary-Kramer]

Referencing https://github.com/DrylandEcology/rSFSW2/blob/355ca39c06ca6e1545c962274664128dd9298557/R/PriorCalculations.R

The loop that adds layers to the input data (lines 97:104) seems to be wrong. It's hard to explain but I will try with an example:

I have some soil data with the depths 5, 33 and clay values 0.10, 0.31. I am trying to interpolate this to 5, 10, 20, 30, 40, .... I am expecting the function to add the values 10, 20, 30 for a total of 5, 10, 20, 30, 33.

Before line 100, our clay values are stored like this:

Clay_L1 Clay_L2
0.10      0.31

After adding layer depth 10 (lines 100:102), we have the following clay values:

         [1]   [2]
Clay_L1: 0.10 0.10
Clay_L2: 0.31 0.31

After adding depth 20:

         [1]  [2]   [3]
Clay_L1: 0.10 0.10 0.10
Clay_L2: 0.31 0.31 0.31

This pattern continues; the rows are now permanently Clay_L1 and L2 and the new layers are being added as columns. I don't understand this. Surely this is not intended? This is also on master. This causes an indexing error on line 150 because we are one column short.

I don't have the time to fully investigate this, but it seems like add_layer_to_soil is at fault. It converts the named vector to a matrix, causing the column names to become row names:

         [1]
Clay_L1: 0.10
Clay_L2: 0.31

Now we only have 1 column or "layer". Depth 10 is trying to get inserted at index 1 (il = 1), but ncol == ilso it is now appended and becomes column 2, per the function's logic. I tried changing the conversion to a matrix to matrix(data = x, nrow = 1, ncol=length(x), dimnames = list(NULL, names(x))), which preserves the number of columns/rows correctly, but this caused other issues.

Could you please look into this and let me know what aspect I am not understanding?

High priority because I am leaving USGS in the coming weeks and this is the main roadblock to finishing SSURGO interpolation. Assigned to @dschlaep because I believe this code was originally written by him.

---

This should be replicable using the following data:

Soil layers:

SoilDepth_cm	depth_L1	depth_L2
33	5	33

Soils:

Matricd_L1	GravelContent_L1	EvapCoeff_L1	Grass_TranspCoeff_L1	Shrub_TranspCoeff_L1	Tree_TranspCoeff_L1	Forb_TranspCoeff_L1	TranspRegion_L1	Sand_L1	Clay_L1	Imperm_L1	SoilTemp_L1	Matricd_L2	GravelContent_L2	EvapCoeff_L2	Grass_TranspCoeff_L2	Shrub_TranspCoeff_L2	Tree_TranspCoeff_L2	Forb_TranspCoeff_L2	TranspRegion_L2	Sand_L2	Clay_L2	Imperm_L2	SoilTemp_L2
1	1	1	0	0	0	0	0	1	1	0	0	1	1	1	0	0	0	0	0	1	1	0	0
1.35	0.1625	0.8111	NA	NA	NA	NA	NA	0.669	0.1	NA	NA	1.33	0.095	0.1889	NA	NA	NA	NA	NA	0.354	0.31	NA	NA

[dschlaep]

My hunch is that this a problem with [ dropping dimensions ( see ?[) because there is one row as the example has one site only. Add , drop = FALSE to appropriate sub-setting calls.

[Zachary-Kramer]

Adding a drop=FALSE to a specific point in calc_AddRequestedSoilLayers did accomplish the same as using matrix(data = x, nrow = 1, ncol=length(x), dimnames = list(NULL, names(x))), but obviously in a much cleaner/better manner.

There were still issues afterwards (as I've gotten to this point before), but they were subsequently solved with another well-placed drop=FALSE in add_layer_to_soil. It appears I was originally only one line off from getting it to work...

The original interpolation now works as expected and I can finish up SSURGO's variation of it. It pays to have another pair of eyes. Thanks.

[dschlaep]

Good. That the function [ has as default drop = TRUE is convenient for interactive data wrangling, but a bane for programming with R...