JsonRegion.java [src/java/m/weps] Revision: Date:
/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package m.weps;
import java.io.File;
import java.io.IOException;
import java.util.logging.Logger;
import javax.ws.rs.Path;
import org.apache.commons.io.FileUtils;
import org.json.JSONArray;
import org.json.JSONException;
import org.json.JSONObject;
import org.json.JSONTokener;
import static util.ErosionConst.*;
/**
*
* @author jplyon
*/
public class JsonRegion {
protected final Logger LOG = Logger.getLogger(JsonRegion.class.getName());
/**
* Center (centroid) coordinates.
*/
public double center_x = 0.0;
public double center_y = 0.0;
// @todo add major and minor axes
// @todo add moment of inertia
// Offset angle [degrees].
public double region_angle = 0.0;
// The computed area, in degrees squared.
// The accessor function getArea() will return square meters.
private double region_area = 0.0;
/**
* Representative height and width. These multiply together to produce the
* area region_area.
*
* @see acccessors getHeight(), getWidth() which return these lengths in
* meters, corrected for geographic location.
*/
private double region_height = 0.0, region_width = 0.0;
/**
* The type of feature this GeoJSON represents. Currently handled are
* "Polygon", "Line".
*/
private String featureType = null;
{
}
/**
* Get the area of the region [meters^2]. This has to convert from the
* internal units of square degrees.
*
* @return The area of the region. This is the same as the area of the
* representative rectangle.
*/
public double getArea() {
// Get the length of one degree of arc for lat and long in meters.
// Note that these depend on the location on the earth.
double lat_meters = getLatMeters(center_x, center_y);
double long_meters = getLongMeters(center_x, center_y);
double mult = lat_meters * long_meters;
double area = region_area * mult;
return Math.abs(area);
}
/**
* Get the height of the representative rectangle for the region [meters].
* This has to convert from the internal units of degrees.
*
* @return The height of the representative rectangle [meters].
*/
public double getHeight() {
// Get the length of one degree of arc for lat and long in meters.
// Note that these depend on the location on the earth.
double lat_meters = getLatMeters(center_x, center_y);
LOG.info("IN GET HEIGHT");
LOG.info("get lat meters=" + lat_meters);
double long_meters = getLongMeters(center_x, center_y);
LOG.info("get long meters=" + long_meters);
// Now rotate these lengths to get the length in the principal direction
// of the representative rectangle (height).
double height_mult = lat_meters * Math.cos(region_angle * (Math.PI / 180)) + long_meters * Math.sin(region_angle * (Math.PI / 180));
LOG.info("height_mult=" + height_mult);
LOG.info("region_angle=" + region_angle);
double height_m = region_height * height_mult;
LOG.info("height_m=" + height_m);
return Math.abs(height_m);
}
/**
* Get the width of the representative rectangle for the region [meters].
* This has to convert from the internal units of degrees.
*
* @return The width of the representative rectangle [meters].
*/
public double getWidth() {
// Get the length of one degree of arc for lat and long in meters.
// Note that these depend on the location on the earth.
double lat_meters = getLatMeters(center_x, center_y);
double long_meters = getLongMeters(center_x, center_y);
// Now rotate these lengths to get the length in the principal direction
// of the representative rectangle (width).
double width_mult = long_meters * Math.cos(region_angle * (Math.PI / 180)) + lat_meters * Math.sin(region_angle * (Math.PI / 180));
double width_m = region_width * width_mult;
return Math.abs(width_m);
}
/**
* Get the length of one degree of arc of latitude, at a given location.
*
* @param center_x The point at which to find the length (longitude).
* @param center_y The point at which to find the length (latitude).
* @return the length in meters.
*/
public double getLatMeters(double center_x, double center_y) {
// Get the length of one degree of arc in meters.
// Note that these depend on the location on the earth.
// http://mathforum.org/library/drmath/view/51879.html
// The Geographic Information Systems FAQ (See Q5.1)
// http://www.faqs.org/faqs/geography/infosystems-faq/
// @note Angles must be expressed in radians.
// We find the distance by taking two points one-half of a degree
// in either direction from the given point.
double R = 636700; // radius of earth [meters]
double delta_lon = 0.0;
double lat_1 = (center_y - 0.5) * (Math.PI / 180.0);
double lat_2 = (center_y + 0.5) * (Math.PI / 180.0);
double delta_lat = lat_2 - lat_1;
double a = Math.sqrt(Math.sin(delta_lat / 2)) + Math.cos(lat_1) * Math.cos(lat_2) * Math.sqrt(Math.sin(delta_lon / 2));
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
double d = R * c;
assert (d > 0);
return d;
}
/**
* Get the length of one degree of arc of longitude, at a given location.
*
* @param center_x The point at which to find the length (longitude).
* @param center_y The point at which to find the length (latitude).
* @return the length in meters.
*/
public double getLongMeters(double center_x, double center_y) {
// Get the length of one degree of arc in meters.
// Note that these depend on the location on the earth.
// http://mathforum.org/library/drmath/view/51879.html
// The Geographic Information Systems FAQ (See Q5.1)
// http://www.faqs.org/faqs/geography/infosystems-faq/
// @note Angles must be expressed in radians.
// We find the distance by taking two points one-half of a degree
// in either direction from the given point.
double R = 636700; // radius of earth [meters]
double delta_lat = 0.0;
double lon_1 = (center_x - 0.5) * (Math.PI / 180.0);
double lon_2 = (center_x + 0.5) * (Math.PI / 180.0);
double delta_lon = lon_2 - lon_1;
double a = Math.sqrt(Math.sin(delta_lat / 2)) + Math.cos(lon_1) * Math.cos(lon_2) * Math.sqrt(Math.sin(delta_lon / 2));
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
double d = R * c;
assert (d > 0);
return d;
}
{
}
/**
* Constructor which takes a GeoJSON object.
*
* @param json GeoJSON object to build a representation for.
* @throws JSONException
*/
public JsonRegion(JSONObject json) throws JSONException {
// @note this function requires a non-null argument
if (json == null) {
throw new IllegalArgumentException("JsonRegion() requires a non-null GeoJSON argument.");
}
JSONArray centroid = null;
// Parse the JSON
JSONArray pointArray = null;
JSONArray features = json.getJSONArray("features");
JSONObject feature = features.getJSONObject(0);
JSONObject geometry = feature.getJSONObject("geometry");
featureType = geometry.getString("type");
if (featureType.equals("Polygon")) {
JSONArray coordinates = geometry.getJSONArray("coordinates");
pointArray = coordinates.getJSONArray(0);
centroid = getPolygonCentroid(pointArray);
region_area = getPolygonArea(pointArray);
} else if (featureType.equals("Multipolygon")) {
// @todo
} else if (featureType.equals("Line")) {
JSONArray coordinates = geometry.getJSONArray("coordinates");
pointArray = coordinates.getJSONArray(0);
centroid = getLineCentroid(pointArray);
region_area = 0.0;
}
// Set the centroid
center_x = centroid.getDouble(0);
center_y = centroid.getDouble(1);
// Set the properties for the representative rectangle.
setRepresentativeRectangle(pointArray);
}
/**
* Calculate the centroid of the JSON array representation of a polygon.
* This will be an array of points, each of which is an array. The first
* point will be repeated at the end of the array.
*
* @param json The JSON representation of the polygon
* @return A JSON array for the centroid point. Returned as [x-coordinate,
* y-coordinate].
*
* @see http://en.wikipedia.org/wiki/Centroid - "Centroid of polygon"
* @see http://en.wikipedia.org/wiki/Polygon_area#Properties
*/
private static JSONArray getPolygonCentroid(JSONArray json) throws JSONException {
/// @note this function requires a non-null argument
if (json == null) {
throw new IllegalArgumentException("JsonRegion() requires a non-null GeoJSON argument.");
}
double centroid_x = 0.0, centroid_y = 0.0;
double polygon_area = 0.0;
// Calculate the centroid.
// @note we skip the last point, which should be a repeat.
// @todo verify that the last point is repeated and skip it conditionally.
int numPoints = json.length() - 1;
for (int i = 0; i < numPoints; i++) {
// first point
double x0 = json.getJSONArray(i).getDouble(0);
double y0 = json.getJSONArray(i).getDouble(1);
double x1 = json.getJSONArray(i + 1).getDouble(0);
double y1 = json.getJSONArray(i + 1).getDouble(1);
// Calculate centroid of the triangle formed by points [[x0,y0], [x1,y1], [0,0]].
centroid_x += (x0 + x1) * (x0 * y1 - x1 * y0);
centroid_y += (y0 + y1) * (x0 * y1 - x1 * y0);
polygon_area += (x0 * y1 - x1 * y0) / 2;
}
// Now divide by the polygon's area to get the true centroid.
centroid_x /= (6 * polygon_area);
centroid_y /= (6 * polygon_area);
JSONArray centroid = new JSONArray();
centroid.put(0, centroid_x);
centroid.put(1, centroid_y);
return centroid;
}
/**
* Calculate the centroid of the JSON array representation a line. This will
* be an array of points, each of which is an array of 2 coordinates.
*
* @param json The JSON representation of the line
* @return A JSON array for the centroid point. Returned as [x-coordinate,
* y-coordinate].
*
*/
private static JSONArray getLineCentroid(JSONArray json) throws JSONException {
/// @note this function requires a non-null argument
if (json == null) {
throw new IllegalArgumentException("getLineCentroid() requires a non-null JSON array.");
}
double centroid_x = 0.0, centroid_y = 0.0;
double total_length = 0.0;
// @note We are cheating, in that we are not using earth coordinates.
// Instead we use the raw latitude and longitude.
// This is sufficiently accurate for regions that are small in extent,
// and far from the poles.
// Calculate the centroid.
// @note we don't skip the last point, which could be a repeat.
int numPoints = json.length() - 1;
for (int i = 0; i < numPoints; i++) {
// first point
double x0 = json.getJSONArray(i).getDouble(0);
double y0 = json.getJSONArray(i).getDouble(1);
double x1 = json.getJSONArray(i + 1).getDouble(0);
double y1 = json.getJSONArray(i + 1).getDouble(1);
double length = Math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0));
// Calculate center of mass of the segment formed by points [[x0,y0], [x1,y1], [0,0]].
centroid_x += length * (x0 + x1) / 2.0;
centroid_y += length * (y0 + y1) / 2.0;
total_length += length;
}
// Now divide by the polygon's area to get the true centroid.
centroid_x /= total_length;
centroid_y /= total_length;
JSONArray centroid = new JSONArray();
centroid.put(0, centroid_x);
centroid.put(1, centroid_y);
return centroid;
}
/**
* Calculate the area of the JSON array representation of a polygon. This
* will be an array of points, each of which is an array. The first point
* will be repeated at the end of the array.
*
* @param json The JSON representation of the polygon
* @return The signed area of the polygon.
*/
private static double getPolygonArea(JSONArray json) throws JSONException {
/// @note this function requires a non-null argument
if (json == null) {
throw new IllegalArgumentException("JsonRegion() requires a non-null GeoJSON argument.");
}
double polygon_area = 0.0;
// @see http://www.faqs.org/faqs/graphics/algorithms-faq/ - "Subject 2.01: How do I find the area of a polygon?"
// @note we skip the last point, which should be a repeat.
// @todo verify that the last point is repeated and skip it conditionally.
int numPoints = json.length() - 1;
for (int i = 0; i < numPoints; i++) {
// first point
double x0 = json.getJSONArray(i).getDouble(0);
double y0 = json.getJSONArray(i).getDouble(1);
double x1 = json.getJSONArray(i + 1).getDouble(0);
double y1 = json.getJSONArray(i + 1).getDouble(1);
polygon_area += (x0 * y1 - x1 * y0) / 2;
}
return polygon_area;
}
/**
* Get the compass direction from this field_region to another field_region.
* Typically, this region is a wind barrier, and the argument region is a
* field boundary.
*
* @param field_region
* @return Character in set {N, E, W, S}.
*/
public String getCompassDir(JsonRegion field_region) {
return getCompassDir(field_region.center_x, field_region.center_y, field_region.region_angle);
}
/**
* Get the compass direction of this region, relative to a given center and
* angular offset. This takes into account the rotation of the rectangle for
* this region.
*
* @param origin_x Origin x coordinate of region to get direction in
* relation to (longitude, WGS84 degrees).
* @param origin_y Origin y coordinate of region to get direction in
* relation to (latitude, WGS84 degrees).
* @param origin_angle Angular rogation of region to get direction in
* relation to (degrees).
* @return Character in set {N, E, W, S}.
*/
public String getCompassDir(double origin_x, double origin_y, double origin_angle) {
// Create vector which is the offset of this region from this origin.
double delta_x = center_x - origin_x;
double delta_y = center_y - origin_y;
if (delta_x == 0.0 || delta_y == 0.0) {
throw new IllegalArgumentException("getCompassDir() requires the origin be different from the region center.");
}
// Get the angle of this delta vector relative to the (rotated) origin.
double vector_angle = Math.atan2(delta_y, delta_x) * (180 / Math.PI) - 90;
double offset_angle = vector_angle - origin_angle;
if (offset_angle < 0) {
offset_angle += 360;
}
if (offset_angle > 360) {
offset_angle -= 360;
}
// Get the compass direction.
// (-45, +45] -> N -45 == 315
// (+45, 135] -> E
// (135, 225] -> S
// (225, 315] -> W
// First rotate by 45 degrees.
// Then it is a simple matter of figuring out what quadrant we are in.
// @todo fix or document the open-closed rounding - calculation doesn't match.
double quandrant_angle = offset_angle + 45;
if (quandrant_angle > 360) {
quandrant_angle -= 360;
}
int quadarant_index = (int) (quandrant_angle / 90.0);
assert (quadarant_index >= 0 && quadarant_index < 4) : "illegal quadrant index";
String compass_dirs[] = {WEPS_BARRIER_NORTH, WEPS_BARRIER_EAST,
WEPS_BARRIER_SOUTH, WEPS_BARRIER_WEST};
String compass_char = compass_dirs[quadarant_index];
return compass_char;
}
// Floating point value used as a tolerance for closeness to zero.
// These algorithms are unstable for thin or zero areas.
static final double double_epsilon = 1e-6;
/**
* Set the representative rectangle by calculating the second moments of the
* area. This produces characteristic vectors which are used to orient the
* characteristic rectangle.
*
* @param json The JSON representation of the polygon
*
* @see http://en.wikipedia.org/wiki/Second_moment_of_area
*/
private void setRepresentativeRectangle(JSONArray json) throws JSONException {
/// @note this function requires a non-null argument
if (json == null) {
throw new IllegalArgumentException("setRepresentativeRectangle() requires a non-null GeoJSON argument.");
}
// don't calculate a representative rectangle for a "Line" geometry type.
if (region_area == 0.0) {
return;
}
// The components of the moment of inertia tensor (2x2 matrix).
double J_xx = 0.0, J_xy = 0.0, J_yx = 0.0, J_yy = 0.0;
// double polygon_area = 0.0;
// Calculate the second moment.
// @note we skip the last point, which should be a repeat.
// @todo verify that the last point is repeated and skip it conditionally.
int numPoints = json.length() - 1;
for (int i = 0; i < numPoints; i++) {
// first point
double x0 = json.getJSONArray(i).getDouble(0) - center_x;
double y0 = json.getJSONArray(i).getDouble(1) - center_y;
double x1 = json.getJSONArray(i + 1).getDouble(0) - center_x;
double y1 = json.getJSONArray(i + 1).getDouble(1) - center_y;
// Calculate contribution of the triangle formed by points [[x0,y0], [x1,y1], [0,0]].
// a_i = twice the area of this triangle.
// @note the tensor is symmetric so we don't explicitly calculate J_yx.
double ai = (x0 * y1 - x1 * y1);
J_xx += (1.0 / 12) * (y0 * y0 + y0 * y1 + y1 * y1) * ai;
J_yy += (1.0 / 12) * (x0 * x0 + x0 * x1 + x1 * x1) * ai;
//
J_xy += (1.0 / 24) * (x0 * y1 + 2 * x0 * y0 + 2 * x1 + y1 + x1 * y0) * ai;
J_yx = J_xy;
}
double area_sign = Math.signum(region_area);
LOG.info ("area_sign=" + area_sign);
// Fix problems with negatively oriented polygon.
// @todo decide orientation based on signs of J_xx, J_yy.
// This allows fixing the region_area sign before this point.
if (area_sign < 0) {
LOG.warning("negatively oriented polygon !!!");
J_xx = -J_xx;
J_yy = -J_yy;
J_xy = -J_xy;
J_yx = -J_yx;
}
// Test for whether the moment of inertia is rotationally symmetric.
// In this case eigenvalues will be the same, and we should assign it an orientation angle of zero.
// This will happen if-and-only-if the discriminant is zero.
// @note We must first convert from degrees to feet or a similar unit.
// 1.0 [degree^2 (lat x long) = 1.0205193309696E+11 [feet^2] (depends on latitude).
// We multiply by the square because the determinant uses squares - its units are [degree^4].
LOG.info("region_area=" + region_area);
double D = J_xx * J_yy - J_xy * J_yx;
LOG.info("THE D VALUE=" + D);
LOG.info("THE D VALUE * 1E+22 vs. double_epsilon=" + (D * 1E+22) + " vs double_epsilon=" + double_epsilon);
if (Math.abs(D) * 1E+22 < double_epsilon) {
region_angle = 0.0;
double area_root = Math.sqrt(Math.abs(region_area));
region_width = area_root;
region_height = area_root;
LOG.info("region width and height set to area_root=" + area_root);
} else {
// We have to solve for the eigenvectors.
// This can be done in closed form for a 2x2 symmetric matrix.
// See e.g. "Applications of Numerical Techniques With C" by Suresh Chandra.
// Note: this is my notation, not from a book or online reference.
double a = J_xx, b = J_xy, c = J_yy;
double delta = c - a;
double trace = c + a;
double R = Math.sqrt(delta * delta + 4.0 * b * b);
// Now the terms of the denominator
double D1 = 8.0 * b * b + 2.0 * delta * delta;
double D2 = 2.0 * delta * R;
// Now the eigenvectors
double J1x = 2.0 * b / (D1 + D2);
double J1y = (delta + R) / (D1 + D2);
double J2x = 2.0 * b / (D1 - D2);
double J2y = (delta - R) / (D1 - D2);
// eigenvalues
double lambda1 = (trace + R) / 2.0;
double lambda2 = (trace - R) / 2.0;
assert (lambda1 > 0 && lambda2 > 0) : "eigenvalues must be positive in setRepresentativeRectangle().";
LOG.info("eigenvalues must be positive in setRepresentativeRectangle().");
LOG.info("lambda1=" + lambda1);
LOG.info("lambda2=" + lambda2);
if (lambda1 <= 0)
{
LOG.warning("lambda1 eigenvalue is not positive!!! -- convert to postive???");
//lambda1 = -lambda1;
}
if (lambda2 <= 0)
{
LOG.warning("lambda2 eigenvalue is not positive!!! -- convert to positive???");
//lambda2 = -lambda2;
}
// Next find the eigenvector which is most northward.
// We first normalize the vectors
double J1 = Math.sqrt(J1x * J1x + J1y * J1y);
double u1x = J1x / J1;
double u1y = J1y / J1;
double J2 = Math.sqrt(J2x * J2x + J2y * J2y);
double u2x = J2x / J2;
double u2y = J2y / J2;
// Verify vector normalization.
assert (Math.abs(u1x * u1x + u1y * u1y - 1.0) < double_epsilon);
assert (Math.abs(u2x * u2x + u2y * u2y - 1.0) < double_epsilon);
/////////////////////////////////////////////////////
// Find the angles of the eigenvectors from "north"
double angle1 = Math.asin(u1x) * (180.0 / Math.PI); // [degrees]
double angle2 = Math.asin(u2x) * (180.0 / Math.PI); // [degrees]
// Reverse the direction of the eigenvector if it is more than 90 degrees.
if (angle1 > +90) {
angle1 -= 180;
}
if (angle1 < -90) {
angle1 += 180;
}
if (angle2 > +90) {
angle2 -= 180;
}
if (angle2 < -90) {
angle2 += 180;
}
//
assert (Math.abs(angle1) <= 90 && Math.abs(angle2) <= 90);
/////////////////////////////////////////////////////
// Find multipliers for each eigenvector to convert from
// moment of inertia to area.
// We need the constraints H*W = A and m1*m2 = 1 to be satisfied.
// Then we can assign H = m1 * sqrt(A), W = m2 * sqrt(A), or the
// reverse if the opposite orientation holds.
LOG.info("lambda1=" + lambda1);
LOG.info("lambda2=" + lambda2);
double mult1 = Math.sqrt(lambda1 / lambda2);
double mult2 = Math.sqrt(lambda2 / lambda1);
//
assert (Math.abs(mult1 / mult2 - lambda1 / lambda2) < double_epsilon);
/////////////////////////////////////////////////////
LOG.info("angle1=" + angle1);
LOG.info("angle2=" + angle2);
LOG.info("region_area=" + region_area);
LOG.info("mult1=" + mult1);
LOG.info("mult2=" + mult2);
// Now identify the most "nortward" eigenvector.
// Use the most northward eigenvector mult for the "height".
// Use the other eigenvector mult for the "width".
if (Math.abs(angle1) < Math.abs(angle2)) {
// Eivenvector J1 is most northward - use it for the "height"
region_angle = angle1;
region_height = mult1 * Math.sqrt(Math.abs(region_area));
region_width = mult2 * Math.sqrt(Math.abs(region_area));
LOG.info("Eigenvector J1 is most northward region_height=" + region_height);
LOG.info("Eigenvector J1 is most northward region_width=" + region_width);
} else {
// Eivenvector J2 is most northward - use it for the "height"
region_angle = angle2;
region_height = mult2 * Math.sqrt(Math.abs(region_area));
region_width = mult1 * Math.sqrt(Math.abs(region_area));
LOG.info("Eigenvector J2 is most northward region_height=" + region_height);
LOG.info("Eigenvector J2 is most northward region_width=" + region_width);
}
}
LOG.info("the region_height=" + region_height);
LOG.info("the region_width=" + region_width);
}
/**
* @param args the command line arguments
* @throws IOException
* @throws JSONException
*/
public static void main(String[] args) throws IOException, JSONException {
// Load a field GeoJSON file.
File field_json_file = new File("test/field.json");
String field_json_in = FileUtils.readFileToString(field_json_file, "UTF-8");
JSONTokener field_json_tok = new JSONTokener(field_json_in);
JSONObject field_json_object = new JSONObject(field_json_tok);
JsonRegion field_region = new JsonRegion(field_json_object);
// Load a barrier GeoJSON file.
File barrier_json_file = new File("test/barrier.json");
String barrier_json_in = FileUtils.readFileToString(barrier_json_file, "UTF-8");
JSONTokener barrier_json_tok = new JSONTokener(barrier_json_in);
JSONObject barrier_json_object = new JSONObject(barrier_json_tok);
JsonRegion barrier_region = new JsonRegion(barrier_json_object);
// Find the compass direction of the barrier in relation to the field.
String barrier_dir = barrier_region.getCompassDir(field_region);
// Log information about the field.
System.out.println("== Field information ==");
System.out.println("Area: " + field_region.getArea() + " [meter^2]");
System.out.println("Angle: " + field_region.region_angle + " [degrees]");
System.out.println("Height: " + field_region.getHeight() + " [meters]");
System.out.println("Width: " + field_region.getWidth() + " [meters]");
System.out.println("Center: (" + field_region.center_x + ", " + field_region.center_y + ") [coords WGS84]");
System.out.println("");
// Log information about the barrier.
System.out.println("== Barrier information ==");
System.out.println("Area: " + barrier_region.getArea() + " [meter^2]");
System.out.println("Angle: " + barrier_region.region_angle + " [degrees]");
System.out.println("Height: " + barrier_region.getHeight() + " [meters]");
System.out.println("Width: " + barrier_region.getWidth() + " [meters]");
System.out.println("Center: (" + barrier_region.center_x + ", " + barrier_region.center_y + ") [coords WGS84]");
System.out.println("Direct: " + barrier_dir);
}
}