SWATCData.java [src/java/d/dataNodes] 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 d.dataNodes;
import java.io.IOException;
import java.io.PrintWriter;
import java.nio.file.Files;
import java.nio.file.Paths;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.List;
import java.util.Iterator;
import java.util.stream.Collectors;
import soils.Component;
import soils.Horizon;
import soils.exceptions.SDMException;
import soils.utils.SoilUtils;
/**
*
* @author geter
*/
public class SWATCData {
static DecimalFormat df = new DecimalFormat("0.00");
public void computeAndWrite(String workspace, String filename, Component comp) throws IOException {
checkForNaN(comp);
String header = buildHeader(comp);
List<SoilLayerData> l = getSoilLayerData(comp);
writeFile(header, l, workspace, filename);
}
private void writeFile(String header, List<SoilLayerData> l, String workspace, String filename) throws IOException {
try (PrintWriter pw = new PrintWriter(
Files.newBufferedWriter(Paths.get(workspace, filename)))) {
pw.println(header);
// pw.println(singleLine);
outputSoilLayerData(l, pw);
}
}
private double computePorosity(Component comp) {
List<Porosity> weightedPorosity = new ArrayList();
comp.horizons.keySet().stream()
.map((key) -> comp.horizons.get(key))
.forEachOrdered((h) -> {
double thickness = h.hzthk_r();
double bulkdensity = h.dbthirdbar_r();
double organicmatter = h.om_r();
double iron = h.getTableHorizon().freeiron_r();
weightedPorosity.add(new Porosity(thickness, bulkdensity, organicmatter, iron));
});
double thickness_sum = 0;
thickness_sum = comp.horizons.entrySet().stream()
.map((map) -> map.getValue().hzthk_r())
.reduce(thickness_sum, (accumulator, _item) -> accumulator + _item);
double porosity_sum = 0;
porosity_sum = weightedPorosity.stream()
.map(wp -> wp.getWeightedPorosity())
.reduce(porosity_sum, (accumulator, _item) -> accumulator + _item);
return (Double.isNaN(porosity_sum)) ? 0.5 : porosity_sum / thickness_sum;
}
private String buildHeader(Component comp) {
String name = comp.compname();
int nly = comp.horizons.values().size();
String hyd_grp = comp.hydgrp();
double dp_tot = comp.profileDepth();
double anion_excl = computePorosity(comp); // OPTIONAL porosity fraction from which anions are excluded
double perc_crk; // OPTIONAL crack volume potential of soil
List<String> texture_list = new ArrayList(); // texture
comp.horizons.keySet().stream()
.map((key) -> comp.horizons.get(key))
.forEachOrdered((h) -> {
try {
texture_list.add(h.getRepresentativeTextureName());
} catch (SDMException ex) {
throw new RuntimeException(ex.getMessage());
}
});
return String.format("%-15s", "swatc sol file\n")
+ String.format("%-12s", "Soil Name : ") + String.format("%-16s", name) + "\n"
+ String.format("%-24s", "Soil Hydrologic Group : ") + String.format("%-1s", hyd_grp) + "\n"
+ String.format("%-28s", "Maximum rooting depth [mm] : ") + String.format("%-12s", df.format(cm2mm(dp_tot))) + "\n"
+ String.format("%-51s", "Porosity fraction from which anions are excluded : ") + String.format("%-5.3f", anion_excl) + "\n"
+ String.format("%-33s", "Crack volume potential of soil : ") + "\n"
+ String.format("%-28s", "Texture 1 : ") + String.format("%-18s", texture_list.stream().collect(Collectors.joining("-")));
/* + String.format("%18s", "nly")
+ String.format("%18s", "hyd_grp")
+ String.format("%12s", "dp_tot")
+ String.format("%16s", "anion_excl")
+ String.format("%14s", "perc_crk")
+ String.format("%2s", " ")
+ String.format("%-18s", "texture")
+ String.format("%21s", "dp")
+ String.format("%14s", "bd")
+ String.format("%14s", "awc")
+ String.format("%14s", "soil_k")
+ String.format("%14s", "carbon")
+ String.format("%14s", "clay")
+ String.format("%14s", "silt")
+ String.format("%14s", "sand")
+ String.format("%14s", "rock")
+ String.format("%14s", "alb")
+ String.format("%14s", "usle_k")
+ String.format("%14s", "ec")
+ String.format("%14s", "caco3")
+ String.format("%14s", "ph");
*/
}
private void outputSoilLayerData(List<SoilLayerData> l, PrintWriter pw){
Iterator l_iterator = l.iterator();
String line = "Depth [mm] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.dp);
}
line += "\n";
l_iterator = l.iterator();
line += "Bulk Density Moist [g/cc] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.bd);
}
line += "\n";
l_iterator = l.iterator();
line += "Ave. AW Incl. Rock Frag :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.awc);
}
line += "\n";
l_iterator = l.iterator();
line += "Ksat. (est.) [mm/hr] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.soil_k);
}
line += "\n";
l_iterator = l.iterator();
line += "Organic Carbon [weight %] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.carbon);
}
line += "\n";
l_iterator = l.iterator();
line += "Clay [weight %] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.clay);
}
line += "\n";
l_iterator = l.iterator();
line += "Silt [weight %] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.silt);
}
line += "\n";
l_iterator = l.iterator();
line += "Sand [weight %] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.sand);
}
line += "\n";
l_iterator = l.iterator();
line += "Rock Fragments [vol. %] :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.rock);
}
line += "\n";
l_iterator = l.iterator();
line += "Soil Albedo (Moist) :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.alb);
}
line += "\n";
l_iterator = l.iterator();
line += "USLE Soil Erosion K :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.usle_k);
}
line += "\n";
l_iterator = l.iterator();
line += "Salinity (EC, Form 5) :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.ec);
}
line += "\n";
l_iterator = l.iterator();
line += "Soil pH :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.ph);
}
line += "\n";
l_iterator = l.iterator();
line += "Soil CACO3 :";
while(l_iterator.hasNext()){
SoilLayerData layer = (SoilLayerData) l_iterator.next();
line += String.format("%12.2f", layer.caco3);
}
pw.print(line);
}
private List<SoilLayerData> getSoilLayerData(Component comp) {
double albedo = comp.albedodry_r();
return comp.horizons.values().stream()
.map(h -> new SoilLayerData(h, albedo)).collect(Collectors.toList());
}
class Porosity {
double thickness;
double bulkdensity;
double iron;
double organic_carbon;
double porosity;
Porosity(double thickness, double bulkdensity, double organic_matter, double iron) {
this.thickness = thickness;
this.bulkdensity = bulkdensity;
this.organic_carbon = omatter2ocarbon(organic_matter);
this.iron = iron;
compute();
}
// ch 618.43 - p. 618-A.43
private double omatter2ocarbon(double organic_matter) {
return organic_matter / 1.724;
}
private void compute() {
double dp1 = 1.4;
double dp2 = 4.2;
double dp3 = 2.65;
double denominator1 = (1.7 * organic_carbon) / dp1;
double denominator2 = (1.6 * iron) / dp2;
double denominator3 = (100 - ((1.7 * organic_carbon) + (1.6 * iron))) / dp3;
// ch 618.45 - p. 618-A.44
double particledensity = 100 / (denominator1 + denominator2 + denominator3);
porosity = 100 * (1 - (bulkdensity / particledensity));
}
public double getWeightedPorosity() {
return thickness * porosity;
}
}
class SoilLayerData {
public double dp; // depth [mm]
public double bd; // bulk density [g/cc] -> check
public double awc; // Available water content
public double soil_k; // [mm/hr] -> check
public double carbon; // [%]
public double clay; // [%]
public double silt; // [%]
public double sand; // [%]
public double rock; // [%]
public double alb;
public double usle_k;
public double ec; // salinity is often measured by electrical conductivity
public double caco3;
public double ph;
public SoilLayerData(Horizon h, double albedo) {
dp = cm2mm(h.hzdepb_r()); // depth [mm]
bd = h.dbthirdbar_r(); // bulk density [g/cc] -> check
awc = availableWaterCapacity(h);
soil_k = h.ksat_r(); // [mm/hr] -> check
carbon = om2oc(h.om_r()); // [%]
clay = h.claytotal_r(); // [%]
silt = h.silttotal_r(); // [%]
sand = h.sandtotal_r(); // [%]
rock = h.fragvol_r(); // [%]
alb = albedo;
usle_k = h.kffact();
ec = h.ec_r(); // [dS/m]
caco3 = h.caco3_r(); // [%]
ph = h.ph1to1h2o_r(); // [%]
}
// SDM Handbook - 618.6 Available Water Capacity
private double availableWaterCapacity(Horizon h) {
return (h.wthirdbar_r() - h.wfifteenbar_r()) * h.dbthirdbar_r() / 100;
}
}
private static double cm2mm(double val) {
return val * 10;
}
private static double micropersec2mmperhour(double val) {
return val * 3.6;
}
// SDM Handbook - p. 618-A.43
private static double om2oc(double om) {
return om / 1.724;
}
private void checkForNaN(Component comp) {
if (Double.isNaN(comp.albedodry_r())) {
comp.albedodry_r(SoilUtils.averageMissingValue(comp.albedodry_h(), comp.albedodry_l(), comp.albedodry_r()));
}
comp.horizons.keySet().stream()
.map((key) -> comp.horizons.get(key))
.forEachOrdered((h) -> {
if (Double.isNaN(h.dbthirdbar_r())) {
h.dbthirdbar_r(SoilUtils.averageMissingValue(
h.dbthirdbar_h(),
h.dbthirdbar_l(),
h.dbthirdbar_r()));
}
if (Double.isNaN(h.om_r())) {
h.om_r(SoilUtils.averageMissingValue(
h.om_h(),
h.om_l(),
h.om_r()));
}
if (Double.isNaN(h.getTableHorizon().freeiron_r())) {
h.getTableHorizon().freeiron_r(SoilUtils.averageMissingValue(
h.getTableHorizon().freeiron_h(),
h.getTableHorizon().freeiron_l(),
h.getTableHorizon().freeiron_r()));
}
if (Double.isNaN(h.hzdepb_r())) {
h.hzdept_r(SoilUtils.averageMissingValue(
h.hzdepb_h(),
h.hzdepb_l(),
h.hzdepb_r()));
}
if (Double.isNaN(h.ksat_r())) {
h.ksat_r(SoilUtils.averageMissingValue(
h.ksat_h(),
h.ksat_l(),
h.ksat_r()));
}
if (Double.isNaN(h.claytotal_r())) {
h.claytotal_r(SoilUtils.averageMissingValue(
h.claytotal_h(),
h.claytotal_l(),
h.claytotal_r()));
}
if (Double.isNaN(h.silttotal_r())) {
h.silttotal_r(SoilUtils.averageMissingValue(
h.silttotal_h(),
h.silttotal_l(),
h.silttotal_r()));
}
if (Double.isNaN(h.sandtotal_r())) {
h.sandtotal_r(SoilUtils.averageMissingValue(
h.sandtotal_h(),
h.sandtotal_l(),
h.sandtotal_r()));
}
if (Double.isNaN(h.ph1to1h2o_r())) {
h.ph1to1h2o_r(SoilUtils.averageMissingValue(
h.ph1to1h2o_h(),
h.ph1to1h2o_l(),
h.ph1to1h2o_r()));
}
if (Double.isNaN(h.caco3_r())) {
h.caco3_r(SoilUtils.averageMissingValue(
h.caco3_h(),
h.caco3_l(),
h.caco3_r()));
}
if (Double.isNaN(h.wthirdbar_r())) {
h.wthirdbar_r(SoilUtils.averageMissingValue(
h.wthirdbar_h(),
h.wthirdbar_l(),
h.wthirdbar_r()));
}
if (Double.isNaN(h.wfifteenbar_r())) {
h.wfifteenbar_r(SoilUtils.averageMissingValue(
h.wfifteenbar_h(),
h.wfifteenbar_l(),
h.wfifteenbar_r()));
}
});
}
}