WindWaterErosion.java [src/java/d/util] 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.util;

import csip.SessionLogger;
import csip.api.server.ServiceException;
import csip.utils.JSONUtils;
import gisobjects.GISObject;
import static gisobjects.GISObject.DEFAULT_ASSUMED_SRID;
import gisobjects.GISObjectException;
import gisobjects.raster.GISRaster;
import gisobjects.raster.interpolation.InterpolationMethodDoubleMean;
import gisobjects.raster.interpolation.InterpolationMethodResult;
import java.io.File;
import java.io.IOException;
import java.sql.Connection;
import java.sql.ResultSet;
import java.sql.Statement;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Map;
import java.util.concurrent.Callable;
import org.codehaus.jettison.json.JSONException;
import org.codehaus.jettison.json.JSONObject;
import org.geotools.coverage.grid.GridCoordinates2D;
import soils.Component;
import soils.MapUnit;
import soils.db.SOILS_DATA;
import soils.db.tables.TableComponent;
import soils.db.tables.TableComponentCalculations;
import soils.db.tables.TableMapUnit;
import soils.db.tables.TableMapUnitCalculations;
import utils.EvalResult;
import utils.MDS_Result;

/**
 *
 * @author <a href="mailto:shaun.case@colostate.edu">Shaun Case</a>
 */
public class WindWaterErosion extends soils.AoA {

  protected boolean outputSoils = false;

  public static final String WWE_TIFF_FILE = "tif_file";
  public static final String BUFFER_SIZE_KEY = "buffer_size";
  static GISRaster gisRaster;
  static final Object GIS_RASTER_LOCK = new Object();

  static final double MINIMUM_PERCENTAGE = 10.0;  // comppct_r values are stored as "percentages" not ratios, so 10% is stored as "10", not "0.10" 
  static final double AREA_SIZE_CUTOFF = 25.0; // The AoI area in acres, below which the minimum_percentage must be reset to 20%
  static final double MINIMUM_PERCENTAGE_2 = 20.0;
  static final int ALLOWED_SRID = DEFAULT_ASSUMED_SRID;
  static final double DEFAULT_WEI = 134.0;
  static final double DEFAULT_POINT_BUFFER_SIZE = 113.5; // 133.5m buffer = 10 acres.

  static final double[] i_factor_list = {
    0,
    38,
    48,
    56,
    86,
    134,
    160,
    180,
    220,
    250,
    310};

  private double bufferSize = DEFAULT_POINT_BUFFER_SIZE;

  /**
   * Singleton.
   *
   * @param rasterFile
   * @throws Exception
   */
  public static synchronized void loadGisRaster(Callable<File> rasterFile) throws Exception {
    if (gisRaster == null) {
      gisRaster = new GISRaster(rasterFile.call());
    }
  }
  double minPercentage = MINIMUM_PERCENTAGE;
  double[] centroid;
  WWErosionResultData resultData = new WWErosionResultData();

