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
}