V2_2.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 java.util.Map;
import javax.ws.rs.Path;
import m.ghg.ApplicationResources;
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.utils.EvalResult;

/**
 *
 * @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/2.2")
@Resource(from = ApplicationResources.class)
@Resource(from = soils.db.DBResources.class)
public class V2_2 extends ModelDataService {

  protected static final String DEFAULT_DAYCENT_SOIL_FILE = "soils.in";
  protected static final String STREAM_FILE_KEY = "stream_output_file";
  protected static final String SOIL_FILE_KEY = "soil_file";

  protected String aoaId;
  protected int cokey;
  protected JSONObject aoa_geometry;
  protected boolean fileAsJSON = false;
  boolean ignoreBadCaco3;
  boolean calibonly;

  protected String soilIN;

  @Override
  protected void preProcess() throws ServiceException {
    aoaId = parameter().getString(soils.AoA.AOA_ID, "");
    cokey = parameter().getInt("cokey", 0);
    //aoa_geometry = parameter().getJSON(soils.AoA.AOA_GEOMETRY);
    fileAsJSON = parameter().getBoolean(STREAM_FILE_KEY, false);
    ignoreBadCaco3 = parameter().getBoolean("ignore_bad_caco3", false);
    calibonly = parameter().getBoolean("calibration", 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)) {

      V2_2.AoA soilData = new AoA(soilsDb, gisDb, crDb, cokey, LOG, ignoreBadCaco3);
      soilIN = getSoilInFile(soilData.getSoilINComp(), soilData.getSoilInCompMapUnit());

    }
  }

  @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, MapUnit mapUnit) throws Exception {
    DaycentStrata daycentStrata = new DaycentStrata(mapUnit, comp);
    Map<Double, List<Profile>> strata = daycentStrata.computeDaycentStrata();

    SoilFileBuilder2 fileBuild = new SoilFileBuilder2(strata, calibonly);
    return fileBuild.build();
  }

  /**
   *
   */
  static class AoA extends soils.AoA {

    Connection crDb;
    double minimumPercentage = 0.0;
    List<String> removedAreaToSmall = new ArrayList<>();
    String soilsIN;
    Component soilINComponent = null;
    MapUnit mapUnit = null;

    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) && (shape.getType() != GISObjectType.multipolygon)) {//  Also returns this for multipolygons.
        throw new ServiceException("The GHG soils service only accepts Polygon or MultiPolygon shapes");
      }

      shape.makeValid(GISObject.UsePurpose.all_purposes, GISObject.GISType.all_types);
      area = shape.areaInAcres();
    }

    AoA(SOILS_DATA soilsDb, Connection gisDb, Connection crDb,
        int cokey, SessionLogger LOG, boolean ignoreBadCaco3) throws Exception {

      super(soilsDb, LOG);

      this.crDb = crDb;

      MapUnit tMapUnit = null;
      Component comp = new Component();
      if (!soilsDb.validateComponent(cokey)) {
        throw new ServiceException("Invalid cokey, (" + cokey + "), specified.");
      }
      comp.cokey(cokey);

      if (soilsDb.validateComponent(Integer.parseInt(comp.cokey()))) {
        try {
          if (null == (tMapUnit = soilsDb.findWEPSDataByCokey(comp, false, false, false))) {
            throw new ServiceException("Could not retrieve the proper SDM soils data for this cokey");
          }
        } catch (ServiceException ex) {
          if (!ex.getMessage().contains("caco3")) {
            throw ex;
          } else {
            if (ex.getMessage().indexOf("is empty") != ex.getMessage().lastIndexOf("is empty")) // More than one missing value{
            {
              throw ex;
            } else {
              if (!ignoreBadCaco3) {
                throw ex;
              }
            }
          }
        }
        if (null == tMapUnit) {
          throw new ServiceException("Could not retrieve the proper SDM soils data for this cokey");
        }
        map_units.put(tMapUnit.mukey(), tMapUnit);
        soilINComponent = comp;
        this.mapUnit = tMapUnit;
      } else {
        throw new ServiceException("The cokey " + cokey + " provided does not exist in the SDM database.  Cannot continue.");
      }

    }

    public Component getSoilINComp() {
      return soilINComponent;
    }

    public MapUnit getSoilInCompMapUnit() {
      return mapUnit;
    }

    /**
     *
     * @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;
          }

          //  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;
    }

  }

}