  public WindWaterErosion(SOILS_DATA soilsDb, Connection gisDb, Map<String, JSONObject> params, SessionLogger Log) throws Exception {
    super(soilsDb, gisDb, params, Log);
    if (params.containsKey("minimum_component_percentage")) {
      minPercentage = JSONUtils.getDoubleParam(params, "minimum_component_percentage", MINIMUM_PERCENTAGE);
    }

    if (params.containsKey(BUFFER_SIZE_KEY)) {
      bufferSize = JSONUtils.getDoubleParam(params, BUFFER_SIZE_KEY, DEFAULT_POINT_BUFFER_SIZE);
    }

    if (ALLOWED_SRID != getSRID()) {
      throw new ServiceException("Wind and Water Erosion, currently, only supports an input CRS of WGS-84 (4326 SRID).  You requested service for a geometry that was SRID: " + getSRID());
    }

    if (shape.getGeometry().getGeometryType() == "Point") {
      shape = shape.buffer(bufferSize);
    }

    if (shape.getType() != GISObject.GISObjectType.polygon) {
      throw new ServiceException("Wind and Water Erosion services cannot process the shape type specified.  Please use a polygon or multipolygon only.");
    }
    centroid = shape.getCentroid();
    resultData.aoa_area(shape.areaInAcres());

    if (resultData.aoa_area() <= AREA_SIZE_CUTOFF) {
      minPercentage = MINIMUM_PERCENTAGE_2;
    }

    resultData.setOutputColumnOrdering(new ArrayList<String>(Arrays.asList(
        WWErosionResultData.AOA_AREA, WWErosionResultData.AOA_IFACTOR, WWErosionResultData.AOA_KFACTOR, WWErosionResultData.AOA_TFACTOR,
        WWErosionResultData.CFACTOR, WWErosionResultData.CFACTOR_PCT, WWErosionResultData.CFACTOR_SOURCE,
        WWErosionResultData.WATER_KFACTOR, WWErosionResultData.WATER_LSFACTOR, WWErosionResultData.WATER_TFACTOR,
        WWErosionResultData.WIND_EP, WWErosionResultData.WIND_IFACTOR, WWErosionResultData.WIND_TFACTOR,
        WWErosionResultData.DOM_WATER_COKEY, WWErosionResultData.DOM_WATER_COMPNAME, WWErosionResultData.DOM_WATER_COMPONENT_AREA, WWErosionResultData.DOM_WATER_COMPONENT_AREA_PERCENT,
        WWErosionResultData.DOM_WIND_COKEY, WWErosionResultData.DOM_WIND_COMPNAME, WWErosionResultData.DOM_WIND_COMPONENT_AREA, WWErosionResultData.DOM_WIND_COMPONENT_AREA_PERCENT
    )));
  }

  public void correctForIrrigation(boolean useIrrigation) {
    if (useIrrigation) {
      if (resultData.wind_ifactor() <= 180) {
        for (int i = 0; i < i_factor_list.length; i++) {
          if (i_factor_list[i] == resultData.wind_ifactor()) {
            if (i > 0) {
              if (resultData.wind_ifactor() == 48) {
                resultData.wind_ifactor(21);
              } else {
                resultData.wind_ifactor(i_factor_list[i - 1]);
              }
              break;
            } else {
              break;
            }
          }
        }
      }
    }
  }

  public double correctForIrrigation(boolean adjust, double ifactor) {
    double ifact = ifactor;
    if (adjust) {
      if (ifact <= 180) {
        for (int i = 0; i < i_factor_list.length; i++) {
          if (i_factor_list[i] == ifact) {
            if (i > 0) {
              if (ifact == 48) {
                ifact = 21;
              } else {
                ifact = i_factor_list[i - 1];
              }
              break;
            } else {
              break;
            }
          }
        }
      }
    }
    return ifact;
  }

  protected double fixWEI(double wei) {
    if (EvalResult.testDefaultDouble(wei)) {
      return DEFAULT_WEI;
    }

    return wei;
  }

