V1_0.java [src/java/m/watershed/average_aspect] 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.average_aspect;
import csip.Check;
import csip.ModelDataService;
import csip.api.server.ServiceException;
import csip.annotations.Description;
import csip.annotations.Name;
import csip.annotations.Resource;
import static csip.annotations.ResourceType.OUTPUT;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.util.HashMap;
import java.util.Map;
import javax.ws.rs.Path;
import m.watershed.Resources;
import static m.watershed.Resources.GDAL;
import static m.watershed.Resources.GDALDEM;
import static m.watershed.Resources.LOCATION_MNT_DATA1;
import static m.watershed.Resources.LOCATION_MNT_CSIP_DEM;
import static m.watershed.Resources.MPI;
import static m.watershed.Resources.OGR;
import static m.watershed.Resources.PITREMOVE;
import static m.watershed.Resources.PYTHON;
import static m.watershed.Resources.ZONALSTATS_MEAN;
import static m.watershed.Resources.ZONALSTATS_MAX;
import static m.watershed.Utils.exec;
import static m.watershed.Utils.checkExist;
import org.apache.commons.io.FileUtils;
import org.codehaus.jettison.json.JSONObject;
/**
* average_aspect
*
* @author @author HK (using JK's template (using OD's template))
*/
@Name("average_aspect")
@Description("Example of average_aspect")
@Path("m/average_aspect/1.0")
@Resource(from = Resources.class)
//@Resource(file = "*.tif *.csv *.kml *.kmz *.geojson *.asc *.shp *.txt *.dbf *.prj *.shx", type = OUTPUT)
@Resource(file = "*.csv", type = OUTPUT)
public class V1_0 extends ModelDataService {
static final String BOUNDARY = "boundary";
static final String RESULTID = "result_ID";
static final String DEM_RES = "DEM_RES";
static final String FILL = "FILL";
static final String SERVICETHREADS = "Threads";
boolean shp_file = false;
String shp_name = "";
Integer resolution = 10;
Map<String, String> GDALslope = new HashMap<>();
Map<String, String> GDALaspect = new HashMap<>();
Map<String, String> DEMelevation = new HashMap<>();
Map<String, String> maxDEMelevation = new HashMap<>();
@Override
protected void doProcess() throws Exception {
resolution = parameter().getInt(DEM_RES, 10);
new Check().positive().forNumber(resolution);
if (parameter().has(BOUNDARY) && parameter().has(RESULTID)) {
String b = parameter().getString(BOUNDARY);
if (b.endsWith("shp")) {
shp_file = true;
File input_file = attachments().getFile(b);
shp_name = input_file.getName();
} else {
JSONObject boundary = parameter().getJSONArray(BOUNDARY).getJSONObject(0);
FileUtils.writeStringToFile(workspace().getFile("boundary.geojson"), boundary.toString(), "UTF-8");
}
String ID = parameter().getString(RESULTID);
// optional parameter
boolean filling = parameter().getBoolean(FILL, false);
int threads = parameter().getInt(SERVICETHREADS, 2);
run(ID, filling, resolution, threads);
} else {
throw new ServiceException("Geometry Input File not found and/or Result-ID not found !! ");
}
}
/**
*
* @param boundary
* @param ID
* @param resolution
* @throws Exception
*
*/
private void run(String ID, boolean filling, int resolution, int threads) throws Exception {
String dem30 = LOCATION_MNT_DATA1 + "/DEM.vrt";
String dem10 = LOCATION_MNT_CSIP_DEM + "/DEM.vrt";
String dem3 = LOCATION_MNT_CSIP_DEM + "/DEM_1_9_arc/DEM_1_9.vrt";
String dem1 = LOCATION_MNT_CSIP_DEM + "/DEM_1m/DEM_1m.vrt";
// String selected_dem = dem10;
// Integer pixel_width = 10;
// String buffer = "2100";
//
// if (resolution < 10) {
// selected_dem = dem3;
// pixel_width = 3;
// }
// if (resolution == 1) {
// selected_dem = dem1;
// pixel_width = 1;
// }
// if (resolution > 10) {
// selected_dem = dem30;
// pixel_width = 30;
// }
//
// setProgress("Clip DEM...");
// exec(resources(), LOG,
// GDAL,
// selected_dem,
// workspace().getFile("DEM_clip0.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"
// );
// checkExist(workspace().getDir(), "DEM_clip0.tif");
//
// setProgress("Reproject 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_clip0.tif"),
// workspace().getFile("DEM_clip.tif")
// );
//
// checkExist(workspace().getDir(), "DEM_clip.tif");
if (!shp_file) {
checkExist(workspace().getDir(), "boundary.geojson");
exec(resources(), LOG,
OGR,
"-t_srs", "EPSG:5070",
"-f", "ESRI Shapefile",
workspace().getFile("boundary_pre.shp"),
workspace().getFile("boundary.geojson")
);
} else {
checkExist(workspace().getDir(), shp_name);
exec(resources(), LOG,
OGR,
"-t_srs", "EPSG:5070",
"-f", "ESRI Shapefile",
workspace().getFile("boundary_pre.shp"),
workspace().getFile(shp_name)
);
}
checkExist(workspace().getDir(), "boundary_pre.shp");
if (resolution > 10) {
exec(resources(), LOG,
OGR,
"-update",
"-append",
"-dialect", "sqlite",
"-sql", "select ST_buffer(extent(geometry), " + (resolution + 5) + ")," + ID + " from boundary_pre",
workspace().getFile("boundary.shp"),
workspace().getFile("boundary_pre.shp")
);
checkExist(workspace().getDir(), "boundary.shp");
exec(resources(), LOG,
GDAL,
LOCATION_MNT_DATA1 + "/DEM.vrt",
workspace().getFile("DEM_clip0.tif"),
"-crop_to_cutline",
"-cutline",
workspace().getFile("boundary.shp"),
"-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", "1000",
"-wm", "1000"
);
checkExist(workspace().getDir(), "DEM_clip0.tif");
exec(resources(), LOG,
GDAL,
"-t_srs", "EPSG:5070",
"-tr", "30", "30",
"-r", "cubicspline",
workspace().getFile("DEM_clip0.tif"),
workspace().getFile("DEM_clip.tif")
);
checkExist(workspace().getDir(), "DEM_clip.tif");
} else {
exec(resources(), LOG,
OGR,
"-update",
"-append",
"-dialect", "sqlite",
"-sql", "select ST_buffer(extent(geometry), 1)," + ID + " from boundary_pre",
workspace().getFile("boundary.shp"),
workspace().getFile("boundary_pre.shp")
);
checkExist(workspace().getDir(), "boundary.shp");
exec(resources(), LOG,
GDAL,
LOCATION_MNT_CSIP_DEM + "/DEM.vrt",
workspace().getFile("DEM_clip0.tif"),
"-crop_to_cutline",
"-cutline",
workspace().getFile("boundary.shp"),
"-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", "2000",
"-wm", "2000"
);
checkExist(workspace().getDir(), "DEM_clip0.tif");
exec(resources(), LOG,
GDAL,
"-t_srs", "EPSG:5070",
"-tr", "10", "10",
"-r", "cubicspline",
workspace().getFile("DEM_clip0.tif"),
workspace().getFile("DEM_clip.tif"),
"-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", "2000",
"-wm", "2000"
);
checkExist(workspace().getDir(), "DEM_clip.tif");
}
String dem_usage = "DEM_clip.tif";
if (filling) {
exec(resources(), LOG,
MPI,
"--allow-run-as-root",
"--path", resources().getFile(PITREMOVE).getParent(),
"-wdir", workspace().getDir(),
"--oversubscribe",
"-np", threads,
resources().getFile(PITREMOVE).getName(),
"-z", workspace().getFile("DEM_clip.tif"),
"-fel", workspace().getFile("dem_filled.tif")
);
checkExist(workspace().getDir(), "dem_filled.tif");
dem_usage = "dem_filled.tif";
}
exec(resources(), LOG,
GDALDEM,
"slope",
workspace().getFile(dem_usage),
workspace().getFile("GDALslope_deg.tif"),
//"-p",
"-of",
"GTiff"//,
//"-compute_edges"
);
exec(resources(), LOG,
PYTHON,
resources().getFile(ZONALSTATS_MEAN),
workspace().getFile("boundary_pre.shp"),
workspace().getFile("GDALslope_deg.tif"),
ID,
"avg_GDAL_slope_deg_",
"False"
);
exec(resources(), LOG,
GDALDEM,
"aspect",
workspace().getFile(dem_usage),
workspace().getFile("GDALaspect_deg.tif"),
//"-alg", "ZevenbergenThorne",
"-zero_for_flat",
"-of",
"GTiff"//,
//"-compute_edges"
);
exec(resources(), LOG,
PYTHON,
resources().getFile(ZONALSTATS_MEAN),
workspace().getFile("boundary_pre.shp"),
workspace().getFile("GDALaspect_deg.tif"),
ID,
"avg_GDAL_aspect_deg_",
"False"
);
exec(resources(), LOG,
PYTHON,
resources().getFile(ZONALSTATS_MEAN),
workspace().getFile("boundary_pre.shp"),
workspace().getFile(dem_usage),
ID,
"avg_DEM_elevation_m_",
"False"
);
exec(resources(), LOG,
PYTHON,
resources().getFile(ZONALSTATS_MAX),
workspace().getFile("boundary_pre.shp"),
workspace().getFile(dem_usage),
ID,
"max_DEM_elevation_m_",
"False"
);
checkExist(workspace().getDir(),
"avg_GDAL_aspect_deg_results.csv",
"avg_DEM_elevation_m_results.csv",
"max_DEM_elevation_m_results.csv",
"avg_GDAL_slope_deg_results.csv");
try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_GDAL_aspect_deg_results.csv")))) {
int c = 0;
String line;
while ((line = br.readLine()) != null) {
if (c > 0) {
String[] infos = line.split(",");
GDALaspect.put(infos[0], infos[1]);
}
c++;
}
}
try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_GDAL_slope_deg_results.csv")))) {
int c = 0;
String line;
while ((line = br.readLine()) != null) {
if (c > 0) {
String[] infos = line.split(",");
GDALslope.put(infos[0], infos[1]);
}
c++;
}
}
try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("avg_DEM_elevation_m_results.csv")))) {
int c = 0;
String line;
while ((line = br.readLine()) != null) {
if (c > 0) {
String[] infos = line.split(",");
DEMelevation.put(infos[0], infos[1]);
}
c++;
}
}
try (BufferedReader br = new BufferedReader(new FileReader(workspace().getFile("max_DEM_elevation_m_results.csv")))) {
int c = 0;
String line;
while ((line = br.readLine()) != null) {
if (c > 0) {
String[] infos = line.split(",");
maxDEMelevation.put(infos[0], infos[1]);
}
c++;
}
}
}
@Override
protected void postProcess() throws Exception {
//results().put(workspace().getFile("boundary.geojson"), "boundary");
results().put(workspace().getFile("avg_GDAL_aspect_deg_results.csv"), "results");
results().put(workspace().getFile("avg_GDAL_slope_deg_results.csv"), "results");
for (Map.Entry<String, String> entry : GDALslope.entrySet()) {
String key = (String) entry.getKey();
results().put("#ID" + key, key);
results().put("avg_GDAL_slope_deg_" + key, entry.getValue());
results().put("avg_GDAL_aspect_deg_" + key, GDALaspect.get(key));
results().put("avg_DEM_elevation_m_" + key, DEMelevation.get(key));
results().put("max_DEM_elevation_m_" + key, maxDEMelevation.get(key));
}
}
}