V4_0.java [src/java/m/watershed/average_slope] 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 m.watershed.average_slope;

import csip.ModelDataService;
import csip.api.server.PayloadParameter;
import csip.api.server.ServiceException;
import csip.annotations.Description;
import csip.annotations.Name;
import csip.annotations.Resource;
import static csip.annotations.ResourceType.OUTPUT;
import java.io.File;
import javax.ws.rs.Path;
import m.watershed.Resources;
import static m.watershed.Resources.D8;
import static m.watershed.Resources.DINF;
import static m.watershed.Resources.GDAL;
import static m.watershed.Resources.LOCATION_MNT_DATA1;
import static m.watershed.Resources.LOCATION_MNT_CSIP_DEM;
import static m.watershed.Resources.MPI;
import static m.watershed.Resources.OGR;
import static m.watershed.Resources.PITREMOVE;
import static m.watershed.Resources.PYTHON;
import static m.watershed.Resources.ZONALSTATS_MEAN;
import static m.watershed.Utils.exec;
import static m.watershed.Utils.checkExist;
import org.apache.commons.io.FileUtils;
import org.codehaus.jettison.json.JSONObject;

/**
 * average_slope
 *
 * @author @author HK (using JK's template (using OD's template))
 */
@Name("average_slope")
@Description("Example of average_slope")
@Path("m/average_slope/4.0")
@Resource(from = Resources.class)
//@Resource(file = "*.tif *.csv *.kml *.kmz *.geojson *.asc *.shp *.txt *.dbf *.prj *.shx", type = OUTPUT)
@Resource(file = "*.csv", type = OUTPUT)
public class V4_0 extends ModelDataService {

    static final String BOUNDARY = "boundary";
    static final String RESULTID = "result_ID";
    static final String DINFIN = "DInf";
    static final String DEM_RES = "DEM_RES";
    static final String SERVICETHREADS = "Threads";
    static final String DEM_RES1 = "DEM_res1";
    static final String AOI_BUF = "AOI_buf";

    boolean shp_file = false;
    String shp_name = "";

    @Override
    protected void doProcess() throws Exception {
        PayloadParameter inputPayload = parameter();

        if (inputPayload.has(BOUNDARY) && inputPayload.has(RESULTID)) {

            String b = parameter().getString(BOUNDARY);

            if (b.endsWith("shp")) {
                shp_file = true;
                File input_file = attachments().getFile(b);
                shp_name = input_file.getName();
            } else {
                JSONObject boundary = parameter().getJSONArray(BOUNDARY).getJSONObject(0);
                FileUtils.writeStringToFile(workspace().getFile("boundary.geojson"), boundary.toString(), "UTF-8");
            }

            String ID = parameter().getString(RESULTID);

            // optional parameter
            boolean TauDEMInf = parameter().getBoolean(DINFIN, false);
            int resolution = parameter().getInt(DEM_RES, 10);
            int threads = parameter().getInt(SERVICETHREADS, 2);
            int resample = parameter().getInt(DEM_RES1, 10);
            int buffer = parameter().getInt(AOI_BUF, 25);

            run(ID, TauDEMInf, resolution, threads, resample, buffer);
        } else {
            throw new ServiceException("Geometry Input File not found and/or Result-ID not found !! ");
        }
    }