  protected void sumComponentData(boolean adjustForIrrig) throws ServiceException {
    String wind_dom_comp = "";
    String wind_dom_compname = "";
    double wind_dom_sand = 0.0;
    double wind_dom_clay = 0.0;
    double wind_ifact = 0.0;
    double wind_tfact = 0.0;
    double wind_sandtotal = 0.0;
    double wind_area = 0.0;
    double wind_area_pct = 0.0;

    String water_dom_comp = "";
    String water_dom_compname = "";
    double water_ifact = 0.0;
    double water_dom_sand = 0.0;
    double water_dom_clay = 0.0;
    double waterlsfactor = 0.0;
    double watertfactor = 0.0;
    double water_kffact = 0.0;
    double water_area = 0.0;
    double water_area_pct = 0.0;

    //Weighted averages variables
    double totalArea = 0.0;
    double weighted_water_ep = 0.0;
    double heighest_water_ep = 0.0;
    Component heighest_ep_comp = null;

    //Weighted averages for wind variables
    double totalWindArea = 0.0;
    double weighted_wind_ep = 0.0;
    double heighest_wind_ep = 0.0;
    Component heighest_wind_ep_comp = null;

    for (MapUnit mapUnit : map_units.values()) {
      if (!mapUnit.isExcluded()) {
        for (Component component : mapUnit.components().values()) {
          if (!component.isExcluded()) {

            //Calculate AoA Level Wind Values  
            if ((component.calculated_sandtotal_r() >= wind_sandtotal) && (component.calculated_sandtotal_r() > 0) && component.area_pct() >= minPercentage) { //Sand content is the main determining factor here..
              if (wind_sandtotal == component.calculated_sandtotal_r()) { // compare sizes of this component area to others since sandtotal matches greatest value found so far...
                if (component.calculated_area() > wind_area) {
                  wind_sandtotal = component.calculated_sandtotal_r();
                  wind_ifact = correctForIrrigation(adjustForIrrig, fixWEI(component.wei()));
                  wind_tfact = component.tfact();
                  wind_dom_comp = component.cokey();
                  wind_dom_compname = component.compname();
                  wind_dom_sand = component.calculated_sandtotal_r();
                  wind_dom_clay = component.calculated_claytotal_r();
                  wind_area = component.calculated_area();
                  wind_area_pct = component.area_pct();
                }
              } else { // Use this compoenent because its sandtotal is greater than the last...
                wind_sandtotal = component.calculated_sandtotal_r();
                wind_ifact = correctForIrrigation(adjustForIrrig, fixWEI(component.wei()));
                wind_tfact = component.tfact();
                wind_dom_comp = component.cokey();
                wind_dom_compname = component.compname();
                wind_dom_sand = component.calculated_sandtotal_r();
                wind_dom_clay = component.calculated_claytotal_r();
                wind_area = component.calculated_area();
                wind_area_pct = component.area_pct();
              }
            }

            //  Calculate AoA Level Water Values
            if ((component.calculated_kffact() >= water_kffact) && (component.calculated_kffact() > 0) && component.area_pct() >= minPercentage) {
              if (water_kffact == component.calculated_kffact()) {
                if (component.calculated_area() > water_area) {
                  water_kffact = component.calculated_kffact();
                  waterlsfactor = component.calculated_lsfact();
                  watertfactor = component.tfact();
                  water_ifact = correctForIrrigation(adjustForIrrig, fixWEI(component.wei()));
                  water_dom_comp = component.cokey();
                  water_dom_compname = component.compname();
                  water_dom_sand = component.calculated_sandtotal_r();
                  water_dom_clay = component.calculated_claytotal_r();
                  water_area = component.calculated_area();
                  water_area_pct = component.area_pct();
                }
              } else {
                water_kffact = component.calculated_kffact();
                waterlsfactor = component.calculated_lsfact();
                watertfactor = component.tfact();
                water_ifact = correctForIrrigation(adjustForIrrig, fixWEI(component.wei()));
                water_dom_comp = component.cokey();
                water_dom_compname = component.compname();
                water_dom_sand = component.calculated_sandtotal_r();
                water_dom_clay = component.calculated_claytotal_r();
                water_area = component.calculated_area();
                water_area_pct = component.area_pct();
              }
            }

            // SAVE wind_ep for each component not exclued
            if ((component.tfact() > 0)
                && (!EvalResult.testDefaultDouble(component.tfact()))) {
              component.calculated_wind_ep((resultData.cfactor() * correctForIrrigation(adjustForIrrig, fixWEI(component.wei()))) / component.tfact());
              if (null == heighest_wind_ep_comp) {
                heighest_wind_ep = component.calculated_wind_ep();
                heighest_wind_ep_comp = component;
              } else {
                if ((component.calculated_wind_ep() > heighest_wind_ep) && (component.area_pct() >= minPercentage)) {
                  heighest_wind_ep = component.calculated_wind_ep();
                  heighest_wind_ep_comp = component;
                } else {
                  if ((component.calculated_wind_ep() == heighest_wind_ep) && (component.area_pct() >= minPercentage) && (component.calculated_area() > heighest_wind_ep_comp.calculated_area())) {
                    heighest_wind_ep = component.calculated_wind_ep();
                    heighest_wind_ep_comp = component;
                  }
                }
              }
              if (component.area_pct() >= minPercentage) {
                weighted_wind_ep += (component.calculated_wind_ep() * component.calculated_area());
                totalWindArea += component.calculated_area();
              }
            }

            // SAVE water_ep for each component not excluded
            if ((component.tfact() > 0)
                && (!EvalResult.testDefaultDouble(component.tfact())
                && (!EvalResult.testDefaultDouble(component.calculated_kffact()))
                && (!EvalResult.testDefaultDouble(component.calculated_lsfact())))) {
              component.calculated_water_ep((component.calculated_kffact() * component.calculated_lsfact()) / component.tfact());
              if ((component.calculated_water_ep() > heighest_water_ep) && (component.area_pct() >= minPercentage)) {
                heighest_water_ep = component.calculated_water_ep();
                heighest_ep_comp = component;
              } else {
                if ((component.calculated_water_ep() == heighest_water_ep) && (component.area_pct() >= minPercentage) && (component.calculated_area() > heighest_ep_comp.calculated_area())) {
                  heighest_water_ep = component.calculated_water_ep();
                  heighest_ep_comp = component;
                }
              }
              if (component.area_pct() >= minPercentage) {
                weighted_water_ep += (component.calculated_water_ep() * component.calculated_area());
                totalArea += component.calculated_area();
              }
            }
          }
        }
      }
    }

    if (totalWindArea > 0) {
      resultData.weighted_avg_wind_ep(weighted_wind_ep / totalWindArea);
    }
    if (null != heighest_wind_ep_comp) {
      resultData.heighest_wind_ep(heighest_wind_ep);
      resultData.heighest_wind_ep_cokey(heighest_wind_ep_comp.cokey());
      resultData.heighest_wind_ep_compname(heighest_wind_ep_comp.compname());
      resultData.heighest_wind_ep_area_pct(heighest_wind_ep_comp.area_pct());
      resultData.heighest_wind_ep_area(heighest_wind_ep_comp.calculated_area());
      resultData.heighest_wind_ep_ifactor(correctForIrrigation(adjustForIrrig, fixWEI(heighest_wind_ep_comp.wei())));
      resultData.heighest_wind_ep_tfactor(heighest_wind_ep_comp.tfact());
    }

    if (totalArea > 0) {
      resultData.weighted_avg_water_ep(weighted_water_ep / totalArea);
    }
    if (null != heighest_ep_comp) {
      resultData.heighest_water_ep(heighest_water_ep);
      resultData.heighest_water_ep_cokey(heighest_ep_comp.cokey());
      resultData.heighest_water_ep_compname(heighest_ep_comp.compname());
      resultData.heighest_water_ep_area_pct(heighest_ep_comp.area_pct());
      resultData.heighest_water_ep_area(heighest_ep_comp.calculated_area());
      resultData.heighest_water_ep_lsfactor(heighest_ep_comp.calculated_lsfact());
      resultData.heighest_water_ep_tfactor(heighest_ep_comp.tfact());
      resultData.heighest_water_ep_kfactor(heighest_ep_comp.calculated_kffact());
    }

    resultData.dom_wind_cokey(wind_dom_comp);
    resultData.dom_wind_compname(wind_dom_compname);
    resultData.dom_wind_sand(wind_dom_sand);
    resultData.dom_wind_clay(wind_dom_clay);
    resultData.dom_wind_area_pct(wind_area_pct);
    resultData.dom_wind_area(wind_area);
    resultData.wind_ifactor(wind_ifact);
    resultData.wind_tfactor(wind_tfact);

    resultData.dom_water_cokey(water_dom_comp);
    resultData.dom_water_compname(water_dom_compname);
    resultData.dom_water_sand(water_dom_sand);
    resultData.dom_water_clay(water_dom_clay);
    resultData.dom_water_area_pct(water_area_pct);
    resultData.dom_water_area(water_area);
    resultData.water_lsfactor(waterlsfactor);
    resultData.water_tfactor(watertfactor);
    //resultData.water_ifactor(water_ifact);
    resultData.water_kfactor(water_kffact);

    if (resultData.dom_wind_compname().isEmpty()) {
      throw new ServiceException(" Could not find any wind critical components/horizons in the database meeting the selection criteria to calculate wind erosion potential.");
    }

    if (resultData.dom_water_compname().isEmpty()) {
      throw new ServiceException(" Could not find any water critical components/horizons in the database meeting the selection criteria to calculate water erosion potential.");
    }
  }

