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

import csip.ModelDataService;
import csip.api.server.ServiceException;
import csip.annotations.Resource;
import javax.ws.rs.Path;
import csip.annotations.*;
import static csip.annotations.ResourceType.EXECUTABLE;
import static csip.annotations.ResourceType.REFERENCE;
import csip.utils.ZipFiles;
import java.io.File;
import m.watershed.Resources;
import static m.watershed.Resources.*;
import static m.watershed.Utils.checkExist;
import static m.watershed.Utils.exec;
import org.codehaus.jettison.json.JSONArray;
import org.codehaus.jettison.json.JSONObject;

/**
 * slope
 *
 * @author HK (using OD's template)
 */
@Name("Extent to WEPP input data")
@Description("Example of extent (lat long) to WEPP input data")
@Path("m/merge_geometrie/1.0")
//@Resource(from = Resources.class)
@Resource(file = "${watershed.taudem.d8flowdir}", type = EXECUTABLE, id = D8)
@Resource(file = "${watershed.taudem.pitremove}", type = EXECUTABLE, id = PITREMOVE)
@Resource(file = "${watershed.taudem.dinfflowdir}", type = EXECUTABLE, id = DINF)
@Resource(file = "${watershed.taudem.aread8}", type = EXECUTABLE, id = AREAD8)
@Resource(file = "${watershed.taudem.gridnet}", type = EXECUTABLE, id = GRIDNET)
@Resource(file = "${watershed.taudem.peukerdouglas}", type = EXECUTABLE, id = PDOUG)
@Resource(file = "${watershed.taudem.areadinf}", type = EXECUTABLE, id = AREADINF)
@Resource(file = "${watershed.taudem.dropanalysis}", type = EXECUTABLE, id = DROP)
@Resource(file = "${watershed.taudem.streamnet}", type = EXECUTABLE, id = STREAMNET)
@Resource(file = "${watershed.taudem.d8hdisttostrm}", type = EXECUTABLE, id = D8STREAMDIST)
@Resource(file = "${watershed.taudem.moveoutletstostrm}", type = EXECUTABLE, id = MOVEOUTLET)
@Resource(file = "${watershed.taudem.threshold}", type = EXECUTABLE, id = THRESHOLD)

@Resource(file = LOCATION_PY + "/zonalstats.py", id = ZONALSTATS)
@Resource(file = LOCATION_PY + "/zonalstats_new.py", id = ZONALSTATS_NEW)
@Resource(file = LOCATION_PY + "/zonalstats_mean.py", id = ZONALSTATS_MEAN)
@Resource(file = LOCATION_PY + "/zonalstats_pts.py", id = ZONALSTATS_PTS)
@Resource(file = LOCATION_PY + "/readmm.py", id = READMM)
@Resource(file = LOCATION_PY + "/clean.py", id = CLEAN)
@Resource(file = LOCATION_PY + "/LS.py", id = LS)
@Resource(file = LOCATION_PY + "/slope_length.py", id = SLOPELENGTH)
@Resource(file = LOCATION_PY + "/cfactor.py", id = CFACTOR)
@Resource(file = LOCATION_PY + "/curve_number.py", id = CURVENUMBER)

@Resource(file = "${watershed.mpi}", type = REFERENCE, id = MPI)

@Resource(file = "${watershed.gdal}", type = REFERENCE, id = GDAL)
@Resource(file = "${watershed.gdaltranslate}", type = REFERENCE, id = GDALTRANSLATE)
@Resource(file = "${watershed.gdalraster}", type = REFERENCE, id = GDALRASTER)
@Resource(file = "${watershed.gdaldem}", type = REFERENCE, id = GDALDEM)
@Resource(file = "${watershed.gdalinfo}", type = REFERENCE, id = GDALINF)
@Resource(file = "${watershed.gdalpolygonize}", type = REFERENCE, id = GDALPOLYGONIZE)
@Resource(file = "${watershed.gdalcalc}", type = REFERENCE, id = GDALCALC)
@Resource(file = "${watershed.gdalmerge}", type = REFERENCE, id = GDALMERGE)

@Resource(file = "${watershed.ogr}", type = REFERENCE, id = OGR)
@Resource(file = "${watershed.ogrinfo}", type = REFERENCE, id = OGRINFO)

