V1_1.java [src/java/m/rse/wepot] Revision: 2282d80b024781c7d251f5f85e0b5f14ced1dc88 Date: Fri Jan 08 12:54:25 MST 2016
package m.rse.wepot;
import csip.ModelDataService;
import csip.ServiceException;
import csip.annotations.Polling;
import java.io.IOException;
import java.sql.ResultSet;
import java.sql.SQLException;
import java.util.HashMap;
import java.util.ArrayList;
import java.util.logging.Level;
import java.util.logging.Logger;
import javax.ws.rs.Path;
import m.rse.cfactor.utils.Const;
import m.rse.cfactor.utils.PGTools;
import oms3.annotations.Description;
import oms3.annotations.Name;
import org.codehaus.jettison.json.JSONArray;
import org.codehaus.jettison.json.JSONException;
import org.codehaus.jettison.json.JSONObject;
/**
*
* @author Shaun Case
* @author wl
*/
@Name("wepot")
@Description("This service computes water and wind erodibility potentials for an area of analysis (AoA). The service clips SSURGO soil mapunits with AoA geometry, determines the dominant soil component in the AoA, gets parameters from the SSURGO component table, including a climate factor from the C Factor layer, and computes the following equations: Wind Erosion Potential = C*I/T ; Water Erosion Potential = K*(LS)/T")
@Path("m/wepot/1.1")
@Polling(first = 5000, next = 2000)
public class V1_1 extends ModelDataService {
//Do we need these for RSE-02? (taken from RSE-01)
//static final String ALASKA_REGION_WARNING_MSG = "Coordinate is in Alaska region. C-Factor value determination is pending further design specification from USDA-NRCS";
//static final String SERVICE_DESC = "Climatic erosivity reflecting wind speed and surface soil moisture; value 100 at Garden City, Kansas";
//static final double ALASKA_MIN_LATITUDE = 50.00;
//static final double UNKNOWN_CFACTOR = -9999.99;
static final int JSON_LATITUDE = 1;
static final int JSON_LONGITUDE = 0;
static final String RSE_AOAID = "AoAId";
static final String AOA_GEOMETRY = "aoa_geometry";
private String AoAId;
private JSONObject aoaGeometry;
private PGTools.PolygonLatLon polygonData;
private String error_msg;
//Result variables
private double aoa_cfactor;
private String aoa_dom_comp;
private String aoa_dom_compname;
private double aoa_ifactor;
private double aoa_tfactor;
private double aoa_lsfactor;
private double aoa_kfactor;
private double aoa_wind_ep;
private double aoa_water_ep;
double cfactorraster = -1;
@Override
protected void preProcess() {
this.error_msg = "";
this.AoAId = "Error";
this.polygonData = null;
this.aoaGeometry = null;
this.aoa_cfactor = Const.UNKNOWN_CFACTOR;
this.aoa_dom_comp = "";
this.aoa_dom_compname = "";
this.aoa_ifactor = Const.UNKNOWN_CFACTOR;
this.aoa_tfactor = Const.UNKNOWN_CFACTOR;
this.aoa_lsfactor = Const.UNKNOWN_CFACTOR;
this.aoa_kfactor = Const.UNKNOWN_CFACTOR;
this.aoa_wind_ep = Const.UNKNOWN_CFACTOR;
this.aoa_water_ep = Const.UNKNOWN_CFACTOR;
try {
//These functions will throw a usuable user message if they cannot find these values.
// Just catch it and return it to the user.
this.AoAId = getStringParam(V1_1.RSE_AOAID);
this.aoaGeometry = getJSONParam(V1_1.AOA_GEOMETRY);
} catch (ServiceException ex) {
this.error_msg = "Cannot process request JSON: " + ex.getMessage();
}
}
@Override
protected String process() {
try {
if (this.error_msg.isEmpty()) {
JSONArray features = this.aoaGeometry.optJSONArray("features");
if (this.buildPolygon(features)) {
//Got a valid Polygon now.
//HINT: You can get the WKT from the PolygonLatLon object to build an SQL query via the "toWKT" funciton,
// and then encasing that WKT in the typical "geometry::STPolyFromText(<KWT goes here>,4326)", etc.
LOG.info("THE POLYGON AS WE KNOW IT IS=" + this.polygonData.toWKT());
String cfactor_polygon = "geometry::STPolyFromText(" + this.polygonData.toWKT() + ",4326)";
LOG.info("getting centroid for polygon now");
PGTools.Centroid centroid = PGTools.getCentroid(cfactor_polygon);
LOG.info("About to get C Factor Raster value for lat=" + centroid.lat + " long=" + centroid.lon + " ");
this.aoa_cfactor = PGTools.getCFactorRaster(centroid.lat, centroid.lon, getWorkspaceDir()) / 100.0;
//Call SSURGO gemoetry intersect code here ( ssurgo.soilmu_a table )
//Call SSURGO component lookups here ( ssurgo.component table )
if (this.componentLookup(this.ssurgoIntersect(this.polygonData.toWKT()), this.isPalouseArea(this.polygonData.toWKT()))) {
//Calc wind and water ep here...
this.aoa_wind_ep = this.aoa_cfactor * this.aoa_ifactor / this.aoa_tfactor;
this.aoa_water_ep = this.aoa_kfactor * this.aoa_lsfactor / this.aoa_tfactor;
}
}
}
} catch (JSONException ex) {
this.error_msg += "Cannot proceed with processing the request: " + ex.getMessage();
} catch (SQLException se) {
this.error_msg += "SQL exception getting cfactor: " + se.getMessage();
} catch (ServiceException sx) {
this.error_msg += "Service exception getting cfactor: " + sx.getMessage();
}
catch (IOException ioe) {
this.error_msg += "IO Exception getting cFactor value from raster geotiff layer: " + ioe.getMessage();
}
return (this.error_msg.isEmpty() ? EXEC_OK : this.error_msg);
}
@Override
protected void postProcess() throws Exception {
if (this.error_msg.isEmpty()) {
putResult("aoa_id", this.AoAId);
putResult("aoa_dom_comp", this.aoa_dom_comp);
putResult("aoa_dom_compname", this.aoa_dom_compname);
putResult("aoa_cfactor", this.aoa_cfactor);
putResult("aoa_ifactor", this.aoa_ifactor);
putResult("aoa_tfactor", this.aoa_tfactor);
putResult("aoa_lsfactor", this.aoa_lsfactor);
putResult("aoa_kfactor", this.aoa_kfactor);
putResult("aoa_wind_ep", this.aoa_wind_ep);
putResult("aoa_water_ep", this.aoa_water_ep);
}
}
//Private functions
private Boolean buildPolygon(JSONArray features) throws JSONException {
Boolean ret_val = false;
if (null != features) {
if (features.length() > 0) {
JSONObject feature = features.getJSONObject(0).optJSONObject("geometry");
if (this.isGeometryPolygonType(feature)) {
JSONArray polygon = feature.optJSONArray("coordinates");
if (null != polygon) {
if (polygon.length() > 0) {
this.polygonData = readPolygonCoordinates(polygon.getJSONArray(0));
if (null != this.polygonData) {
if (!this.polygonData.isValid()) {
this.error_msg += "Invalid latitude and/or longitude data contained in this polygon. Cannot proceed with processing of this request. ";
} else {
ret_val = true;
}
}//No else needed here, the error message will be built in readPolygonCoordinates()
} else {
this.error_msg = "No coordinates found associated with the polygon specified in feature collection number: 1 . ";
}
} else {
this.error_msg = "No coordinates found associated with the polygon specified in feature collection number: 1 . ";
}
}//No else needed here, the error message will be built in isGeometryPolygonType()
} else {
this.error_msg = "No geometry found associated with this feature. ";
}
} else {
this.error_msg = "Cannot process request JSON, missing features. ";
}
return ret_val;
}
private Boolean isGeometryPolygonType(JSONObject geometry) {
Boolean ret_val = false;
if (null != geometry) {
String geometryType;
geometryType = geometry.optString("type");
if (!geometryType.isEmpty()) {
if (geometryType.toLowerCase().equals("polygon")) {
ret_val = true;
} else {
this.error_msg = "No valid geometry type found in the feature collection number: 1 . Looking for 'Polygon'. ";
}
} else {
this.error_msg = "No geometry type specified in the feature collection number: 1 . ";
}
} else {
this.error_msg = "No geometry found in the feature collection number: 1 . ";
}
return ret_val;
}
private PGTools.PolygonLatLon readPolygonCoordinates(JSONArray shape) throws JSONException {
PGTools.PolygonLatLon tPolygon = new PGTools.PolygonLatLon();
if (null != shape) {
for (int k = 0; k < shape.length(); k++) {
JSONArray jPoint = shape.getJSONArray(k);
tPolygon.add(jPoint.getDouble(V1_1.JSON_LATITUDE), jPoint.getDouble(V1_1.JSON_LONGITUDE));
}
} else {
this.error_msg = "Cannot find the latitude and longitude values for this polygon. ";
tPolygon = null;
}
return tPolygon;
}
// If the WKTPolygon of interest intersects any polygons in the Palouse, MLRA-9, region, this funciton will return true.
private boolean isPalouseArea(String WKTPolygon) {
boolean ret_val = false;
if (null != WKTPolygon) {
String query;
query = "SELECT id"
+ " FROM ssurgo.mlra_9 "
+ " WHERE geometry::STGeomFromText(component_boundary.STAsText(),4326).STIntersects(geometry::STGeomFromText(" + WKTPolygon + ",4326))=1"
+ " AND geometry::STGeomFromText(" + WKTPolygon + ",4326).STIsValid()=1;";
LOG.log(Level.INFO, "ssurgo Palouse query={0}", query);
try (
java.sql.Connection conn = PGTools.getConnection("cfactor", LOG);
java.sql.Statement statement = conn.createStatement();) {
ResultSet results = statement.executeQuery(query);
LOG.info(" past the query execute");
if (null != results) {
if (results.next()) {
ret_val = true;
}
}
} catch (SQLException | ServiceException ex) {
this.error_msg += "Cannot continue processing this request in the find Palouse region procedure: " + ex.getMessage();
StackTraceElement ste[] = ex.getStackTrace();
LOG.severe(ex.toString());
for (StackTraceElement e : ste) {
LOG.severe(e.toString());
}
}
}
return ret_val;
}
// Returns the intersected mukey with the largest area.
private String ssurgoIntersect(String WKTPolygon) {
HashMap<String, V1_1.ssurgoAttributes> mukeyAttributes;
ArrayList<V1_1.ssurgoAttributes> attributeList;
//String polygonText = " ST_PolygonFromText(" + WKTPolygon + ") ";
String ret_val = null;
if (null != WKTPolygon) {
String query;
mukeyAttributes = new HashMap<>();
attributeList = new ArrayList<>();
// postgis sql
// query = "SELECT m.areasymbol, m.musym, m.mukey,"
// + "st_area(st_transform(st_intersection(ST_PolygonFromText(" + WKTPolygon + ", 4326), m.the_geom::geography::geometry ),3541))/43560 as sizeIntersectionAcres "
// + " FROM ssurgo.soilmu_a as m "
// + " WHERE ST_Intersects(" + polygonText + ", m.the_geom ) "
// + " AND st_isvalid( m.the_geom) AND st_isvalid(" + polygonText + ");";
// sql server spatial sql
query
= /* "SELECT areasymbol, musym, mukey, geography::STPolyFromText(" + WKTPolygon + ",4326).MakeValid().ReorientObject().STIntersection(soilpoly).STArea() / 4046.86 as sizeIntersectionAcres "
+ "FROM "
+ "(SELECT m.areasymbol, m.musym, m.mukey, geography::STPolyFromText(m.the_geom.STAsText(),4326).MakeValid().ReorientObject() as soilpoly "
+ "FROM ssurgo.ssurgo.soilmu_a as m "
+ "WITH (index(geom_sidx)) "
+ "WHERE m.the_geom.STIntersects(geometry::STPolyFromText(" + WKTPolygon + ",0)) = 1 "
+ "AND m.the_geom.STIsValid()=1 "
+ "AND (geometry::STPolyFromText(" + WKTPolygon + ",0).STIsValid())=1) "
+ "as a;";
*/ "SELECT areasymbol, musym, mukey, geography::STGeomFromText(intersectPoly.STAsText(), 4326).MakeValid().STArea() / 4046.86 as sizeIntersectionAcres"
+ " FROM"
+ " (SELECT m.areasymbol, m.musym, m.mukey, "
+ " m.the_geom.STIntersection(geometry:: STGeomFromText (" + WKTPolygon + ",0)).MakeValid() as intersectPoly"
+ " FROM ssurgo.soilmu_a as m "
+ " WITH (index(geom_sidx))"
+ " WHERE m.the_geom.STIntersects(geometry:: STGeomFromText (" + WKTPolygon + ",0)) = 1"
+ " AND m.the_geom.STIsValid()=1"
+ " AND (geometry:: STGeomFromText (" + WKTPolygon + ",0).STIsValid())=1)"
+ " as a;";
LOG.info("ssurgo query=" + query);
try (
java.sql.Connection conn = PGTools.getConnection("cfactor", LOG);
java.sql.Statement statement = conn.createStatement();) {
ResultSet results = statement.executeQuery(query);
LOG.info(" past the query execute");
if (null != results) {
while (results.next()) {
V1_1.ssurgoAttributes tAttributes;
String tKey = results.getString("areasymbol") + ":" + results.getString("mukey");
tAttributes = mukeyAttributes.get(tKey);
//Merge, i.e. sum, areas of duplicate AoA:mukey combinations
if (null != tAttributes) {
tAttributes.gid_area += results.getDouble("sizeIntersectionAcres");
} else {
tAttributes = new V1_1.ssurgoAttributes();
tAttributes.gid = results.getString("musym");
tAttributes.AoA_Id = results.getString("areasymbol");
tAttributes.mukey = results.getString("mukey");
tAttributes.gid_area = results.getDouble("sizeIntersectionAcres");
mukeyAttributes.put(tKey, tAttributes);
attributeList.add(tAttributes);
}
}
// Have a unique set of mukeys now, let's find the largest
double largestArea = 0.0;
for (V1_1.ssurgoAttributes tAttribute : attributeList) {
if (largestArea < tAttribute.gid_area) {
largestArea = tAttribute.gid_area;
ret_val = tAttribute.mukey;
}
}
} else {
this.error_msg += "No results from the intersect query for this geometry. ";
}
} catch (SQLException | csip.ServiceException ex) {
this.error_msg += "Cannot continue processing this request in the intersect procedure: " + ex.getMessage();
StackTraceElement ste[] = ex.getStackTrace();
LOG.severe(ex.toString());
for (StackTraceElement e : ste) {
LOG.severe(e.toString());
}
ret_val = null;
}
}
return ret_val;
}
private Boolean componentLookup(String mukey, boolean isPalouseRegion) {
String query;
ResultSet results;
Boolean ret_val = false;
double slope_r = 0.0;
double theta = 0.0;
double lambda = 0.0;
double length_m = 0.0;
double[] slopeDistance = {200.0, 300.0, 200.0, 180.0, 160.0, 150.0, 140.0, 130.0, 125.0, 120.0, 110.0, 100.0, 90.0, 80.0, 70.0};
double[] palouseSlopeDistance = {350.0, 275.0, 225.0, 175.0, 150.0, 125.0, 100.0};
boolean palouseRegion = isPalouseRegion;
if ((null != mukey) && !mukey.isEmpty()) {
try (
java.sql.Connection conn = PGTools.getConnection("cfactor", LOG);
java.sql.Statement statement = conn.createStatement();) {
query = "SELECT component.cokey, component.compname, component.wei as ifact, component.slope_r, component.slopelenusle_r as lsfact, component.tfact, component.comppct_r "
+ " FROM ssurgo.component "
+ " WHERE component.mukey='" + mukey + "' "
+ " AND component.slope_r is not null"
+ " ORDER BY component.comppct_r DESC;";
LOG.info(" ssurgo query=" + query);
results = statement.executeQuery(query);
LOG.info(" past the query execute");
if (null != results) {
results.next(); //Only need the first result row, since it will contain the larges comppct_r value, due to the order by clause.
this.aoa_dom_comp = results.getString("cokey");
this.aoa_dom_compname = results.getString("compname");
this.aoa_ifactor = results.getDouble("ifact");
slope_r = results.getDouble("slope_r");
if (!palouseRegion) {
if (slope_r >= 18) {
lambda = 50;
} else if (slope_r >= 16) {
lambda = 60;
} else if (slope_r >= 0 && slope_r < 1.0) { //maybe should just be <1.0
lambda = 100;
} else {
lambda = slopeDistance[((int) slope_r) - 1]; // Maybe should be ((int) round(slope_r)) - 1
}
} else if (slope_r < 6) {
lambda = palouseSlopeDistance[0];
} else if (slope_r < 11) {
lambda = palouseSlopeDistance[1];
} else if (slope_r < 16) {
lambda = palouseSlopeDistance[2];
} else if (slope_r < 21) {
lambda = palouseSlopeDistance[3];
} else if (slope_r < 26) {
lambda = palouseSlopeDistance[4];
} else if (slope_r < 36) {
lambda = palouseSlopeDistance[5];
} else if (slope_r < 46) {
lambda = palouseSlopeDistance[6];
}
length_m = ((slope_r >= 5.0) ? 0.5 : (slope_r >= 3.5) ? 0.4 : (slope_r >= 1.0) ? 0.3 : 0.2);
theta = Math.atan(slope_r / 100.0);
this.aoa_lsfactor = Math.pow((lambda / 72.6), length_m) * (65.41 * Math.pow(Math.sin(theta), 2.0) + 4.56 * Math.sin(theta) + 0.065); //results.getDouble("lsfact");
this.aoa_tfactor = results.getDouble("tfact");
results.close();
query = "SELECT chorizon.cokey, chorizon.kffact, chorizon.kwfact "
+ "FROM ssurgo.chorizon "
+ "WHERE ssurgo.chorizon.cokey='" + this.aoa_dom_comp + "'";
LOG.info(" ssurgo query=" + query);
results = statement.executeQuery(query);
LOG.info(" past the query execute");
if (null != results) {
results.next();
this.aoa_kfactor = results.getDouble("kffact");
if (results.wasNull()) {
this.aoa_kfactor = results.getDouble("kwfact");
if (results.wasNull()) {
this.aoa_kfactor = 0.02;
}
}
ret_val = true;
} else {
this.error_msg += "Could not find a matching cokey for this intersected mapunit in the cHorizon table.";
}
} else {
this.error_msg += "Could not find a matching cokey for this intersected mapunit.";
}
} catch (SQLException | csip.ServiceException ex) {
this.error_msg += "Could not continue processing the component lookups for this mapunit intersect: " + ex.getMessage();
}
}
return ret_val;
}
//Inner Classes/Structures go here...
/**
*
*/
public static class ssurgoAttributes {
public String gid;
public String AoA_Id;
public String mukey;
public double gid_area;
}
}