  public void outputSoilsData(boolean value) {
    this.outputSoils = value;
  }

  public final int getSRID() {
    return shape.getGeometry().getSRID();
  }

  public void getErosion(Callable<Connection> r2gis, boolean adjustForIrrig) throws Exception {
    getCFactor(r2gis);
    getSoilsData(adjustForIrrig);
  }

  protected void getSoilsData(boolean adjustForIrrig) throws ServiceException, GISObjectException {

    findIntersectedMapUnitsWithAreasAndShapes();
    if (soilsDb.findAllBasicComponentHorizonFragTextureData(map_units)) {
      //findAllTextureData();
      for (MapUnit mapunit : map_units.values()) {
        for (Component tComponent : mapunit.components().values()) {
          tComponent.applyHorizonFilter("WIND");
          tComponent.applyHorizonFilter("WATER");
          tComponent.updateCalculations();
        }
      }
      sumComponentData(adjustForIrrig);

      resultData.wind_ep(((resultData.cfactor()) * resultData.wind_ifactor()) / resultData.wind_tfactor());
      resultData.water_ep((resultData.water_kfactor() * resultData.water_lsfactor()) / resultData.water_tfactor());
    } else {
      throw new ServiceException("There was an error retrieving all of the wind and water related soil component data");
    }

//    if (findIntersectedMapUnitsWithAreas()) {
//      if (findAllComponentsApplyWindWaterFilters()) {
//        sumComponentData();
//        resultData.wind_ep(((resultData.cfactor()) * resultData.wind_ifactor()) / resultData.wind_tfactor());
//        resultData.water_ep((resultData.water_kfactor() * resultData.water_lsfactor()) / resultData.water_tfactor());
//      } else {
//        throw new ServiceException("There was an error retrieving all of the wind and water related soil component data");
//      }
//    } else {
//      throw new ServiceException("No soil mapunits were found that intersected with the shape provided");
//    }
  }

