V1_0.java [src/java/m/rhem/rhem01_runmodel] Revision: a8b6f0ed4f1ed07f29b3c6281689247b1fbc142c  Date: Thu Sep 29 11:17:34 MDT 2016
/*
 * 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 m.rhem.rhem01_runmodel;

import csip.Executable;
import csip.ModelDataService;
import static csip.ModelDataService.EXEC_FAILED;
import static csip.ModelDataService.EXEC_OK;
import csip.ServiceException;
import csip.annotations.Polling;
import csip.annotations.Resource;
import static csip.annotations.ResourceType.EXECUTABLE;
import static csip.annotations.ResourceType.FILE;
import static csip.annotations.ResourceType.OUTPUT;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;
import java.io.PrintWriter;
import java.sql.Connection;
import java.sql.ResultSet;
import java.sql.SQLException;
import java.sql.Statement;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.logging.Level;
import javax.ws.rs.Path;
import oms3.annotations.Description;
import oms3.annotations.Name;
import rhem.utils.DBResources;

/**
 * @version 1.0
 * @author Rumpal Sidhu
 */
@Name("RHEM-01:Run RHEM Model")
@Description("Run RHEM Model utilizing parameters including climate station, "
        + "surface soil texture class, slope percent and length, and vegetative "
        + "cover characteristics (bunchgrass, forbs/annual, shrub, and sodgrass "
        + "foliar cover; plant basal cover; rock cover; litter cover; and "
        + "cryptogam cover.")
@Path("m/rhem/runrhem/1.0")
@Polling(first = 10000, next = 2000)
@Resource(from = DBResources.class)
@Resource(file = "/bin/win-x86/rhem_v2.3.exe", wine = true, id = "rhem", type = EXECUTABLE)
@Resource(file = "*.par *.pre *.sum *.out", type = OUTPUT)
@Resource(file = "/bin/win-x86/storm_input_865_5517.pre", type = FILE)

public class V1_0 extends ModelDataService {

    private AoA aoa;
    private String parameterFileName;
    private String stormFileName;
    private String summaryFileName;

    @Override
    public void preProcess() throws ServiceException {
        try {
            int aoaId = getIntParam("AoAID");
            int rhemSiteId = getIntParam("rhem_site_id");
            String scenarioName = getStringParam("scenarioname");
            String scenarioDescription = getStringParam("scenariodescription");
            int unit = getIntParam("units");
            String stateId = getStringParam("stateid");
            String climatestationId = getStringParam("climatestationid");
            String soilTexture = getStringParam("soiltexture");
            double slopeLength = getIntParam("slopelength");
            String slopeShape = getStringParam("slopeshape");
            double slopeSteepness = getIntParam("slopesteepness");
            double bunchGgrassCanopyCover = getIntParam("bunchgrasscanopycover");
            double forbsCanopyCover = getIntParam("forbscanopycover");
            double shrubsCanopyCover = getIntParam("shrubscanopycover");
            double sodGrassCanopyCover = getIntParam("sodgrasscanopycover");
            double basalCover = getIntParam("basalcover");
            double rockCover = getIntParam("rockcover");
            double litterCover = getIntParam("littercover");
            double cryptogamsCover = getIntParam("cryptogamscover");
            this.aoa = new AoA(aoaId, rhemSiteId, scenarioName,
                    scenarioDescription, unit, stateId, climatestationId,
                    soilTexture, slopeLength, slopeShape, slopeSteepness,
                    bunchGgrassCanopyCover, forbsCanopyCover, shrubsCanopyCover,
                    sodGrassCanopyCover, basalCover, rockCover, litterCover,
                    cryptogamsCover);
        } catch (ServiceException ex) {
            LOG.log(Level.SEVERE, "RHEM-01: Error in processing the request JSON.", ex);
            throw new ServiceException("Error in processing the request JSON.", ex);
        }
    }

