V1_0.java [src/java/m/watershed/raster_stats] 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.raster_stats;
import csip.ModelDataService;
import csip.api.server.PayloadParameter;
import csip.api.server.ServiceException;
import csip.annotations.Description;
import csip.annotations.Name;
import csip.annotations.Resource;
import javax.ws.rs.Path;
import csip.utils.ZipFiles;
import java.io.*;
import java.nio.charset.Charset;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import m.watershed.Resources;
import static m.watershed.Resources.GDALCALC;
import static m.watershed.Resources.GDALTRANSLATE;
import static m.watershed.Resources.OGR;
import static m.watershed.Resources.PYTHON;
import static m.watershed.Resources.ZONALSTATS_MAMC;
import static m.watershed.Utils.checkExist;
import static m.watershed.Utils.exec;
import org.apache.commons.io.FileUtils;
import org.codehaus.jettison.json.JSONObject;
/**
* raster_stats
*
* @author HK (using OD's template)
*/
@Name("Raster Stats by Vector data")
@Description("Example ... input data")
@Path("m/raster_stats/1.0")
@Resource(from = Resources.class)
public class V1_0 extends ModelDataService {
static final String BOUNDARY = "boundary";
static final String RESULTID = "result_ID";
static final String STATS = "stats";
static final String RASTER = "raster";
static final String UTMFILE = "utm_file";
boolean shp_file = false;
String shp_name = "";
String used_stats = "";
String used_raster = "";
String utm_string = "";
Map<String, String> Results = new HashMap<>();
JSONObject boundary = null;
double Result = Double.NaN;
@Override
protected void doProcess() throws Exception {
PayloadParameter inputPayload = parameter();
if (inputPayload.has(BOUNDARY) && inputPayload.has(RESULTID) && inputPayload.has(UTMFILE)) {
String b = parameter().getString(BOUNDARY);
String ID = parameter().getString(RESULTID);
used_stats = parameter().getString(STATS, "ALL");
used_raster = parameter().getString(RASTER, "help.asc");
if (b.endsWith(".shp")) {
shp_file = true;
File input_file = attachments().getFile(b);
checkExist(workspace().getDir(), shp_name);
shp_name = input_file.getName();
} else {
boundary = parameter().getJSONArray(BOUNDARY).getJSONObject(0);
FileUtils.writeStringToFile(workspace().getFile("boundary.geojson"), boundary.toString(), "UTF-8");
checkExist(workspace().getDir(), "boundary.geojson");
shp_name = "field_extent_utm.shp";
}
run(ID);
} else {
throw new ServiceException("No result-id or No UTM-file or No Extent or no GeoJSON Geometry input listed!");
}
}
private void run(String ID) throws Exception {
int threads = Resources.mpiThreads(parameter());
setProgress("Started...");
String condition_string = "--calc=A*(A>=-900)";
if (used_stats.equals("ALL")) {
condition_string = "--calc=A*(A>=-900)";
}
if (used_stats.equals("NEGATIVE")) {
condition_string = "--calc=A*(A<0)";
}
if (used_stats.equals("POSITIVE")) {
condition_string = "--calc=A*(A>0)";
}
if (used_stats.equals("ZERO")) {
condition_string = "--calc=A*(A==0.0)";
}
checkExist(workspace().getDir(), parameter().getString(UTMFILE));
BufferedReader br = null;
String utmline = "";
br = new BufferedReader(new FileReader(workspace().getFile(parameter().getString(UTMFILE))));
utmline = br.readLine();
if (!shp_file) {
setProgress("No shp file, generate one...");
workspace().writeString("boundary.geojson", boundary.toString());
exec(resources(), LOG,
OGR,
"-t_srs", "EPSG:4326",
"-f", "ESRI Shapefile",
workspace().getFile("field_extent0.shp"),
workspace().getFile("boundary.geojson")
);
checkExist(workspace().getDir(), "field_extent0.shp");
exec(resources(), LOG,
OGR,
"-explodecollections",
"-skipfailures",
"-f", "ESRI Shapefile",
"-dialect", "sqlite",
"-sql", "select ST_buffer(geometry,0), * from field_extent0",
workspace().getFile("field_extent.shp"),
workspace().getFile("field_extent0.shp")
);
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");
setProgress("Reproject ...");
exec(resources(), LOG,
OGR,
"-t_srs", "" + utmline,
"-f", "ESRI Shapefile",
workspace().getFile("field_extent_utm.shp"),
workspace().getFile("field_extent.shp")
);
checkExist(workspace().getDir(), "field_extent_utm.shp");
} else {
setProgress("Use provided shp file...");
}
setProgress("ASCII to UTM GeoTIFF...");
setProgress("Translate FlowPathLoss ASCII to GeoTIFF...");
exec(resources(), LOG,
GDALTRANSLATE,
"-of",
"GTiff",
"-ot", "Float32",
"-a_srs", "" + utmline,
workspace().getFile(used_raster),
workspace().getFile("output0.tif"),
"-a_nodata", "-999"
);
checkExist(workspace().getDir(), "output0.tif");
setProgress("Filter flowpath raster...");
exec(resources(), LOG,
GDALCALC,
"-A", workspace().getFile("output0.tif"),
"--overwrite",
"--outfile=" + workspace().getFile("output.tif"),
"--NoDataValue=0",
condition_string
);
checkExist(workspace().getDir(), "output.tif");
String resultfile = "";
setProgress("Stats...");
exec(resources(), LOG,
PYTHON,
resources().getFile(ZONALSTATS_MAMC),
workspace().getFile(shp_name),
workspace().getFile("output.tif"),
ID,
"stats_",
"False",
"-999"
);
resultfile = "stats_results.csv";
checkExist(workspace().getDir(), "stats_results.csv");
br = new BufferedReader(new FileReader(workspace().getFile(resultfile)));
Integer c = 0;
String line;
while ((line = br.readLine()) != null) {
if (c > 0) {
String[] infos = line.split(",");
Results.put(infos[0], infos[1]);
}
c++;
}
// pullResult();
setProgress("Finished.");
}
// 3) provide the temperature as a result.
@Override
protected void postProcess() throws Exception {
// results().put("area_km2", area);
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"));
}
}
protected void pullResult(String prev) throws IOException, ServiceException {
File resultFile = workspace().getFile(prev + "_results.csv");
List<String> results = FileUtils.readLines(resultFile, Charset.defaultCharset());
if (results.size() >= 2) {
String average = results.get(1);
String[] values = average.split(",");
Result = Double.parseDouble(values[1]);
} else {
throw new ServiceException("No values are found in the resulting csv file.");
}
}
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
}