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

import csip.Check;
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.*;
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.JSONObject;

/**
 * slope
 *
 * @author HK (using OD's template)
 */
@Name("checking if any DEM cover is available")
@Description("Example of extent (lat long) to WEPP input data")
@Path("m/checkDEM_cover/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 = LOCATION_PY + "/grid_info.py", id = GRID_INFO)

@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.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 EXTENT = "extent";
    static final String JSONGEOM = "jsongeom";
    static final String RESOLUTION = "resolution";

    JSONObject boundary = null;

    boolean jsongeometry = false;
    boolean extent = false;

    boolean m1_check = false;
    boolean m3_check = false;
    boolean output = false;

    double area = 0;
    String latlong = "";
    Integer resolution = 10;
    double grid_cover_threshold = -1.0;
    double grid_value1 = -2.0;
    double grid_value3 = -2.0;

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

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

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

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

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

        setProgress("Running the DEM cover check ...");

        if (!jsongeometry && extent) {
            String s = "Name,wkt,\n"
                    + "Box,\"POLYGON((" + latlong + "))\"\n";
            workspace().writeString("field_extent.csv", s);

            String[] all_latlong = latlong.split(",");
            double[] llat_rad = new double[all_latlong.length - 1];
            double[] llong_rad = new double[all_latlong.length - 1];

            double d = Math.PI / 180;
            for (int i = 0; i < all_latlong.length - 1; i++) {
                String[] spl_string = all_latlong[i].split("\\s+");
                llat_rad[i] = Double.parseDouble(spl_string[1]) * d;
                llong_rad[i] = Double.parseDouble(spl_string[0]) * d;
            }

            area = sphericalPolygonArea(llat_rad, llong_rad);

            StringBuilder b = new StringBuilder();
            b.append("<OGRVRTDataSource>\n");
            b.append("    <OGRVRTLayer name=\"field_extent\">\n");
            b.append("        <SrcDataSource>field_extent.csv</SrcDataSource>\n");
            b.append("        <GeometryType>wkbPolygon</GeometryType>\n");
            b.append("        <LayerSRS>WGS84</LayerSRS>\n");
            b.append("        <Field name=\"Name\" src=\"Name\" />\n");
            b.append("        <GeometryField encoding=\"WKT\" field=\"wkt\" />\n");
            b.append("    </OGRVRTLayer>\n");
            b.append("</OGRVRTDataSource>\n");
            workspace().writeString("field_extent.vrt", b.toString());

            exec(resources(), LOG,
                    OGR,
                    "-f", "KML",
                    workspace().getFile("field_extent.kml"),
                    workspace().getFile("field_extent.vrt")
            );
            checkExist(workspace().getDir(), "field_extent.kml");

            exec(resources(), LOG,
                    OGR,
                    "-f", "ESRI Shapefile",
                    workspace().getFile("field_extent.shp"),
                    workspace().getFile("field_extent.kml")
            );
        } else {
            workspace().writeString("boundary.geojson", boundary.toString());
            exec(resources(), LOG,
                    OGR,
                    "-t_srs", "EPSG:4326",
                    "-f", "ESRI Shapefile",
                    workspace().getFile("field_extent.shp"),
                    workspace().getFile("boundary.geojson")
            );
        }
        checkExist(workspace().getDir(), "field_extent.shp");

        setProgress("Get field extend ...");
        exec(resources(), LOG,
                OGR,
                "-f", "CSV",
                "-dialect", "sqlite",
                "-sql", "select AsText(St_Convexhull(geometry)) AS geom, * from field_extent",
                workspace().getFile("field_boundary.csv"),
                workspace().getFile("field_extent.shp")
        );
        checkExist(workspace().getDir(), "field_boundary.csv");

        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("field_boundary.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                if (c == 1) {
                    String[] infos = line.split("\\)\\)\",");
                    String[] coord_array = infos[0].replaceAll("\"POLYGON\\(\\(", "").split(",\\s+");

                    double[] llat_rad = new double[coord_array.length - 1];
                    double[] llong_rad = new double[coord_array.length - 1];
                    for (int i = 0; i < coord_array.length - 1; i++) {
                        llat_rad[i] = Double.parseDouble((coord_array[i].split("\\s+")[1])) * Math.PI / 180;
                        llong_rad[i] = Double.parseDouble((coord_array[i].split("\\s+")[0])) * Math.PI / 180;
                    }
                    area = sphericalPolygonArea(llat_rad, llong_rad);
                }
                c++;
            }
        }

        setProgress("Get centroid ...");
        exec(resources(), LOG,
                OGR,
                "-f", "CSV",
                "-dialect", "sqlite",
                "-sql", "select AsText(ST_Centroid(geometry)) AS geom, * from field_extent",
                workspace().getFile("field_centroid.csv"),
                workspace().getFile("field_extent.shp")
        );
        checkExist(workspace().getDir(), "field_centroid.csv");

        String llat = "";
        String llong = "";
        try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("field_centroid.csv")))) {
            int c = 0;
            String line;
            while ((line = br.readLine()) != null) {
                // use comma as separator
                String[] infos = line.split(",");
                if (c == 1) {
                    llat = infos[0].split("\\s+")[1].replaceAll("\\)", "");
                    llong = infos[0].split("\\s+")[0].replaceAll("POINT\\(", "");
                }
                c++;
            }
        }

        int zone = 16;
        String hem = "north";
        if (!"".equals(llat)) {
            zone = (int) Math.floor(Double.parseDouble(llong) / 6 + 31);
            hem = ((int) Math.floor(Double.parseDouble(llat)) > 0 ? "north" : "south");
        }

        setProgress("Reproject ...");
        exec(resources(), LOG,
                OGR,
                "-t_srs", "+proj=utm +zone=" + zone + " +" + hem + " +datum=WGS84 +units=m +no_defs",
                "-f", "ESRI Shapefile",
                workspace().getFile("field_extent_utm.shp"),
                workspace().getFile("field_extent.shp")
        );
        checkExist(workspace().getDir(), "field_extent_utm.shp");

        setProgress("Buffering ...");
        exec(resources(), LOG,
                OGR,
                "-t_srs", "crs:84",
                "-f", "GeoJSON",
                "-update",
                "-append",
                "-dialect",
                "sqlite",
                "-sql",
                "select ST_buffer(extent(geometry), " + 10 + ") from field_extent_utm",
                workspace().getFile("extent.geojson"),
                workspace().getFile("field_extent_utm.shp")
        );
        checkExist(workspace().getDir(), "extent.geojson");

        setProgress("Update field extend ...");
        exec(resources(), LOG,
                OGRINFO,
                workspace().getFile("field_extent_utm.shp"),
                "-sql", "ALTER TABLE field_extent_utm ADD COLUMN IDX integer(4)"
        );

        exec(resources(), LOG,
                OGRINFO,
                workspace().getFile("field_extent_utm.shp"),
                "-dialect", "SQLite", "-sql", "UPDATE field_extent_utm SET IDX=1"
        );

        setProgress("Create Field Buffer...");

        exec(resources(), LOG,
                OGR,
                "-update",
                "-append",
                "-dialect",
                "sqlite",
                "-sql",
                "select extent(ST_buffer(geometry, " + "10" + ")) from field_extent_utm",
                //              "select ST_buffer(extent(geometry), " + buffer + ") from field_extent_utm",
                workspace().getFile("index.shp"),
                workspace().getFile("field_extent_utm.shp")
        );
        checkExist(workspace().getDir(), "index.shp");

        setProgress("Get extend ...");
        exec(resources(), LOG,
                OGR,
                "-f", "CSV",
                "-dialect", "sqlite",
                "-sql", "select AsText(ST_envelope(geometry)) as geom from 'index.shp'.index",
                workspace().getFile("extent.csv"),
                workspace().getFile("index.shp")
        );

        exec(resources(), LOG,
                OGR,
                "-f", "GeoJSON",
                "-t_srs", "crs:84",
                workspace().getFile("indexFD.geojson"),
                workspace().getFile("index.shp")
        );
        checkExist(workspace().getDir(), "indexFD.geojson");

        String dem3 = LOCATION_MNT_CSIP_DEM + "/DEM_1_9_arc/DEM_1_9.vrt";
        String dem1 = LOCATION_MNT_CSIP_DEM + "/DEM_1m/1m.vrt";

        String selected_dem = dem1;
        Integer pixel_width = 1;
        output = parameter().getBoolean("out_zip", false);

        setProgress("Checking 1m DEM ...");

        try {
            exec(resources(), LOG,
                    GDAL,
                    selected_dem,
                    workspace().getFile("DEM_clip1.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"
            );

            if (new File(workspace().getDir(), "DEM_clip1.tif").exists()) {
                exec(resources(), LOG,
                        PYTHON,
                        resources().getFile(GRID_INFO),
                        workspace().getFile("DEM_clip1.tif"),
                        workspace().getFile("pixel_test1.txt")
                );
                if (new File(workspace().getDir(), "pixel_test1.txt").exists()) {
                    try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("pixel_test1.txt")))) {
                        int c = 0;
                        String line;
                        while ((line = br.readLine()) != null) {
                            if (c == 0) {
                                grid_value1 = Double.parseDouble(line);
                                break;
                            }
                        }
                    }
                }
                if (grid_value1 > grid_cover_threshold) {
                    m1_check = true;
                }
            }

        } catch (FileNotFoundException e) {
        } catch (IOException e) {
        } finally {

            setProgress("1m DEM! ...");

        }

        if (m1_check && output) {
            setProgress("still checking 1m 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_clip1.tif"),
                    workspace().getFile("DEM_clip1re.tif")
            );

            checkExist(workspace().getDir(), "DEM_clip1re.tif");

            setProgress("DEM 1m hillshade ...");

            exec(resources(), LOG,
                    GDALDEM,
                    "hillshade",
                    workspace().getFile("DEM_clip1re.tif"),
                    workspace().getFile("GDAL_1mDEM_hillshade.tif"),
                    "-of",
                    "GTiff"
            );
            checkExist(workspace().getDir(), "GDAL_1mDEM_hillshade.tif");

        }

        selected_dem = dem3;
        pixel_width = 3;

        setProgress("Checking 3m DEM ...");

        try {
            exec(resources(), LOG,
                    GDAL,
                    selected_dem,
                    workspace().getFile("DEM_clip3.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"
            );

            if (new File(workspace().getDir(), "DEM_clip3.tif").exists()) {
                exec(resources(), LOG,
                        PYTHON,
                        resources().getFile(GRID_INFO),
                        workspace().getFile("DEM_clip3.tif"),
                        workspace().getFile("pixel_test3.txt")
                );
                if (new File(workspace().getDir(), "pixel_test3.txt").exists()) {
                    try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("pixel_test3.txt")))) {
                        int c = 0;
                        String line;
                        while ((line = br.readLine()) != null) {
                            if (c == 0) {
                                grid_value3 = Double.parseDouble(line);
                                break;
                            }
                        }
                    }
                }
                if (grid_value3 > grid_cover_threshold) {
                    m3_check = true;
                }
            }

        } catch (FileNotFoundException e) {
        } catch (IOException e) {
        } finally {
            setProgress("3m DEM! ...");
        }

        if (m3_check && output) {
            setProgress("still checking 3m 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_clip3.tif"),
                    workspace().getFile("DEM_clip3re.tif")
            );

            checkExist(workspace().getDir(), "DEM_clip3re.tif");

            setProgress("DEM 3m hillshade ...");

            exec(resources(), LOG,
                    GDALDEM,
                    "hillshade",
                    workspace().getFile("DEM_clip3re.tif"),
                    workspace().getFile("GDAL_3mDEM_hillshade.tif"),
                    "-of",
                    "GTiff"
            );
            checkExist(workspace().getDir(), "GDAL_3mDEM_hillshade.tif");
        }

        setProgress("Finished 1m & 3m DEM cover checking.");
    }

    @Override
    protected void postProcess() throws Exception {
        results().put("area_km2", area);
        results().put("1m_DEM", m1_check);
        results().put("3m_DEM", m3_check);

        if (output) {
            ZipFiles.zip(workspace().getFile("output.zip"),
                    (Object) workspace().getFiles("*.tif", "*.csv", "*.kml", "*.kmz", "*.json",
                            "*.geojson", "*.asc", "*.shp", "*.dbf", "*.prj", "*.shx", "*.dat"));
            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

}