@Resource(file = "${watershed.python}", type = REFERENCE, id = PYTHON)
@Resource(file = "${watershed.python2}", type = REFERENCE, id = PYTHON2)
@Resource(file = "${watershed.python3}", type = REFERENCE, id = PYTHON3)

@Gzip
public class V1_0 extends ModelDataService {

    static final String JSONGEOM = "jsongeom";

    JSONObject boundary = null;

    boolean jsongeometry = false;

    double area = 0;
    int max_field = 0;

    @Override
    protected void doProcess() throws Exception {
        if (parameter().has(JSONGEOM)) {
            boundary = parameter().getJSONArray(JSONGEOM).getJSONObject(0);
            jsongeometry = true;
        }

        if (jsongeometry) {
            LOG.info("===>  Processing ");
            run(boundary, jsongeometry);
        } else {
            throw new ServiceException("No Extent or no GeoJSON Geometry input listed!");
        }
    }

    private void run(JSONObject boundary, boolean jsongeometry) throws Exception {

        int threads = Resources.mpiThreads(parameter());

        setProgress("Running Merging ...");

        JSONArray field_json = boundary.getJSONObject("geometries").getJSONArray("field");
        System.out.println(" # fields: " + field_json.length());

        for (int i = 0; i < field_json.length(); i++) {
            workspace().writeString("field_" + i + ".geojson", field_json.getJSONObject(i).getJSONArray("value").getJSONObject(0).toString());
            exec(resources(), LOG,
                    OGR,
                    "-explodecollections",
                    "-skipfailures",
                    "-f", "KML",
                    workspace().getFile("field_" + i + ".kml"),
                    workspace().getFile("field_" + i + ".geojson")
            );
            checkExist(workspace().getDir(), "field_" + i + ".kml");
        }

        exec(resources(), LOG,
                OGR,
                "-explodecollections",
                "-skipfailures",
                "-f", "KML",
                workspace().getFile("merged.kml"),
                workspace().getFile("field_0.kml"),
                "-nln", "field_0"
        );
        checkExist(workspace().getDir(), "merged.kml");

        JSONArray wsd_json = boundary.getJSONObject("geometries").getJSONArray("watersheds");
        System.out.println(" # wsds: " + wsd_json.length());

        for (int i = 0; i < wsd_json.length(); i++) {
            JSONObject wsd = wsd_json.getJSONObject(i).getJSONArray("value").getJSONObject(0);
//          //wsd.put("name", "wsd_" + i);
//          //wsd.append("name", "wsd_" + i);
            //workspace().writeString("wsd_" + i + ".geojson", wsd_json.getJSONObject(i).getJSONArray("value").getJSONObject(0).toString());
            workspace().writeString("wsd_" + i + ".geojson", wsd.toString());
            exec(resources(), LOG,
                    OGR,
                    "-explodecollections",
                    "-skipfailures",
                    "-f", "KML",
                    workspace().getFile("wsd_" + i + ".kml"),
                    workspace().getFile("wsd_" + i + ".geojson")
            );
            checkExist(workspace().getDir(), "wsd_" + i + ".kml");

//            if (new File(workspace().getDir(), "wsd_" + i + ".geojson").delete()) {
//                System.out.println("File deleted successfully");
//            } else {
//                throw new ServiceException("File not found: " + new File(workspace().getDir(), "wsd_" + i + ".geojson"));
//            }
//            exec(resources(), LOG,
//                    OGR,
//                    "-explodecollections",
//                    "-skipfailures",
//                    "-f", "GeoJSON",
//                    workspace().getFile("wsd_" + i + ".geojson"),
//                    workspace().getFile("wsd_" + i + ".kml"),
//                    "-nln", "wsd_" + i
//            );
//            checkExist(workspace().getDir(), "wsd_" + i + ".geojson");
            if (i < 1) {
//                exec(resources(), LOG,
//                        OGR,
//                        "-explodecollections",
//                        "-skipfailures",
//                        "-f", "GeoJSON",
//                        workspace().getFile("union_" + i + ".geojson"),
//                        workspace().getFile("field_0.geojson"),
//                        "-dialect", "sqlite",
//                        "-sql", "select ST_Buffer(ST_Union(field.geometry, wsd.geometry),0) from field_0 as field, 'wsd_0.geojson'.wsd_0 as wsd"
//                );
//                checkExist(workspace().getDir(), "union_" + i + ".geojson");
            } else {
                exec(resources(), LOG,
                        OGR,
                        "-explodecollections",
                        "-skipfailures",
                        "-f", "KML",
                        "-update",
                        "-append",
                        workspace().getFile("merged.kml"),
                        workspace().getFile("wsd_" + i + ".geojson"),
                        "-nln", "field_0"
                );
//                System.out.println(" " + i);
//                exec(resources(), LOG,
//                        OGR,
//                        "-explodecollections",
//                        "-skipfailures",
//                        "-f", "GeoJSON",
//                        workspace().getFile("union_" + i + ".geojson"),
//                        workspace().getFile("union_" + (i - 1) + ".geojson"),
//                        "-dialect", "sqlite",
//                        "-sql", "select ST_Buffer(ST_Union(field.geometry, wsd.geometry),0) from 'union_" + (i - 1) + ".geojson'.SELECT as field, 'wsd_" + i + ".geojson'.wsd_" + i + " as wsd"
//                );
//                checkExist(workspace().getDir(), "union_" + i + ".geojson");
            }
            max_field += 1;

        }

        exec(resources(), LOG,
                OGR,
                "-explodecollections",
                "-skipfailures",
                "-f", "GeoJSON",
                workspace().getFile("merged.geojson"),
                workspace().getFile("merged.kml"),
                "-nln", "field_0"
        );
        checkExist(workspace().getDir(), "merged.geojson");

//        exec(resources(), LOG,
//                OGR,
//                "-t_srs", "crs:84",
//                "-f", "GeoJSON",
//                "-update",
//                "-append",
//                "-dialect",
//                "sqlite",
//                "-sql",
//                "select ST_buffer(geometry,0) from 'merged.geojson'.field_0",
//                workspace().getFile("extent.geojson"),
//                workspace().getFile("merged.geojson")
//        );
//        checkExist(workspace().getDir(), "extent.geojson");
        setProgress("Finished Delineation.");
    }