  public void setOutputColumns(ArrayList<String> outputList) {
    resultData.setOutputColumnOrdering(outputList);
  }

  public void toJSON(MDS_Result output) throws JSONException, IOException, ServiceException {
    resultData.setNonOutputColumns(null);

    //resultData.setOutputColumnOrdering(new ArrayList<>(Arrays.asList(AOA_AREA, CFACTOR, CFACTOR_PCT, CFACTOR_SOURCE,         )));
    resultData.toJSON(output);
    if (outputSoils) {
      for (MapUnit mapUnit : map_units.values()) {
        mapUnit.setOutputColumnOrdering(new ArrayList<String>(Arrays.asList(
            TableMapUnit.MUKEY, TableMapUnit.MUSYM, TableMapUnit.MUACRES,
            TableMapUnit.AREASYMBOL_NAME, TableMapUnitCalculations.AREA_NAME, TableMapUnitCalculations.AREA_PCT
        )));
        for (Component component : mapUnit.components().values()) {
          component.setOutputColumnOrdering(new ArrayList<String>(Arrays.asList(
              TableComponent.COKEY, TableComponent.COMPNAME, TableComponent.COMPPCT_R_NAME,
              TableComponentCalculations.AREA_PCT_NAME, TableComponentCalculations.COMP_AREA_NAME,
              TableComponentCalculations.LSFACT_NAME, TableComponentCalculations.KFFACT_NAME,
              TableComponent.TFACT_NAME, TableComponent.WEI_NAME, TableComponentCalculations.WATER_EP, TableComponentCalculations.WIND_EP
          )));
        }
      }

      JSONObject soilsData = toJSON(false);
      output.putResult("detailed_soils_data", soilsData, null, null);
    }
    output.putResult(EXCLUDED_LIST, getExcluded(), null, null);

    if (hasGeomChanged()) {
      output.putResult(CORRECTED_GEOMETRY, shapeWKT(), "If this section is present, the input geometry was invalid.  This service attempted to correct it, and was successful in creating a new geometry that could be utilized.  Please check this corrected geometry to be sure it represents what you originally intended.  If it does, please contact the source of your original geometry to have it corrected.  You may use this WKT to do so.", null);
    }
  }