    /**
     *
     * @param boundary
     * @param ID
     * @param TauDEMInf
     * @param resolution
     * @throws Exception
     *
     */
    private void run(String ID, boolean TauDEMInf, int resolution, int threads, int resample, int buffer) throws Exception {
        checkExist(new File(LOCATION_MNT_DATA1), "DEM.vrt");

        if (!shp_file) {
            checkExist(workspace().getDir(), "boundary.geojson");
            exec(resources(), LOG,
                    OGR,
                    "-t_srs", "EPSG:5070",
                    "-f", "ESRI Shapefile",
                    workspace().getFile("boundaryPre.shp"),
                    workspace().getFile("boundary.geojson")
            );

        } else {
            checkExist(workspace().getDir(), shp_name);
            exec(resources(), LOG,
                    OGR,
                    "-t_srs", "EPSG:5070",
                    "-f", "ESRI Shapefile",
                    workspace().getFile("boundaryPre.shp"),
                    workspace().getFile(shp_name)
            );
        }
        checkExist(workspace().getDir(), "boundaryPre.shp");
        if (resolution > 10) {
            exec(resources(), LOG,
                    OGR,
                    "-update",
                    "-append",
                    "-dialect", "sqlite",
                    "-sql", "select ST_buffer(extent(geometry), " + buffer + ")," + ID + " from boundaryPre",
                    workspace().getFile("boundaryBuf.shp"),
                    workspace().getFile("boundaryPre.shp")
            );
            checkExist(workspace().getDir(), "boundaryBuf.shp");
            exec(resources(), LOG,
                    GDAL,
                    LOCATION_MNT_DATA1 + "/DEM.vrt",
                    workspace().getFile("DEMclip0.tif"),
                    "-crop_to_cutline",
                    "-cutline",
                    workspace().getFile("boundaryBuf.shp"),
                    "-multi",
                    "-of", "GTiff",
                    "-co", "TILED=YES",
                    "-co", "COMPRESS=LZW",
                    "-co", "BIGTIFF=YES",
                    "-ot", "Float32",
                    "-co", "NUM_THREADS=" + threads,
                    "-wo", "NUM_THREADS=" + threads,
                    "-wo", "CUTLINE_ALL_TOUCHED=TRUE",
                    "--config", "GDAL_CACHEMAX", "1000",
                    "-wm", "1000"
            );
            checkExist(workspace().getDir(), "DEMclip0.tif");
            exec(resources(), LOG,
                    GDAL,
                    "-t_srs", "EPSG:5070",
                    "-tr", "30", "30",
                    "-r", "cubicspline",
                    workspace().getFile("DEMclip0.tif"),
                    workspace().getFile("DEMclip.tif")
            );
            checkExist(workspace().getDir(), "DEMclip.tif");
        } else {
            exec(resources(), LOG,
                    OGR,
                    "-update",
                    "-append",
                    "-dialect", "sqlite",
                    "-sql", "select ST_buffer(extent(geometry), " + buffer + ")," + ID + " from boundaryPre",
                    workspace().getFile("boundaryBuf.shp"),
                    workspace().getFile("boundaryPre.shp")
            );
            checkExist(workspace().getDir(), "boundaryBuf.shp");
            exec(resources(), LOG,
                    GDAL,
                    LOCATION_MNT_CSIP_DEM + "/DEM.vrt",
                    workspace().getFile("DEMclip0.tif"),
                    "-crop_to_cutline",
                    "-cutline",
                    workspace().getFile("boundaryBuf.shp"),
                    "-multi",
                    "-of", "GTiff",
                    "-co", "TILED=YES",
                    "-co", "COMPRESS=LZW",
                    "-co", "BIGTIFF=YES",
                    "-ot", "Float32",
                    "-co", "NUM_THREADS=" + threads,
                    "-wo", "NUM_THREADS=" + threads,
                    "-wo", "CUTLINE_ALL_TOUCHED=TRUE",
                    "--config", "GDAL_CACHEMAX", "2000",
                    "-wm", "2000"
            );
            checkExist(workspace().getDir(), "DEMclip0.tif");

            if (resample < 1) {
                resample = 1;
            }
            if (resample > 90) {
                resample = 90;
            }

            exec(resources(), LOG,
                    GDAL,
                    "-t_srs", "EPSG:5070",
                    "-tr", "" + resample, "" + resample,
                    "-r", "cubicspline",
                    //"-r", "near",
                    workspace().getFile("DEMclip0.tif"),
                    workspace().getFile("DEMclip.tif"),
                    "-multi",
                    "-of", "GTiff",
                    "-co", "TILED=YES",
                    "-co", "COMPRESS=LZW",
                    "-co", "BIGTIFF=YES",
                    "-ot", "Float32",
                    "-co", "NUM_THREADS=" + threads,
                    "-wo", "NUM_THREADS=" + threads,
                    "-wo", "CUTLINE_ALL_TOUCHED=TRUE",
                    "--config", "GDAL_CACHEMAX", "2000",
                    "-wm", "2000"
            );
            checkExist(workspace().getDir(), "DEMclip.tif");
        }
        exec(resources(), LOG,
                MPI,
                "--allow-run-as-root",
                "--path", resources().getFile(PITREMOVE).getParent(),
                "-wdir", workspace().getDir(),
                "--oversubscribe",
                "-np", threads,
                resources().getFile(PITREMOVE).getName(),
                "-z", workspace().getFile("DEMclip.tif"),
                "-fel", workspace().getFile("demfill.tif")
        );
        checkExist(workspace().getDir(), "demfill.tif");

        if (!TauDEMInf) {
            exec(resources(), LOG,
                    MPI,
                    "--allow-run-as-root",
                    "--path", resources().getFile(D8).getParent(),
                    "-wdir", workspace().getDir(),
                    "--oversubscribe",
                    "-np", threads,
                    resources().getFile(D8).getName(),
                    "-fel", workspace().getFile("demfill.tif"),
                    "-sd8", workspace().getFile("slope.tif"),
                    "-p", workspace().getFile("p.tif")
            );
        } else {
            exec(resources(), LOG,
                    MPI,
                    "--allow-run-as-root",
                    "--path", resources().getFile(DINF).getParent(),
                    "-wdir", workspace().getDir(),
                    "--oversubscribe",
                    "-np", threads,
                    resources().getFile(DINF).getName(),
                    "-fel", workspace().getFile("demfill.tif"),
                    "-slp", workspace().getFile("slope.tif"),
                    "-ang", workspace().getFile("flow_direction_inf_pre.tif")
            );
        }
        checkExist(workspace().getDir(), "slope.tif");

        // zonal stats
        exec(resources(), LOG,
                PYTHON,
                resources().getFile(ZONALSTATS_MEAN),
                workspace().getFile("boundaryPre.shp"),
                workspace().getFile("slope.tif"),
                ID,
                "avg_slope_",
                "False"
        );
        checkExist(workspace().getDir(), "avg_slope_results.csv");
    }

    @Override
    protected void postProcess() throws Exception {
        //results().put(workspace().getFile("boundary.geojson"), "boundary");
        results().put(workspace().getFile("avg_slope_results.csv"), "results");
    }
}