V1_0.java [src/java/d/soils] 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.soils;
import csip.Config;
import csip.ModelDataService;
import csip.api.server.PayloadResults;
import csip.api.server.ServiceException;
import csip.SessionLogger;
import csip.annotations.Description;
import csip.annotations.Name;
import csip.annotations.Resource;
import gisobjects.GISObject;
import gisobjects.GISObject.GISObjectType;
import gisobjects.GISObjectException;
import static gisobjects.GISObjectFactory.createGISObject;
import static gisobjects.db.GISEngineFactory.createGISEngine;
import java.sql.Connection;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import javax.ws.rs.Path;
import static m.ghg.ApplicationResources.GISDB_SQLSVR;
import static m.ghg.ApplicationResources.STEP_CRDB_SVR;
import org.codehaus.jettison.json.JSONArray;
import org.codehaus.jettison.json.JSONException;
import org.codehaus.jettison.json.JSONObject;
import soils.Component;
import soils.Horizon;
import soils.MapUnit;
import soils.SoilsData;
import soils.db.SOILS_DATA;
import soils.db.SOILS_DB_Factory;
import soils.db.tables.TableComponent;
import soils.db.tables.TableComponentCalculations;
import soils.db.tables.TableHorizon;
import soils.db.tables.TableHorizonCalculations;
import soils.db.tables.TableMapUnit;
import soils.db.tables.TableMapUnitCalculations;
import soils.exceptions.WEPPException;
import soils.utils.EvalResult;
import soils.utils.SoilUtils;
/**
*
* @author <a href="mailto:shaun.case@colostate.edu">Shaun Case</a>
*/
/*
SOILS_IN File specifications:
==============================================================================
SOILS.IN example:
0.0 2.0 1.44 0.31092 0.13578 0.80 0.01 0.39 0.28 0.01 0.11 0.00027 5.00
2.0 5.0 1.44 0.31092 0.13578 0.20 0.04 0.39 0.28 0.01 0.08 0.00027 5.00
5.0 10.0 1.44 0.31092 0.13578 0.00 0.25 0.39 0.28 0.01 0.05 0.00027 5.00
10.0 20.0 1.44 0.31092 0.13578 0.00 0.30 0.39 0.28 0.01 0.01 0.00027 5.00
20.0 30.0 1.44 0.31092 0.13578 0.00 0.10 0.39 0.28 0.01 0.00 0.00027 5.00
30.0 45.0 1.44 0.31092 0.13578 0.00 0.05 0.39 0.28 0.01 0.00 0.00027 5.00
45.0 60.0 1.44 0.31092 0.13578 0.00 0.04 0.39 0.28 0.01 0.00 0.00027 5.00
60.0 75.0 1.44 0.31092 0.13578 0.00 0.03 0.39 0.28 0.01 0.00 0.00027 5.00
75.0 90.0 1.44 0.31092 0.13578 0.00 0.02 0.39 0.28 0.01 0.00 0.00027 5.00
90.0 105.0 1.44 0.31092 0.13578 0.00 0.01 0.39 0.28 0.01 0.00 0.00027 5.00
105.0 120.0 1.44 0.31092 0.13578 0.00 0.00 0.39 0.28 0.01 0.00 0.00027 5.00
120.0 150.0 1.44 0.31092 0.13578 0.00 0.00 0.39 0.28 0.01 0.00 0.00027 5.00
Column 1 - Minimum depth of soil layer (cm)
Column 2 - Maximum depth of soil layer (cm)
Column 3 - Bulk density of soil layer (g/cm^3)
Column 4 - Field capacity of soil layer, volumetric
Column 5 - Wilting point of soil layer, volumetric
Column 6 - Evaporation coefficient for soil layer (currently not being used) ??????
Column 7 - Percentage of roots in soil layer, these values must sum to 1.0 ??????
Column 8 - Fraction of sand in soil layer, 0.0 - 1.0
Column 9 - Fraction of clay in soil layer, 0.0 - 1.0
Column 10 - Organic matter in soil layer, fraction 0.0 - 1.0
Column 11 - Minimum volumetric soil water content below wilting point for soil ??????
layer, soil water content will not be allowed to drop below this
value
Column 12 - Saturated hydraulic conductivity of soil layer in centimeters per ??? CONVERTED FROM mm/h ??? (/10/60)
second
Column 13 - pH of soil layer
NOTES:
Percentage of silt for soil layer is computed as follows:
percent silt = (1.0 - (percent sand + percent clay))
For the trace gas subroutines it is currently recommended to use the following
layering structure for the top 3 soil layers in your soils.in file:
layer 1 - 0.0 cm to 2.0 cm
layer 2 - 2.0 cm to 5.0 cm
layer 3 - 5.0 cm to 10.0 cm
The depth structure in this file should match the ADEP values in the FIX.100
file in such a way that the boundaries for the soil layer depths can be
matched with the ADEP values. For example, using the file above and ADEP
values of 10, 20, 15, 15, 30, 30, 30, 30, 30, and 30:
layers 1, 2 and 3 match the first 10 centimeter ADEP value
layers 4 and 5 match the second 20 centimeter ADEP value
layer 6 matches the third 15 centimeter ADEP value
layer 7 matches the fourth 15 centimeter ADEP value
layers 8 and 9 match the first 30 centimeter ADEP value
layers 10 and 11 match the second 30 centimeter ADEP value
layer 12 matches the third 30 centimeter ADEP value
The value for NLAYER in the <site>.100 file should be set to match the number
of ADEP values that you are using when you match the layering to the soils.in
file. For the example above NLAYER should be set to 7.
==============================================================================
*/
@Name("DayCent IN soil file")
@Description("")
@Path("d/insoilfile/1.0")
@Resource(from = m.ghg.ApplicationResources.class)
@Resource(from = soils.db.DBResources.class)
public class V1_0 extends ModelDataService {
protected static final String DEFAULT_DAYCENT_SOIL_FILE = "dayCentSoil.IN";
protected static final String STREAM_FILE_KEY = "stream_output_file";
protected static final String SOIL_FILE_KEY = "soil_file";
protected String aoaId;
protected JSONObject aoa_geometry;
protected boolean fileAsJSON = false;
protected String soilIN;
protected JSONArray excluded;
@Override
protected void preProcess() throws ServiceException {
aoaId = parameter().getString(soils.AoA.AOA_ID, "");
aoa_geometry = parameter().getParamJSON(soils.AoA.AOA_GEOMETRY);
fileAsJSON = parameter().getBoolean(STREAM_FILE_KEY, false);
}
@Override
protected void doProcess() throws Exception {
try (SOILS_DATA soilsDb = SOILS_DB_Factory.createEngine(getClass(), LOG, Config.getString("soils.gis.database.source"));
Connection gisDb = resources().getJDBC(GISDB_SQLSVR);
Connection crDb = resources().getJDBC(STEP_CRDB_SVR)) {
AoA soilData = new AoA(soilsDb, gisDb, crDb, aoa_geometry, LOG);
soilData.getIntersectionsAndComponents();
soilIN = getSoilInFile(soilData.getMajorComponent());
}
}
@Override
protected void postProcess() throws Exception {
if (!fileAsJSON) {
workspace().writeString(DEFAULT_DAYCENT_SOIL_FILE, soilIN);
results().put(workspace().getFile(DEFAULT_DAYCENT_SOIL_FILE), "DayCent soil file");
} else {
results().put(SOIL_FILE_KEY, soilIN, "DayCent soil file data.");
}
}
protected String getSoilInFile(Component comp) {
String outLine = "";
for (Horizon horizon : comp.horizons().values()) {
outLine += buildInLine(horizon);
}
return outLine;
}
/*
Column 1 - Minimum depth of soil layer (cm)
Column 2 - Maximum depth of soil layer (cm)
Column 3 - Bulk density of soil layer (g/cm^3)
Column 4 - Field capacity of soil layer, volumetric
Column 5 - Wilting point of soil layer, volumetric
Column 6 - Evaporation coefficient for soil layer (currently not being used)
Column 7 - Percentage of roots in soil layer, these values must sum to 1.0
Column 8 - Fraction of sand in soil layer, 0.0 - 1.0
Column 9 - Fraction of clay in soil layer, 0.0 - 1.0
Column 10 - Organic matter in soil layer, fraction 0.0 - 1.0
Column 11 - Minimum volumetric soil water content below wilting point for soil
layer, soil water content will not be allowed to drop below this
value
Column 12 - Saturated hydraulic conductivity of soil layer in centimeters per
second
Column 13 - pH of soil layer
*/
protected String buildInLine(Horizon horizon) {
try {
horizon.computeErodibility(); //Set erodability values....may throw an exception if missing values
} catch (WEPPException ex) {
horizon.conductivity(0.0); // Defaulting to 0.0 for now, just for DayCent, if missing values to calculate...discussion?
}
return horizon.hzdept_r() + " " + horizon.hzdepb_r() + " "
+ horizon.dbthirdbar_r() + " " + SoilUtils.getFieldCapacitySWC(horizon) + " "
+ SoilUtils.getWiltingPointSWC(horizon) + " " + "0" + " " + "0" + " "
+ horizon.sandtotal_r() + " " + horizon.claytotal_r() + " "
+ horizon.om_r() + " " + "0" + " " + horizon.getConductivity() / 10.0 / 60.0 + " "//mm/h to cm/s
+ SoilUtils.getSoilpH(horizon) + System.lineSeparator();
}
/**
*
*/
static class AoA extends soils.AoA {
Connection crDb;
double minimumPercentage = 0.0;
List<String> removedAreaToSmall = new ArrayList<>();
AoA(SOILS_DATA soilsDb, Connection gisDb, Connection crDb,
JSONObject geometry, SessionLogger LOG) throws Exception {
super(soilsDb, LOG);
this.crDb = crDb;
shape = createGISObject(geometry, createGISEngine(gisDb));
if (shape.getType() != GISObjectType.polygon) {// Also returns this for multipolygons.
throw new ServiceException("The step thresholds service only accepts Polygon or MultiPolygon shapes");
}
shape.makeValid(GISObject.UsePurpose.all_purposes, GISObject.GISType.all_types);
area = shape.areaInAcres();
}
/**
*
* @param component
*/
public void computeHorizonResults(Component component) {
double profile_thk = 0.0;
boolean haveHzDepth = false;
boolean haveOM_R = false;
double comp_product = 0.0;
int counter = 0;
if (component.horizons.size() > 0) {
double pH = 0.0;
double pH_thickness = 0.0;
if (EvalResult.testDefaultDouble(component.calculated_coarse_frag_vol_total())) {
component.calculated_coarse_frag_vol_total(0.0);
}
for (Horizon horizon : component.horizons.values()) {
double horizonThickness = 0.0;
double horizonProduct = 0.0;
double pH_l = horizon.ph1to1h2o_l();
double pH_r = horizon.ph1to1h2o_r();
double pH_h = horizon.ph1to1h2o_h();
// if (0 == counter) {
// component.calculated_om_r(horizon.om_r());
// }
counter++;
if (EvalResult.testDefaultDouble(horizon.hzthk_r())) {
horizonThickness = horizon.hzdepb_r() - horizon.hzdept_r();
} else {
horizonThickness = horizon.hzthk_r();
}
if ((!haveOM_R) && (!haveHzDepth) && (horizon.selected())) {
component.calculated_om_r(horizon.om_r());
component.calculated_hzdepb_r(horizonThickness);
haveHzDepth = true;
haveOM_R = true;
}
if (EvalResult.testDefaultDouble(pH_r)) {
if (!EvalResult.testDefaultDouble(pH_h) && !EvalResult.testDefaultDouble(pH_l)) {
pH_r = (pH_h + pH_l) / 2.0;
pH_thickness += horizonThickness;
pH += pH_r * horizonThickness;
}//Else ignore this layer in the calculation...don't need to do anything else in this case...
} else {
pH_thickness += horizonThickness;
pH += pH_r * horizonThickness;
}
// if (!haveHzDepth) { //We only want to use the first thickness since these are presorted by the SQL statement in soilsdb
// component.calculated_hzdepb_r(horizonThickness);
// haveHzDepth = true;
// }
// TODO: Is this valid?? What is "coarse" sized fragments?? Shouldn't we be looking at individual fragment records within the horizon for their fragsize_r value ??
component.calculated_coarse_frag_vol_total(component.calculated_coarse_frag_vol_total() + (EvalResult.testDefaultDouble(horizon.fragvol_r()) ? 0 : horizon.fragvol_r()));
profile_thk += horizonThickness;
horizonProduct = horizonThickness * (EvalResult.testDefaultDouble(horizon.fragvol_r()) ? 0 : horizon.fragvol_r());
comp_product += horizonProduct;
}
if (profile_thk > 0.0) {
component.calculated_coarse_frag(comp_product / profile_thk);
}
if (pH_thickness > 0.0) {
component.calculated_pH(pH / pH_thickness);
}
}
}
/**
*
* @throws ServiceException
* @throws GISObjectException
*/
public void getIntersectionsAndComponents() throws ServiceException, GISObjectException {
if (findIntersectedMapUnitsWithAreas()) {
if (soilsDb.findAllBasicComponentHorizonFragTextureData(map_units)) {
for (MapUnit mapUnit : map_units.values()) {
for (Component tComponent : mapUnit.components().values()) {
tComponent.applyHorizonFilter("WATER");
}
}
removeComponentsByAreaSize();
findComponentSoilMoistures();
for (MapUnit mapUnit : map_units.values()) {
if (!mapUnit.isExcluded()) {
for (Component component : mapUnit.components().values()) {
if (!component.isExcluded()) {
computeHorizonResults(component);
}
}
}
}
} else {
throw new ServiceException("Cannot find components for the map units interesected by this AoA");
}
} else {
throw new ServiceException("Cannot find the mapunits intersected by this AoA");
}
}
protected void removeComponentsByAreaSize() {
for (String mukey : map_units.keySet()) {
MapUnit tMapUnit = map_units.get(mukey);
if (!tMapUnit.isExcluded()) {
for (String cokey : tMapUnit.components().keySet()) {
Component tComponent = tMapUnit.components().get(cokey);
if (!tComponent.isExcluded()) {
if (tComponent.area_pct() < this.minimumPercentage) {
removedAreaToSmall.add(cokey);
//tMapUnit.components().remove(cokey);
tComponent.setExclude(true, "Component removed because its area percentage of the AoI is less than the minimum percentage required by the input parameters.");
}
}
}
}
}
}
public void toJSON(PayloadResults output) throws JSONException, ServiceException {
// JSONArray topLevelData = new JSONArray();
// JSONObject map_unit_list = new JSONObject();
// Example: output.putResult(RFACTOR, this.rfactor, "R Factor", null);
JSONArray mapunitData = new JSONArray();
for (MapUnit mapUnit : map_units.values()) {
if (!mapUnit.isExcluded()) {
setOutputs(mapUnit);
//TODO: Maybe just call a toINFile() function...add code to do this.
mapunitData.put(mapUnit.toJSON(true, new ArrayList<>(Arrays.asList("WATER"))));
}
}
output.put(SoilsData.MAP_UNIT_LIST, mapunitData, null, null);
}
// TODO: Remove unneeded outputs from this function.
protected void setOutputs(MapUnit mapUnit) throws ServiceException {
mapUnit.setOutputColumnOrdering(new ArrayList<>(Arrays.asList(TableMapUnit.MUKEY, TableMapUnit.MUNAME,
TableMapUnit.MUSYM, TableMapUnit.AREASYMBOL_NAME, TableMapUnitCalculations.AREA_NAME
)));
mapUnit.setComponentOutputColumnOrdering(new ArrayList<>(Arrays.asList(TableComponent.COKEY, TableComponent.COMPNAME,
TableComponent.COMPPCT_R_NAME, TableComponentCalculations.AREA_PCT_NAME, TableComponentCalculations.COMP_AREA_NAME,
TableComponentCalculations.LONG_NAME,
TableComponent.HYDGRP_NAME, TableComponent.TAXORDER_NAME,
TableComponent.SLOPE_R_NAME,
TableComponentCalculations.NSLP_NAME, TableComponentCalculations.SRP_NAME,
TableComponentCalculations.KFFACT_NAME,
TableComponentCalculations.COARSE_FRAG_TOTAL_NAME, TableComponentCalculations.HWT_LT_24_NAME,
TableComponentCalculations.WTBL_TOP_MIN_NAME
)));
mapUnit.setHorizonOutputColumnOrdering(new ArrayList<>(Arrays.asList(TableHorizon.CHKEY_NAME,
TableHorizon.HZDEPB_R_NAME, TableHorizon.HZDEPT_R_NAME,
TableHorizon.HZTHK_R_NAME, TableHorizon.KFFACT_NAME, TableHorizon.KWFACT_NAME,
TableHorizonCalculations.FRAGVOL_R_NAME
)));
}
private Component getMajorComponent() {
Component maxSize = null;
for (MapUnit mapUnit : map_units.values()) {
if (!mapUnit.isExcluded()) {
for (Component comp : mapUnit.components().values()) {
if (!comp.isExcluded()) {
if (null == maxSize) {
maxSize = comp;
} else {
if (comp.calculated_area() > maxSize.calculated_area()) {
maxSize = comp;
}
}
}
}
}
}
return maxSize;
}
}
}