SoilUtils.java [src/soils/utils] Revision: default Date:
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package soils.utils;
import csip.Config;
import java.text.NumberFormat;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
import java.util.function.Function;
import soils.Component;
import soils.Fragments;
import soils.Horizon;
import soils.MapUnit;
import soils.TextureGroup;
/**
*
* @author <a href="mailto:shaun.case@colostate.edu">Shaun Case</a>
*/
public class SoilUtils {
static final List<String> PALOUSE_AREAS = Arrays.asList("ID607", "ID610", "OR021", "OR049", "OR055", "OR625",
"OR667", "OR670", "OR673", "WA001", "WA021", "WA025",
"WA043", "WA063", "WA071", "WA075", "WA603", "WA605",
"WA613", "WA617", "WA623", "WA639", "WA676", "WA677");
static final double[][] LIGHTLE_WEESIES_SLOPE_LENGTH = {
{0, .75, 100},
{.75, 1.5, 200},
{1.5, 2.5, 300},
{2.5, 3.5, 200},
{3.5, 4.5, 180},
{4.5, 5.5, 160},
{5.5, 6.5, 150},
{6.5, 7.5, 140},
{7.5, 8.5, 130},
{8.5, 9.5, 125},
{9.5, 10.5, 120},
{10.5, 11.5, 110},
{11.5, 12.5, 100},
{12.5, 13.5, 90},
{13.5, 14.5, 80},
{14.5, 15.5, 70},
{15.5, 17.5, 60},
{17.5, -1, 50}
};
static final double[][] PALOUSE_SLOPE_LENGTH = {
{0, 5.5, 350},
{5.5, 10.5, 275},
{10.5, 15.5, 225},
{15.5, 20.5, 175},
{20.5, 25.5, 150},
{25.5, 35.5, 125},
{35.5, -1, 100}
};
static boolean force = Config.getBoolean("force_calculations", false);
static boolean sdmksatOnly = Config.getBoolean("sdm.ksat.only", false);
public static boolean isPalouse(String areasymbol) {
return PALOUSE_AREAS.contains(areasymbol);
}
public static double calcLightleWeesieSlopeSlopeLength(double slope_r, boolean isPalouse) {
double[][] slopeLength = isPalouse ? PALOUSE_SLOPE_LENGTH : LIGHTLE_WEESIES_SLOPE_LENGTH;
// Calculate lambda
for (int i = 0; i < slopeLength.length; i++) {
// Area at the end of the list of possible slope lengths?
if (slopeLength[i][1] != -1) {
if ((slope_r >= slopeLength[i][0]) && (slope_r < slopeLength[i][1])) {
// Found a slope length to use, set and leave.
return slopeLength[i][2];
}
// Else, we have no more data in the array, assign the last known value if possible.
} else if (slope_r >= slopeLength[i][0]) {
return slopeLength[i][2];
}
}
return Double.NaN;
}
public static double calcLsFactor(double slope_r, double slope_length) {
double nn = 0.2;
if ((1.0 <= slope_r) && (3 > slope_r)) {
nn = 0.3;
} else {
if ((3 <= slope_r) && (5 > slope_r)) {
nn = 0.4;
} else {
if (5 <= slope_r) {
nn = 0.5;
}
}
}
return (0.065 + 0.0456 * slope_r + 0.006541 * Math.pow(slope_r, 2)) * (Math.pow((slope_length / 72.5), nn));
}
public static double averageMissingValue(double h, double l, double r) {
return averageMissingValue(h, l, r, Double.NaN);
}
public static double averageMissingValue(double h, double l, double r, double defaultVal) {
double value;
if (EvalResult.testDefaultDouble(r)) {
if (EvalResult.testDefaultDouble(h) || EvalResult.testDefaultDouble(l)) {
value = defaultVal;
} else {
value = (h + l) / 2;
}
} else {
value = r;
}
return value;
}
public static double getRockFrags(Horizon horizon) {
double rockFrag = 0;
for (Fragments f : horizon.fragments.values()) {
rockFrag += averageMissingValue(f.fragvol_h(), f.fragvol_l(), f.fragvol_r());
}
return rockFrag / 100;
}
public static double getSoilpH(Horizon horizon) {
double soilPh = !EvalResult.testDefaultDouble(horizon.ph1to1h2o_r()) ? horizon.ph1to1h2o_r() : horizon.ph01mcacl2_r();
if (soilPh == 0.0 || EvalResult.testDefaultDouble(soilPh)) {
soilPh = 7;
}
return soilPh;
}
public static double getCationExchangeCapacity(Horizon horizon) {
double cec7 = !EvalResult.testDefaultDouble(horizon.cec7_r()) ? horizon.cec7_r() : horizon.ecec_r();
if (cec7 == 0 || EvalResult.testDefaultDouble(cec7)) {
cec7 = horizon.claytotal_r() * 100 * 0.5 + horizon.om_r() * 100 * 2.0;
cec7 = checkBoundaries(cec7, 0.0, 400);
}
horizon.cec7_r(cec7);
return cec7;
}
public static double checkBoundaries(double val, double low, double high) {
if (val < low) {
val = low;
}
if (val > high) {
val = high;
}
return val;
}
public static String formatDouble(double val, int numdig) {
if (Double.isNaN(val)) {
return "NaN";
}
if (val < 0.001 && val > -0.001 && val != 0.0) {
String tmp = Double.toString(val);
int pdx = tmp.indexOf('.');
int edx = tmp.indexOf('E');
if (edx == -1) {
edx = tmp.indexOf('e');
}
if ((pdx + numdig) >= edx) {
return tmp;
}
return tmp.substring(0, pdx + numdig + 1) + tmp.substring(edx);
}
NumberFormat nf = NumberFormat.getNumberInstance(Locale.US);
java.text.FieldPosition fp = new java.text.FieldPosition(NumberFormat.INTEGER_FIELD);
nf.setMaximumFractionDigits(numdig);
nf.setMinimumFractionDigits(numdig);
nf.setGroupingUsed(false);
return nf.format(val, new StringBuffer(), fp).toString();
}
public static String formatIFCValue(String nam, String val) {
return println("# " + nam) + println(val);
}
public static String formatIFCValue(String nam, int val) {
return println("# " + nam) + println(Integer.toString(val));
}
public static String formatIFCValue(String nam, double val, int numdig) {
String ret_val = println("# " + nam);
if (EvalResult.testDefaultDouble(val)) {
throw new ArithmeticException("Soil property not found in one or more layers - " + nam);
} else {
ret_val += println(formatDouble(val, numdig));
}
return ret_val;
}
public static String formatIFCArray(int numLayers, String nam, double arr[], int numdig) {
String ret_val = println("# " + nam);
for (int ldx = 0; ldx < numLayers; ldx++) { // check to make sure nsl is correct
if (Double.isNaN(arr[ldx])) {
throw new ArithmeticException("Soil property not found in one or more layers - " + nam);
} else {
ret_val += formatDouble(arr[ldx], numdig) + " ";
}
}
ret_val += println("");
return ret_val;
}
public static String formatIFCArray(String name, Component comp, int numdig, Function<Horizon, Double> getter) {
String ret_val = println("# " + name);
for (Horizon horizon : comp.horizons.values()) {
double value = getter.apply(horizon);
if (EvalResult.testDefaultDouble(value)) {
throw new ArithmeticException("Soil property not found in one or more layers - " + name);
} else {
ret_val += formatDouble(value, numdig) + " ";
}
}
ret_val += println("");
return ret_val;
}
public static String println(String line) {
return line + System.lineSeparator();
}
public static void convertUnits(MapUnit mapUnit, Component comp) {
mapUnit.brockdepmin(mapUnit.brockdepmin() * 10);
comp.calculated_resdeptmin_r(comp.calculated_resdeptmin_r() * 10);
comp.slope_r(comp.slope_r() / 100);
for (Horizon h : comp.horizons.values()) {
h.hzdepb_r(metricConversion(h.hzdepb_r(), Metric.CM, Metric.MM));
h.hzdept_r(metricConversion(h.hzdept_r(), Metric.CM, Metric.MM));
h.hzthk_r(metricConversion(h.hzthk_r(), Metric.CM, Metric.MM));
h.sandtotal_r(pctToFraction(h.sandtotal_r()));
h.claytotal_r(pctToFraction(h.claytotal_r()));
h.silttotal_r(pctToFraction(h.silttotal_r()));
h.sandvc_r(pctToFraction(h.sandvc_r()));
h.sandco_r(pctToFraction(h.sandco_r()));
h.sandmed_r(pctToFraction(h.sandmed_r()));
h.sandfine_r(pctToFraction(h.sandfine_r()));
h.sandvf_r(pctToFraction(h.sandvf_r()));
h.om_r(pctToFraction(h.om_r()));
h.caco3_r(pctToFraction(h.caco3_r()));
}
}
public static double getSurfaceAlbedo(double surfaceAlbedo, double organicMaterial) {
if (surfaceAlbedo == 0.0 || EvalResult.testDefaultDouble(surfaceAlbedo) || surfaceAlbedo < 0 || surfaceAlbedo > 1) {
surfaceAlbedo = 0.6 / (Math.exp(0.4 * organicMaterial * 100.0));
}
return surfaceAlbedo;
}
/**
* This variable in the WEPS code is only set by the GUI but I created this in
* case we hear otherwise
*
* @return 0
*/
public static double getSurfaceFragmentCover() {
return 0;
}
public static double getLinearExtensibility(Horizon horizon) {
double lep = horizon.lep_r();
if (lep == 0 || EvalResult.testDefaultDouble(lep)) {
double bsd = WEPSUtils.estimateSettledBulkDensity(horizon.claytotal_r(), horizon.sandtotal_r(), horizon.om_r());
if (horizon.wthirdbar_r() > bsd) {
bsd = horizon.wthirdbar_r();
}
lep = WEPSUtils.estimateLinearExtensibility(bsd, horizon.wthirdbar_r());
lep = checkBoundaries(lep, 0, 30);
}
horizon.lep_r(lep);
return lep;
}
public static double getAggregateMeanDiameter(Horizon horizon) {
double mean = Math.exp(1.343 - 2.235 * horizon.sandtotal_r()
- // aggregate geometric mean diameter
1.226 * horizon.silttotal_r() - 0.0238 * horizon.sandtotal_r() / horizon.claytotal_r()
+ 33.6 * horizon.om_r() + 6.85
* horizon.caco3_r()) * (1 + 0.006 * horizon.hzdepb_r());
return checkBoundaries(mean, 0.03, 30.0);
}
public static double getAggregateGeometricMeanDiameter(Horizon horizon){
return (1.226 * horizon.silttotal_r() - 0.0238 * horizon.sandtotal_r() / horizon.claytotal_r()
+ 33.6 * horizon.om_r() + 6.85
* horizon.caco3_r());
}
public static double getAggregateStdDeviation(Horizon horizon, double mean) {
if (EvalResult.testDefaultDouble(mean)) {
mean = getAggregateMeanDiameter(horizon);
}
double deviation = 1 / (0.012448 + 0.002463 * mean + (0.093467 / (Math.pow(mean, 0.5))));
return checkBoundaries(deviation, 1.0, 20.0);
}
public static double getAggregateStdDeviation(Horizon horizon) {
return getAggregateStdDeviation(horizon, Double.NaN);
}
public static double getMaxAggregateSize(Horizon horizon) {
double mean = getAggregateMeanDiameter(horizon);
double dev = getAggregateStdDeviation(horizon, mean);
double maxsize = Math.pow(dev, (1.52 * (Math.pow(mean, -0.449)))) * mean;
return checkBoundaries(maxsize, 1.0, 1000.0);
}
/**
* The WEPS code set this value to 0.01.
*
* @return
*/
public static double getMinAggregateSize() {
return 0.01;
}
/**
* The WEPS code sets this to 1.8 It then checks to make sure it is between
* 0.6 and 2.5?
*
* @return
*/
public static double getAggregateDensity(Horizon horizon) {
return 1.8;
}
public static double getAggregateStability(Horizon horizon) {
double stability = (horizon.claytotal_r() > 0.5) ? 2.73 : 0.83 + 15.7 * horizon.claytotal_r()
- 23.8 * Math.pow(horizon.claytotal_r(), 2);
return checkBoundaries(stability, 0.1, 7.0);
}
/**
* The WEPS code sets this to 0.01
*
* @return
*/
public static double getCrustThickness() {
return 0.01;
}
public static double getCrustDensity(Horizon horizon) {
return getAggregateDensity(horizon);
}
public static double getCrustStability(Horizon horizon) {
return getAggregateStability(horizon);
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getCrustFraction() {
return 0.0;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getCrustLooseMaterialMass() {
return 0.0;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getCrustLooseMaterialFraction() {
return 0.0;
}
/**
* WEPS code sets this to 4
*
* @return
*/
public static double getRandomRoughness() {
return 4;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getRoughnessOrientation() {
return 0.0;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getRoughnessHeight() {
return 0.0;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getRoughnessSpacing() {
return 10.0;
}
/**
* Only set by WEPS GUI
*
* @return
*/
public static double getRoughnessWidth() {
return 10.0;
}
public static double getInitialSWC(Horizon horizon) {
return (getFieldCapacitySWC(horizon) + getWiltingPointSWC(horizon)) / 2;
}
public static double getSaturatedSWC(Horizon horizon) {
double sax_x = 0.332, sax_y = -7.251e-4, sax_z = 0.1276; // coef used to calc THETAS
double saturated = sax_x + sax_y * horizon.sandtotal_r() * 100.
+ // saturated soil water content
sax_z * (Math.log(horizon.claytotal_r() * 100.) / Math.log(10.));
return checkBoundaries(saturated, 0.208, 0.550);
}
public static double getSoilCB(Horizon horizon) {
double sax_e = -3.140, sax_f = -2.22e-3, sax_g = -3.484e-5; // coef used to calc CB
double cb = -(sax_e
+ // cb(k) == ifc3_layer[][23]
(sax_f * Math.pow(horizon.claytotal_r() * 100., 2.0))
+ (sax_g * Math.pow(horizon.sandtotal_r() * 100., 2.0) * horizon.claytotal_r() * 100.));
return checkBoundaries(cb, 0.917, 27.027);
}
public static double getFieldCapacitySWC(Horizon horizon) {
double capacity = pctToFraction(horizon.wthirdbar_r());
if (force || capacity == 0.0 || EvalResult.testDefaultDouble(capacity)) {
capacity = Math.pow((33.0 / get_sax_a(horizon)), (1 / -getSoilCB(horizon)));
capacity = checkBoundaries(capacity, 0.012, .0335);
}
return capacity;
}
public static double getWiltingPointSWC(Horizon horizon) {
double wilting = pctToFraction(horizon.wfifteenbar_r());
if (force || wilting == 0.0 || EvalResult.testDefaultDouble(wilting)) {
wilting = Math.pow((1500.0 / get_sax_a(horizon)), (1 / -getSoilCB(horizon)));
wilting = checkBoundaries(wilting, 0.005, 0.242);
}
return wilting;
}
public static double getAirEntryPotential(Horizon horizon) {
double aep = -get_sax_a(horizon) * Math.pow(getSaturatedSWC(horizon), -getSoilCB(horizon));
return checkBoundaries(aep, -17.91, 0.0);
}
public static double getSaturatedHydraulicConductivity(Horizon horizon) {
double shc = horizon.ksat_r() / 1000000; // cvt from micro_meter/sec to m/sec
if (!sdmksatOnly) {
if (force || shc == 0.0 || EvalResult.testDefaultDouble(shc)) {
double sk = 9.6 - 0.81 * Math.log(horizon.silttotal_r() * 100.) / Math.log(10.)
- 1.09 * Math.log(horizon.claytotal_r() * 100.) / Math.log(10.) - 4.64 * horizon.dbthirdbar_r();
shc = Math.pow(10, sk) * 0.24 * (1.0 / 86400);
shc = checkBoundaries(shc, 0.0, 0.001);
}
}
return shc;
}
private static double get_sax_a(Horizon horizon) {
return 100 * Math.exp(-4.396
- // calc A factor for pote
0.0715 * (horizon.claytotal_r() * 100.)
- // using two steps because that is how
(4.880e-4) * Math.pow(horizon.sandtotal_r() * 100., 2)
- // saxton95 does it
(4.285e-5) * Math.pow(horizon.sandtotal_r() * 100., 2) * (horizon.claytotal_r() * 100.));
}
public static void averageStratifiedLayers(Component comp, StratifiedTextures st) {
boolean first = true;
Layers:
for (Horizon h : comp.horizons.values()) {
for (TextureGroup t : h.textureGroups.values()) {
if (t.stratextsflag() && t.rvindicator()
&& t.texture().trim().toUpperCase().startsWith("SR")) { // stratified key
String[] textures = t.texture().split(" ", -1);
first = true;
double sand = 0;
double clay = 0;
double vsf = 0;
int count = 0;
for (String subTexture : textures) {
subTexture = subTexture.trim();
if (first) {
//skip the sr marker
first = false;
continue;
}
if (st.isTextureValid(subTexture)) {
count++;
sand += st.getSand(subTexture);
sand /= count;
clay += st.getClay(subTexture);
clay /= count;
vsf += st.getVeryFineSand(subTexture);
vsf /= count;
} else {
//can't fix up this layer
continue Layers;
}
}
if (sand + clay > 1) {
//TODO: log error
continue Layers;
}
double silt = 1 - sand - clay;
//if(EvalResult.testDefaultDouble(h.sandtotal_r()))
h.sandtotal_r(sand * 100);
//if(EvalResult.testDefaultDouble(h.silttotal_r()))
h.silttotal_r(silt * 100);
//if(EvalResult.testDefaultDouble(h.claytotal_r()))
h.claytotal_r(clay * 100);
//if(EvalResult.testDefaultDouble(h.sandvf_r()))
h.sandvf_r(vsf * 100);
}
}
}
}
public static double testNumericSoilElement(double high, double low, double r, double defaultValue) {
if (EvalResult.testDefaultDouble(r)) {
r = defaultValue;
}
if (r < low) {
r = defaultValue;
}
if (r > high) {
r = defaultValue;
}
return r;
}
public static int testNumericSoilElement(int high, int low, int r, int defaultValue) {
if (EvalResult.testDefaultInteger(r)) {
r = defaultValue;
}
if (r < low) {
r = defaultValue;
}
if (r > high) {
r = defaultValue;
}
return r;
}
public static double metricConversion(double value, Metric from, Metric to) {
return value / (to.value / from.value);
}
public static double pctToFraction(double value) {
return value / 100;
}
public enum Metric {
MM(.001), CM(.01), DM(.1), M(1), DKM(10), HM(100), KM(1000);
private final double value;
private Metric(double value) {
this.value = value;
}
public double getValue() {
return value;
}
}
}