V1_0.java [src/java/m/watershed/checkDEM3m_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.checkDEM3m_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/checkDEM3m_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 m3_check = false;
boolean output = false;
double area = 0;
String latlong = "";
Integer resolution = 10;
double grid_cover_threshold = -1.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 selected_dem = dem3;
Integer pixel_width = 3;
output = parameter().getBoolean("out_zip", false);
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("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
}