    @Override
    protected void postProcess() throws Exception {
        results().put("area_km2", area);

        ZipFiles.zip(workspace().getFile("output.zip"),
                (Object) workspace().getFiles(
                        //"union_" + (max_field - 1) + ".geojson"
                        "*.geojson",
                        "*.kml"
                )
        );
        results().put(workspace().getFile("output.zip"));

        if (parameter().getBoolean("std_zip", false)) {
            ZipFiles.zip(workspace().getFile("std.zip"), (Object) workspace().getFiles("*.txt"));
            results().put(workspace().getFile("std.zip"));
        }
    }

    public static double haversine(double x) {
        return (1.0 - Math.cos(x)) / 2.0;
    }

    public double sphericalPolygonArea(double[] lat, double[] lon) {
        double lam1 = 0, lam2 = 0, beta1 = 0, beta2 = 0, cosB1 = 0, cosB2 = 0;
        double sum = 0;

        for (int j = 0; j < lat.length; j++) {
            int k = j + 1;
            if (j == 0) {
                lam1 = lon[j];
                beta1 = lat[j];
                lam2 = lon[j + 1];
                beta2 = lat[j + 1];
                cosB1 = Math.cos(beta1);
                cosB2 = Math.cos(beta2);
            } else {
                k = (j + 1) % lat.length;
                lam1 = lam2;
                beta1 = beta2;
                lam2 = lon[k];
                beta2 = lat[k];
                cosB1 = cosB2;
                cosB2 = Math.cos(beta2);
            }
            if (lam1 != lam2) {
                double hav = haversine(beta2 - beta1) + cosB1 * cosB2 * haversine(lam2 - lam1);
                double a = 2 * Math.asin(Math.sqrt(hav));
                double b = Math.PI / 2 - beta2;
                double c = Math.PI / 2 - beta1;
                double s = 0.5 * (a + b + c);
                double t = Math.tan(s / 2) * Math.tan((s - a) / 2) * Math.tan((s - b) / 2) * Math.tan((s - c) / 2);
                double excess = Math.abs(4 * Math.atan(Math.sqrt(Math.abs(t))));
                if (lam2 < lam1) {
                    excess = -excess;
                }

                sum += excess;
            }
        }
        return Math.abs(sum) * R * R;
    }

    static final double R = 6371.001; // kilometers

}