    @Override
    public void doProcess() throws ServiceException, FileNotFoundException, IOException {
        try {
            this.parameterFileName = aoa.scenarioName + ".par";
            this.stormFileName = "storm_input_865_5517.pre";
            this.summaryFileName = aoa.scenarioName + ".sum";

            String line = null;
            FileReader fileReader = new FileReader(this.stormFileName);
            try (BufferedReader bufferedReader
                    = new BufferedReader(fileReader)) {
                try (PrintWriter writer
                        = new PrintWriter(new File(getWorkspaceDir(), this.stormFileName))) {
                    while ((line = bufferedReader.readLine()) != null) {
                        writer.println(line);
//                    System.out.println(line);
                    }
                } catch (FileNotFoundException ex) {
                    LOG.log(Level.SEVERE, "RHEM-01: FileNotFoundException.", ex);
                    throw new ServiceException("Problem in writing .par file", ex);
                }
            } catch (Exception e) {
                LOG.log(Level.SEVERE, "RHEM-01: ", e);
                throw new ServiceException("Exception: ", e);
            }

            this.generateParamFile();
//        this.generateSummaryFile();
            System.out.println("run model = " + this.runModel());
        } catch (Exception e) {
            LOG.log(Level.SEVERE, "RHEM-01: ", e);
            throw new ServiceException("Exception: ", e);
        }
    }

