V1_0.java [src/java/m/watershed/average_aspect] 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_aspect;

import csip.Check;
import csip.ModelDataService;
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.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.util.HashMap;
import java.util.Map;
import javax.ws.rs.Path;
import m.watershed.Resources;
import static m.watershed.Resources.GDAL;
import static m.watershed.Resources.GDALDEM;
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.Resources.ZONALSTATS_MAX;
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_aspect
 *
 * @author @author HK (using JK's template (using OD's template))
 */
@Name("average_aspect")
@Description("Example of average_aspect")
@Path("m/average_aspect/1.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 V1_0 extends ModelDataService {

    static final String BOUNDARY = "boundary";
    static final String RESULTID = "result_ID";
    static final String DEM_RES = "DEM_RES";
    static final String FILL = "FILL";
    static final String SERVICETHREADS = "Threads";

    boolean shp_file = false;
    String shp_name = "";
    Integer resolution = 10;

    Map<String, String> GDALslope = new HashMap<>();
    Map<String, String> GDALaspect = new HashMap<>();
    Map<String, String> DEMelevation = new HashMap<>();
    Map<String, String> maxDEMelevation = new HashMap<>();

    @Override
    protected void doProcess() throws Exception {

        resolution = parameter().getInt(DEM_RES, 10);
        new Check().positive().forNumber(resolution);

        if (parameter().has(BOUNDARY) && parameter().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 filling = parameter().getBoolean(FILL, false);
            int threads = parameter().getInt(SERVICETHREADS, 2);

            run(ID, filling, resolution, threads);

        } else {
            throw new ServiceException("Geometry Input File not found and/or Result-ID not found !! ");
        }
    }

    /**
     *
     * @param boundary
     * @param ID
     * @param resolution
     * @throws Exception
     *
     */
    private void run(String ID, boolean filling, int resolution, int threads) throws Exception {

        String dem30 = LOCATION_MNT_DATA1 + "/DEM.vrt";
        String dem10 = LOCATION_MNT_CSIP_DEM + "/DEM.vrt";
        String dem3 = LOCATION_MNT_CSIP_DEM + "/DEM_1_9_arc/DEM_1_9.vrt";
        String dem1 = LOCATION_MNT_CSIP_DEM + "/DEM_1m/DEM_1m.vrt";

//        String selected_dem = dem10;
//        Integer pixel_width = 10;
//        String buffer = "2100";
//
//        if (resolution < 10) {
//            selected_dem = dem3;
//            pixel_width = 3;
//        }
//        if (resolution == 1) {
//            selected_dem = dem1;
//            pixel_width = 1;
//        }
//        if (resolution > 10) {
//            selected_dem = dem30;
//            pixel_width = 30;
//        }
//        
//                setProgress("Clip DEM...");
//        exec(resources(), LOG,
//                GDAL,
//                selected_dem,
//                workspace().getFile("DEM_clip0.tif"),
//                "-crop_to_cutline",
//                "-cutline",
//                workspace().getFile("indexFD.geojson"),
//                "-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", "2048",
//                "--config",
//                "GDALWARP_IGNORE_BAD_CUTLINE", "YES",
//                "-wm", "2048"
//        );
//        checkExist(workspace().getDir(), "DEM_clip0.tif");
//
//        setProgress("Reproject DEM...");
//        exec(resources(), LOG,
//                GDAL,
//                "-ot", "Float32",
//                "-t_srs", "+proj=utm +zone=" + zone + " +" + hem + " +datum=WGS84 +units=m +no_defs",
//                "-r", "cubicspline",
//                "-tr", pixel_width + "", pixel_width + "",
//                "-tap",
//                "-of", "GTiff",
//                "-co", "COMPRESS=DEFLATE",
//                "-co", "PREDICTOR=1",
//                "-co", "ZLEVEL=6",
//                "-co", "TILED=YES",
//                "-co", "NUM_THREADS=" + threads + "",
//                "-wo", "NUM_THREADS=" + threads + "",
//                "-wo", "OPTIMIZE_SIZE=TRUE",
//                "-wm", "2048",
//                workspace().getFile("DEM_clip0.tif"),
//                workspace().getFile("DEM_clip.tif")
//        );
//
//        checkExist(workspace().getDir(), "DEM_clip.tif");
        if (!shp_file) {
            checkExist(workspace().getDir(), "boundary.geojson");
            exec(resources(), LOG,
                    OGR,
                    "-t_srs", "EPSG:5070",
                    "-f", "ESRI Shapefile",
                    workspace().getFile("boundary_pre.shp"),
                    workspace().getFile("boundary.geojson")
            );
        } else {
            checkExist(workspace().getDir(), shp_name);
            exec(resources(), LOG,
                    OGR,
                    "-t_srs", "EPSG:5070",
                    "-f", "ESRI Shapefile",
                    workspace().getFile("boundary_pre.shp"),
                    workspace().getFile(shp_name)
            );
        }
        checkExist(workspace().getDir(), "boundary_pre.shp");

        if (resolution > 10) {
            exec(resources(), LOG,
                    OGR,
                    "-update",
                    "-append",
                    "-dialect", "sqlite",
                    "-sql", "select ST_buffer(extent(geometry), " + (resolution + 5) + ")," + ID + " from boundary_pre",
                    workspace().getFile("boundary.shp"),
                    workspace().getFile("boundary_pre.shp")
            );
            checkExist(workspace().getDir(), "boundary.shp");

            exec(resources(), LOG,
                    GDAL,
                    LOCATION_MNT_DATA1 + "/DEM.vrt",
                    workspace().getFile("DEM_clip0.tif"),
                    "-crop_to_cutline",
                    "-cutline",
                    workspace().getFile("boundary.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(), "DEM_clip0.tif");

            exec(resources(), LOG,
                    GDAL,
                    "-t_srs", "EPSG:5070",
                    "-tr", "30", "30",
                    "-r", "cubicspline",
                    workspace().getFile("DEM_clip0.tif"),
                    workspace().getFile("DEM_clip.tif")
            );
            checkExist(workspace().getDir(), "DEM_clip.tif");
        } else {
            exec(resources(), LOG,
                    OGR,
                    "-update",
                    "-append",
                    "-dialect", "sqlite",
                    "-sql", "select ST_buffer(extent(geometry), 1)," + ID + " from boundary_pre",
                    workspace().getFile("boundary.shp"),
                    workspace().getFile("boundary_pre.shp")
            );
            checkExist(workspace().getDir(), "boundary.shp");

            exec(resources(), LOG,
                    GDAL,
                    LOCATION_MNT_CSIP_DEM + "/DEM.vrt",
                    workspace().getFile("DEM_clip0.tif"),
                    "-crop_to_cutline",
                    "-cutline",
                    workspace().getFile("boundary.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(), "DEM_clip0.tif");

            exec(resources(), LOG,
                    GDAL,
                    "-t_srs", "EPSG:5070",
                    "-tr", "10", "10",
                    "-r", "cubicspline",
                    workspace().getFile("DEM_clip0.tif"),
                    workspace().getFile("DEM_clip.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(), "DEM_clip.tif");
        }

        String dem_usage = "DEM_clip.tif";
        if (filling) {
            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("DEM_clip.tif"),
                    "-fel", workspace().getFile("dem_filled.tif")
            );

            checkExist(workspace().getDir(), "dem_filled.tif");
            dem_usage = "dem_filled.tif";
        }

        exec(resources(), LOG,
                GDALDEM,
                "slope",
                workspace().getFile(dem_usage),
                workspace().getFile("GDALslope_deg.tif"),
                //"-p",
                "-of",
                "GTiff"//,
        //"-compute_edges"
        );
        exec(resources(), LOG,
                PYTHON,
                resources().getFile(ZONALSTATS_MEAN),
                workspace().getFile("boundary_pre.shp"),
                workspace().getFile("GDALslope_deg.tif"),
                ID,
                "avg_GDAL_slope_deg_",
                "False"
        );
        exec(resources(), LOG,
                GDALDEM,
                "aspect",
                workspace().getFile(dem_usage),
                workspace().getFile("GDALaspect_deg.tif"),
                //"-alg", "ZevenbergenThorne",
                "-zero_for_flat",
                "-of",
                "GTiff"//,
        //"-compute_edges"
        );
        exec(resources(), LOG,
                PYTHON,
                resources().getFile(ZONALSTATS_MEAN),
                workspace().getFile("boundary_pre.shp"),
                workspace().getFile("GDALaspect_deg.tif"),
                ID,
                "avg_GDAL_aspect_deg_",
                "False"
        );

        exec(resources(), LOG,
                PYTHON,
                resources().getFile(ZONALSTATS_MEAN),
                workspace().getFile("boundary_pre.shp"),
                workspace().getFile(dem_usage),
                ID,
                "avg_DEM_elevation_m_",
                "False"
        );

        exec(resources(), LOG,
                PYTHON,
                resources().getFile(ZONALSTATS_MAX),
                workspace().getFile("boundary_pre.shp"),
                workspace().getFile(dem_usage),
                ID,
                "max_DEM_elevation_m_",
                "False"
        );

        checkExist(workspace().getDir(),
                "avg_GDAL_aspect_deg_results.csv",
                "avg_DEM_elevation_m_results.csv",
                "max_DEM_elevation_m_results.csv",
                "avg_GDAL_slope_deg_results.csv");

        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_GDAL_aspect_deg_results.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                if (c > 0) {
                    String[] infos = line.split(",");
                    GDALaspect.put(infos[0], infos[1]);
                }
                c++;
            }
        }

        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_GDAL_slope_deg_results.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                if (c > 0) {
                    String[] infos = line.split(",");
                    GDALslope.put(infos[0], infos[1]);
                }
                c++;
            }
        }

        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_DEM_elevation_m_results.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                if (c > 0) {
                    String[] infos = line.split(",");
                    DEMelevation.put(infos[0], infos[1]);
                }
                c++;
            }
        }

        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("max_DEM_elevation_m_results.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                if (c > 0) {
                    String[] infos = line.split(",");
                    maxDEMelevation.put(infos[0], infos[1]);
                }
                c++;
            }
        }
    }

    @Override
    protected void postProcess() throws Exception {
        //results().put(workspace().getFile("boundary.geojson"), "boundary");
        results().put(workspace().getFile("avg_GDAL_aspect_deg_results.csv"), "results");
        results().put(workspace().getFile("avg_GDAL_slope_deg_results.csv"), "results");
        for (Map.Entry<String, String> entry : GDALslope.entrySet()) {
            String key = (String) entry.getKey();
            results().put("#ID" + key, key);
            results().put("avg_GDAL_slope_deg_" + key, entry.getValue());
            results().put("avg_GDAL_aspect_deg_" + key, GDALaspect.get(key));
            results().put("avg_DEM_elevation_m_" + key, DEMelevation.get(key));
            results().put("max_DEM_elevation_m_" + key, maxDEMelevation.get(key));
        }
    }
}