  protected double lookupAlaskaCFactor(Callable<Connection> r2gis) throws Exception {
    double ret_val = EvalResult.getDefaultDouble();

    String query = "select TOP 1 name, (geometry::STGeomFromText('POINT(" + centroid[1] + " " + centroid[0]
        + ")',4326).STDistance(geom)) as distance  FROM r2gis.alaska_points order by distance;";

    try (Connection conn = r2gis.call(); Statement statement = conn.createStatement();) {
      try (ResultSet results = statement.executeQuery(query)) {
        if (results.next()) {
          ret_val = results.getDouble("name");
        } else {
          throw new ServiceException("Could not find an Alaska C_Factor for that input geography. No results returned from the Alaska Points query.");
        }
      }
    }
    return ret_val;
  }

  protected void getCFactor(Callable<Connection> r2gis) throws Exception {
    resultData.cfactor(EvalResult.getDefaultDouble());
    //  NOTE:  Remember:  the stored value of centroid places latitude in position 0 and longitude in position 1.
    GridCoordinates2D posGrid = null;
    synchronized (GIS_RASTER_LOCK) {
      posGrid = gisRaster.transformToRasterGrid(centroid[0], centroid[1], getSRID());
    }
    resultData.cfactor_source("Centroid");
    if (gisRaster.contains(posGrid)) {
      synchronized (GIS_RASTER_LOCK) {
        resultData.cfactor(gisRaster.pixelValue(posGrid));
        if (gisRaster.isNoData(resultData.cfactor() * 100.0)) {
//                        InterpolationMethodResult<Double> cfactorResult = new InterpolationMethodResult<>();
          InterpolationMethodResult<Double> cfactorResult
              = gisRaster.interpolatedPixelValue(shape, new InterpolationMethodDoubleMean(), true, 3);

          resultData.cfactor(cfactorResult.value);
          resultData.cfactor_source(cfactorResult.methodName);

          if (gisRaster.isNoData(resultData.cfactor())) {
            // Per requests from NRCS we will default to a zero value if no c_factor can be located.
            resultData.cfactor(0.0);
            resultData.cfactor_source("Default Zero");
          }
        }
      }
    } else {
      if (!soilsDb.isCONUS(centroid[0], centroid[1])) {
        resultData.cfactor(lookupAlaskaCFactor(r2gis));
      } else {
        // Per requests from NRCS we will default to a zero value if no c_factor can be located.
        resultData.cfactor(0.0);
        resultData.cfactor_source("Default Zero");
        //throw new ServiceException("The shape, for which you requested service, is not contained within the cfactor raster map (CONUS), nor contained within the State of Alaska.");
      }
    }
  }

}