    private void generateParamFile() throws ServiceException {

        // canopy cover for grass
        //double grasscanopycover = aoa.bunchGgrassCanopyCover + aoa.forbsCanopyCover + aoa.sodGrassCanopyCover;
        // TOTAL CANOPY COVER
        double totalcanopycover = aoa.bunchGgrassCanopyCover + aoa.forbsCanopyCover + aoa.shrubsCanopyCover + aoa.sodGrassCanopyCover;
        // TOTAL GROUND COVER
        double totalgroundcover = aoa.rockCover + aoa.basalCover + aoa.litterCover + aoa.cryptogamsCover;

        if (aoa.unit == 2) {
            aoa.slopeLength = aoa.slopeLength * 0.3048;
        }

        double clen = aoa.slopeLength * 2.5;
        String unit;

        if (aoa.unit == 1) {
            unit = "metric";
        } else {
            unit = "English";
        }

        String diams = null;
        String density = null;
        String fract = null;

        double g, dist, por;
        g = dist = por = 0.0;

        try (Connection connection = getResourceJDBC(DBResources.MSSQL_SSURGO);
                Statement statement = connection.createStatement();) {

            String query = "SELECT clay_diameter, silt_diameter, "
                    + "small_aggregates_diameter, large_aggregates_diameter, "
                    + "sand_diameter, clay_specific_gravity, "
                    + "silt_specific_gravity, small_aggregates_specific_gravity, "
                    + "large_aggregates_specific_gravity, sand_specific_gravity, "
                    + "mean_matric_potential, pore_size_distribution, mean_porosity, "
                    + "clay_fraction, silt_fraction, small_aggregates_fraction, "
                    + "large_aggregates_fraction, sand_fraction "
                    + "FROM rhem.soil_texture_table "
                    + "WHERE class_name ='" + aoa.soilTexture + "';";
            try (ResultSet resultSet = statement.executeQuery(query)) {
                while (resultSet.next()) {
                    double diam1 = resultSet.getDouble("clay_diameter");
                    double diam2 = resultSet.getDouble("silt_diameter");
                    double diam3 = resultSet.getDouble("small_aggregates_diameter");
                    double diam4 = resultSet.getDouble("large_aggregates_diameter");
                    double diam5 = resultSet.getDouble("sand_diameter");
                    diams = diam1 + "   " + diam2 + "   " + diam3 + "   " + diam4 + "   " + diam5;
                    double density1 = resultSet.getDouble("clay_specific_gravity");
                    double density2 = resultSet.getDouble("silt_specific_gravity");
                    double density3 = resultSet.getDouble("small_aggregates_specific_gravity");
                    double density4 = resultSet.getDouble("large_aggregates_specific_gravity");
                    double density5 = resultSet.getDouble("sand_specific_gravity");
                    density = density1 + "   " + density2 + "   " + density3 + "   " + density4 + "   " + density5;
                    g = resultSet.getDouble("mean_matric_potential");
                    dist = resultSet.getDouble("pore_size_distribution");
                    por = resultSet.getDouble("mean_porosity");
                    double fract1 = resultSet.getDouble("clay_fraction");
                    double fract2 = resultSet.getDouble("silt_fraction");
                    double fract3 = resultSet.getDouble("small_aggregates_fraction");
                    double fract4 = resultSet.getDouble("large_aggregates_fraction");
                    double fract5 = resultSet.getDouble("sand_fraction");
                    fract = fract1 + "   " + fract2 + "   " + fract3 + "   " + fract4 + "   " + fract5;
                }
            }

        } catch (ServiceException | SQLException se) {
            LOG.log(Level.SEVERE, "RHEM-01: SQLException.", se);
            throw new ServiceException("SQL problem.", se);
        }

        double ft = (-1 * 0.109) + (1.425 * aoa.litterCover)
                + (0.442 * aoa.rockCover) + (1.764 * (aoa.basalCover
                + aoa.cryptogamsCover)) + (2.068 * aoa.slopeSteepness);
        ft = Math.pow(10, ft);

        double chezy = Math.pow((8 * 9.8) / ft, 0.5);

        String sl = null;
        String sx = null;
        switch (aoa.slopeShape) {

            case "Uniform":
                sl = aoa.slopeSteepness + ", " + aoa.slopeSteepness;
                sx = 0.00 + ", " + 1.00;
                break;

            case "Convex":
                sl = 0.001 + ", " + (aoa.slopeSteepness * 2);
                sx = 0.00 + ", " + 1.00;
                break;

            case "Concave":
                sl = (aoa.slopeSteepness * 2) + ", " + 0.001;
                sx = 0.00 + ", " + 1.00;
                break;

            case "S-shaped":
                sl = 0.001 + ", " + (aoa.slopeSteepness * 2) + ", " + 0.001;
                sx = 0.00 + ", " + 0.50 + ", " + 1.00;
                break;
        }

        double Keb = 0;
        switch (aoa.soilTexture) {
            case "Sand":
                Keb = 24 * Math.exp(0.3483 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Loamy Sand":
                Keb = 10 * Math.exp(0.8755 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Sandy Loam":
                Keb = 5 * Math.exp(1.1632 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Loam":
                Keb = 2.5 * Math.exp(1.5686 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Silt Loam":
                Keb = 1.2 * Math.exp(2.0149 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Silt":
                Keb = 1.2 * Math.exp(2.0149 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Sandy Clay Loam":
                Keb = 0.80 * Math.exp(2.1691 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Clay Loam":
                Keb = 0.50 * Math.exp(2.3026 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Silty Clay Loam":
                Keb = 0.40 * Math.exp(2.1691 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Sandy Clay":
                Keb = 0.30 * Math.exp(2.1203 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Silty Clay":
                Keb = 0.25 * Math.exp(1.7918 * (aoa.basalCover + aoa.litterCover));
                break;
            case "Clay":
                Keb = 0.2 * Math.exp(1.3218 * (aoa.basalCover + aoa.litterCover));
                break;
        }
        double weightedKe = 0;

        // calculate weighted Ke and Kss values for the vegetation types that have non-zero values
        if (totalcanopycover != 0) {
            weightedKe = weightedKe + ((aoa.shrubsCanopyCover / totalcanopycover) * (Keb * 1.2));
            weightedKe = weightedKe + ((aoa.sodGrassCanopyCover / totalcanopycover) * (Keb * 0.8));
            weightedKe = weightedKe + ((aoa.bunchGgrassCanopyCover / totalcanopycover) * (Keb * 1.0));
            weightedKe = weightedKe + ((aoa.forbsCanopyCover / totalcanopycover) * (Keb * 1.0));
        } else {
            weightedKe = Keb;
        }
        double ke = weightedKe;

        // Kss variables
        double Kss_Seg_Bunch, Kss_Seg_Sod, Kss_Seg_Shrub, Kss_Seg_Shrub_0, Kss_Seg_Forbs;
        double Kss_Average, Kss_Final;

        // 1)
        //   a) CALCULATE KSS FOR EACH VEGETATION COMMUNITY USING TOTAL FOLIAR COVER
        //		A)   BUNCH GRASS
        if (totalgroundcover < 0.475) {
            Kss_Seg_Bunch = 4.154 + 2.5535 * aoa.slopeSteepness
                    - 2.547 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Bunch = Math.pow(10, Kss_Seg_Bunch);

            Kss_Seg_Sod = 4.2169 + 2.5535 * aoa.slopeSteepness
                    - 2.547 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Sod = Math.pow(10, Kss_Seg_Sod);

            Kss_Seg_Shrub = 4.2587 + 2.5535 * aoa.slopeSteepness - 2.547 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Shrub = Math.pow(10, Kss_Seg_Shrub);

            Kss_Seg_Forbs = 4.1106 + 2.5535 * aoa.slopeSteepness - 2.547 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Forbs = Math.pow(10, Kss_Seg_Forbs);

            Kss_Seg_Shrub_0 = 4.2587 + 2.5535 * aoa.slopeSteepness - 2.547 * totalgroundcover;
            Kss_Seg_Shrub_0 = Math.pow(10, Kss_Seg_Shrub_0);

        } else {
            Kss_Seg_Bunch = 3.1726975 + 2.5535 * aoa.slopeSteepness
                    - 0.4811 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Bunch = Math.pow(10, Kss_Seg_Bunch);

            Kss_Seg_Sod = 3.2355975 + 2.5535 * aoa.slopeSteepness
                    - 0.4811 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Sod = Math.pow(10, Kss_Seg_Sod);

            Kss_Seg_Shrub = 3.2773975 + 2.5535 * aoa.slopeSteepness - 0.4811 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Shrub = Math.pow(10, Kss_Seg_Shrub);

            Kss_Seg_Forbs = 3.1292975 + 2.5535 * aoa.slopeSteepness - 0.4811 * totalgroundcover - 0.7822 * totalcanopycover;
            Kss_Seg_Forbs = Math.pow(10, Kss_Seg_Forbs);

            Kss_Seg_Shrub_0 = 3.2773975 + 2.5535 * aoa.slopeSteepness - 0.4811 * totalgroundcover;
            Kss_Seg_Shrub_0 = Math.pow(10, Kss_Seg_Shrub_0);
        }
        if (totalcanopycover > 0 && totalcanopycover < 0.02) {

            Kss_Average = totalcanopycover / 0.02 * ((aoa.shrubsCanopyCover / totalcanopycover) * Kss_Seg_Shrub
                    + (aoa.sodGrassCanopyCover / totalcanopycover) * Kss_Seg_Sod
                    + (aoa.bunchGgrassCanopyCover / totalcanopycover) * Kss_Seg_Bunch
                    + (aoa.forbsCanopyCover / totalcanopycover) * Kss_Seg_Forbs)
                    + (0.02 - totalcanopycover) / 0.02 * Kss_Seg_Shrub_0;

        } else {

            Kss_Average = (aoa.shrubsCanopyCover / totalcanopycover) * Kss_Seg_Shrub
                    + (aoa.sodGrassCanopyCover / totalcanopycover) * Kss_Seg_Sod
                    + (aoa.bunchGgrassCanopyCover / totalcanopycover) * Kss_Seg_Bunch
                    + (aoa.forbsCanopyCover / totalcanopycover) * Kss_Seg_Forbs;

        }

        // 3) CALCULATE KSS USED FOR RHEM
        if (totalcanopycover == 0) {
            if (totalgroundcover < 0.475) {
                Kss_Final = 4.2587 + 2.5535 * aoa.slopeSteepness - 2.547 * totalgroundcover;
                Kss_Final = Math.pow(10, Kss_Final);
            } else {
                Kss_Final = 3.2773975 + 2.5535 * aoa.slopeSteepness - 0.4811 * totalgroundcover;
                Kss_Final = Math.pow(10, Kss_Final);
            }
        } else if (totalgroundcover < 0.475) {
            Kss_Final = totalgroundcover / 0.475 * Kss_Average + (0.475 - totalgroundcover) / 0.475 * Kss_Seg_Shrub;
        } else {
            Kss_Final = Kss_Average;
        }

        Kss_Final = (Kss_Final * 1.3) * 2.0;
        double kss = Kss_Final;

        SimpleDateFormat sdf = new SimpleDateFormat("MMMM dd, YYYY, hh:mm a");
        Date now = new Date();
        String today = sdf.format(now);

        try (PrintWriter writer = new PrintWriter(new File(getWorkspaceDir(), this.parameterFileName))) {
            writer.println("! Parameter file for scenario: " + aoa.scenarioName);
            writer.println("! Date built: " + today + " (Version 2.3)");
            writer.println("! Parameter units: DIAMS(mm), DENSITY(g/cc),TEMP(deg C)");
            writer.println("BEGIN GLOBAL");
            writer.println(String.format("        %-8s=   %s", "CLEN", clen));
            writer.println(String.format("        %-8s=   %s", "UNITS", unit));
            writer.println(String.format("        %-8s=   %s", "DIAMS", diams));
            writer.println(String.format("        %-8s=   %s", "DENSITY", density));
            writer.println(String.format("        %-8s=   %s", "TEMP", 40));
            writer.println(String.format("        %-8s=   %s", "NELE", 1));
            writer.println("BEGIN GLOBAL");
            writer.println("BEGIN PLANE");
            writer.println(String.format("        %-8s=   %s", "ID", 1));
            writer.println(String.format("        %-8s=   %s", "LEN", aoa.slopeLength));
            writer.println(String.format("        %-8s=   %s", "WIDTH", 1.0));
            writer.println(String.format("        %-8s=   %s", "CHEZY", chezy));
            writer.println(String.format("        %-8s=   %s", "RCHEZY", chezy));
            writer.println(String.format("        %-8s=   %s", "SL", sl));
            writer.println(String.format("        %-8s=   %s", "SX", sx));
            writer.println(String.format("        %-8s=   %s", "CV", 1.00));
            writer.println(String.format("        %-8s=   %s", "SAT", 0.25));
            writer.println(String.format("        %-8s=   %s", "PR", 1));
            writer.println(String.format("        %-8s=   %s", "KSS", kss));
            writer.println(String.format("        %-8s=   %s", "KOMEGA", 0.000007747));
            writer.println(String.format("        %-8s=   %s", "KCM", 0.00029936430000));
            writer.println(String.format("        %-8s=   %s", "CA", 1.0));
            writer.println(String.format("        %-8s=   %s", "IN", 0.0));
            writer.println(String.format("        %-8s=   %s", "KE", ke));
            writer.println(String.format("        %-8s=   %s", "G", g));
            writer.println(String.format("        %-8s=   %s", "DIST", dist));
            writer.println(String.format("        %-8s=   %s", "POR", por));
            writer.println(String.format("        %-8s=   %s", "ROCK", 0.00));
            writer.println(String.format("        %-8s=   %s", "SMAX", 1.0));
            writer.println(String.format("        %-8s=   %s", "ADF", 0.0));
            writer.println(String.format("        %-8s=   %s", "ALF", 0.8));
            writer.println(String.format("        %-8s=   %s", "BARE", 0.23));
            writer.println(String.format("        %-8s=   %s", "RSP", 1.0));
            writer.println(String.format("        %-8s=   %s", "SPACING", 1.0));
            writer.println(String.format("        %-8s=   %s", "FRACT", fract));
            writer.println("END PLANE");
        } catch (FileNotFoundException ex) {
            LOG.log(Level.SEVERE, "RHEM-01: FileNotFoundException.", ex);
            throw new ServiceException("Problem in writing .par file", ex);
        }
    }

    private void generateSummaryFile() {

    }

    private String runModel() throws ServiceException {
        int run = -1;
        try {
            Executable rh = getResourceExe("rhem");
            rh.setArguments(this.parameterFileName, this.stormFileName,
                    this.summaryFileName, aoa.scenarioName, 0, 2, "y", "y", "n", "n", "y");
            run = rh.exec();

//        } catch (ServiceException | IOException se) {
        } catch (Exception se) {
            LOG.log(Level.SEVERE, "RHEM-01: ", se);
            throw new ServiceException("Problem in running the model.", se);
        }
        return run == 0 ? EXEC_OK : EXEC_FAILED;
    }

    @Override
    public void postProcess() throws Exception {
        try {
            putReport(new File(getWorkspaceDir(), this.parameterFileName), "Parameter input file");
            putReport(new File(getWorkspaceDir(), this.stormFileName), "Storm input file");
            putReport(new File(getWorkspaceDir(), this.summaryFileName), "Summary file");
        } catch (Exception e) {
            LOG.log(Level.SEVERE, "RHEM-01: ", e);
            throw new Exception("", e);
        }
    }

    static class AoA {

        protected int aoaId;
        protected int rhemSiteId;
        protected String scenarioName;
        protected String scenarioDescription;
        protected int unit;
        protected String stateId;
        protected String climatestationId;
        protected String soilTexture;
        protected double slopeLength;
        protected String slopeShape;
        protected double slopeSteepness;
        protected double moisturecontent;
        protected double bunchGgrassCanopyCover;
        protected double forbsCanopyCover;
        protected double shrubsCanopyCover;
        protected double sodGrassCanopyCover;
        protected double basalCover;
        protected double rockCover;
        protected double litterCover;
        protected double cryptogamsCover;

        public AoA(int aoaId, int rhemSiteId, String scenarioName,
                String scenarioDescription, int unit, String stateId,
                String climatestationId, String soilTexture, double slopeLength,
                String slopeShape, double slopeSteepness,
                double bunchGgrassCanopyCover, double forbsCanopyCover,
                double shrubsCanopyCover, double sodGrassCanopyCover, double basalCover,
                double rockCover, double litterCover, double cryptogamsCover) {
            this.aoaId = aoaId;
            this.rhemSiteId = rhemSiteId;
            this.scenarioName = scenarioName;
            this.scenarioDescription = scenarioDescription;
            this.unit = unit;
            this.stateId = stateId;
            this.climatestationId = climatestationId;
            this.soilTexture = soilTexture;
            this.slopeLength = slopeLength;
            this.slopeShape = slopeShape;
            this.slopeSteepness = slopeSteepness / 100.0;
            this.bunchGgrassCanopyCover = bunchGgrassCanopyCover / 100.0;
            this.forbsCanopyCover = forbsCanopyCover / 100.0;
            this.shrubsCanopyCover = shrubsCanopyCover / 100.0;
            this.sodGrassCanopyCover = sodGrassCanopyCover / 100.0;
            this.basalCover = basalCover / 100.0;
            this.rockCover = rockCover / 100.0;
            this.litterCover = litterCover / 100.0;
            this.cryptogamsCover = cryptogamsCover / 100.0;
            moisturecontent = moisturecontent / 100.0;
        }
    }
}