/*
 * 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.watershed.topography_index;

import csip.api.server.Executable;
import csip.ModelDataService;
import csip.annotations.Resource;
import static csip.annotations.ResourceType.*;
import javax.ws.rs.Path;
import csip.annotations.*;
import java.io.*;

/**
 * slope
 *
 * @author JK (using OD's template)
 */
@Name("topography_index")
@Description("Example of topography_index")
@Path("m/topography_index/1.0")
@Resource(file = "/python/zonalstats.py", type = FILE, id = "zonalstats")
@Resource(file = "gdal_calc.py", type = REFERENCE, id = "gdal_calc")
@Resource(file = "gdalwarp", type = REFERENCE, id = "gdalwarp")
@Resource(file = "gdal_translate", type = REFERENCE, id = "gdal_translate")
@Resource(file = "ogr2ogr", type = REFERENCE, id = "ogr")
//@Resource(file = "mpirun", type = REFERENCE, id = "mpirun")
@Resource(file = "python", type = REFERENCE, id = "python")
@Resource(file = "/bin/lin-amd64/pitremove", type = EXECUTABLE, id = "pitremove")
@Resource(file = "/bin/lin-amd64/dinfflowdir", type = EXECUTABLE, id = "dinfflowdir")
@Resource(file = "/bin/lin-amd64/areadinf", type = EXECUTABLE, id = "areadinf")

public class V1_0 extends ModelDataService {

    File result_file = null;

    // 2) watershed
    @Override
    protected void doProcess() throws Exception {

        String ip = request().getRemoteAddr();

        File ws = workspace().getDir();
        String wd = ws.toString();

        String file1 = parameter().getString("boundary");
        File file_1 = attachments().getFile(file1);
        String boundary_path = file_1.getPath();

        Executable e = resources().getExe("ogr");
        e.addArguments("-t_srs", "EPSG:5070", "-s_srs", "EPSG:4326", wd + "/boundary.shp", boundary_path);
        e.exec();

        String x_min = parameter().getString("l");
        String y_min = parameter().getString("b");
        String x_max = parameter().getString("r");
        String y_max = parameter().getString("t");

        e = resources().getExe("gdalwarp");
        e.addArguments("/mnt/csip-watershed/DEM.vrt", wd + "/DEM_clip.tif", "-te", x_min, y_min, x_max, y_max, "-multi");
        e.exec();

        //e = resources().getExe("mpirun");
        //File e2 = resources().getFile("pitremove");
        e = resources().getExe("pitremove");
        //e.addArguments("-np", "2", "-host", "localhost", e2.getAbsolutePath(), "-z", wd+ "/DEM_clip.tif", "-fel", wd + "/dem_fill.tif");        
        e.addArguments("-z", wd + "/DEM_clip.tif", "-fel", wd + "/dem_fill.tif");
        e.exec();

        //e = resources().getExe("mpirun");
        //e2 = resources().getFile("dinfflowdir");
        e = resources().getExe("dinfflowdir");
        //e.addArguments("-np", "2", "-host", "localhost", e2.getAbsolutePath(), "-fel", wd+ "/dem_fill.tif", "-slp", wd + "/slp.tif", "-ang", wd + "/ang.tif");        
        e.addArguments("-fel", wd + "/dem_fill.tif", "-slp", wd + "/slp.tif", "-ang", wd + "/ang.tif");
        e.exec();

        //e = resources().getExe("mpirun");
        //e2 = resources().getFile("areadinf");
        e = resources().getExe("areadinf");
        //e.addArguments("-np", "2", "-host", "localhost", e2.getAbsolutePath(), "-ang", wd+ "/ang.tif", "-sca", wd + "/sca.tif");        
        e.addArguments("-ang", wd + "/ang.tif", "-sca", wd + "/sca.tif");
        e.exec();

        e = resources().getExe("gdal_calc");
        e.addArguments("-A", wd + "/slp.tif", "--outfile", wd + "/slope.tif", "--overwrite", "--calc", "0.0001*(A==0)+A*(A>0)");
        e.exec();

        e = resources().getExe("gdal_translate");
        e.addArguments("-of", "GTiff", wd + "/sca.tif", wd + "/sca_2.tif", "-a_nodata", -32768);
        e.exec();

        e = resources().getExe("gdal_calc");
        e.addArguments("-A", wd + "/sca_2.tif", "--outfile", wd + "/acc.tif", "--overwrite", "--calc", "1*(A<=0)+(A+1)*(A>0)", "--NoDataValue=0");
        e.exec();

        // 900 is the area of the pixel of NHDPlusV2.0 DEM because the x, y size is each 30m. If difference DEM is used, please update the area of pixel.
        e = resources().getExe("gdal_calc");
        e.addArguments("-A", wd + "/acc.tif", "-B", wd + "/slope.tif", "--outfile", wd + "/TI.tif", "--overwrite", "--calc", "1*(A==0)+log(900*A/B)*(A>0)");
        e.exec();

        File f_1 = resources().getFile("zonalstats");
        e = resources().getExe("python");
        e.addArguments(f_1.getAbsolutePath(), wd + "/boundary.shp", wd + "/TI.tif");
        e.exec();

    }

    // 3) provide the temperature as a result.
    @Override
    protected void postProcess() throws Exception {

        File ws = workspace().getDir();
        results().put(new File(ws, "TI.tif"), "Topography index");
        results().put(new File(ws, "results.csv"), "result");

    }

}