/*
* PolygonArea.js
* Transcription of PolygonArea.[ch]pp into JavaScript.
*
* See the documentation for the C++ class. The conversion is a literal
* conversion from C++.
*
* The algorithms are derived in
*
* Charles F. F. Karney,
* Algorithms for geodesics, J. Geodesy 87, 43-55 (2013);
* https://doi.org/10.1007/s00190-012-0578-z
* Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
*
* Copyright (c) Charles Karney (2011-2022) <karney@alum.mit.edu> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*/
// Load AFTER geodesic/Math.js and geodesic/Geodesic.js
(function(
/**
* @exports geodesic/PolygonArea
* @description Compute the area of geodesic polygons via the
* {@link module:geodesic/PolygonArea.PolygonArea PolygonArea}
* class.
*/
p, g, m, a) {
"use strict";
var transit, transitdirect, AreaReduceA, AreaReduceB;
transit = function(lon1, lon2) {
// Return 1 or -1 if crossing prime meridian in east or west direction.
// Otherwise return zero. longitude = +/-0 considered to be positive.
// This is (should be?) compatible with transitdirect which computes
// exactly the parity of
// int(floor((lon1 + lon12) / 360)) - int(floor(lon1 / 360)))
var lon12 = m.AngDiff(lon1, lon2).d;
lon1 = m.AngNormalize(lon1);
lon2 = m.AngNormalize(lon2);
// edge case lon1 = 180, lon2 = 360->0, lon12 = 180 to give 1
return lon12 > 0 &&
// lon12 > 0 && lon1 > 0 && lon2 == 0 implies lon1 == 180
((lon1 < 0 && lon2 >= 0) || (lon1 > 0 && lon2 === 0)) ? 1 :
// non edge case lon1 = -180, lon2 = -360->-0, lon12 = -180
(lon12 < 0 && lon1 >= 0 && lon2 < 0 ? -1 : 0);
};
// an alternate version of transit to deal with longitudes in the direct
// problem.
transitdirect = function(lon1, lon2) {
// We want to compute exactly
// int(floor(lon2 / 360)) - int(floor(lon1 / 360))
lon1 = lon1 % 720; lon2 = lon2 % 720;
return ((0 <= lon2 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
((0 <= lon1 && lon1 < 360) || lon1 < -360 ? 0 : 1 );
};
// Reduce Accumulator area
AreaReduceA = function(area, area0, crossings, reverse, sign) {
area.Remainder(area0);
if (crossings & 1)
area.Add( (area.Sum() < 0 ? 1 : -1) * area0/2 );
// area is with the clockwise sense. If !reverse convert to
// counter-clockwise convention.
if (!reverse)
area.Negate();
// If sign put area in (-area0/2, area0/2], else put area in [0, area0)
if (sign) {
if (area.Sum() > area0/2)
area.Add( -area0 );
else if (area.Sum() <= -area0/2)
area.Add( +area0 );
} else {
if (area.Sum() >= area0)
area.Add( -area0 );
else if (area.Sum() < 0)
area.Add( +area0 );
}
return 0 + area.Sum();
};
// Reduce double area
AreaReduceB = function(area, area0, crossings, reverse, sign) {
area = m.remainder(area, area0);
if (crossings & 1)
area += (area < 0 ? 1 : -1) * area0/2;
// area is with the clockwise sense. If !reverse convert to
// counter-clockwise convention.
if (!reverse)
area *= -1;
// If sign put area in (-area0/2, area0/2], else put area in [0, area0)
if (sign) {
if (area > area0/2)
area -= area0;
else if (area <= -area0/2)
area += area0;
} else {
if (area >= area0)
area -= area0;
else if (area < 0)
area += area0;
}
return 0 + area;
};
/**
* @class
* @property {number} a the equatorial radius (meters).
* @property {number} f the flattening.
* @property {bool} polyline whether the PolygonArea object describes a
* polyline or a polygon.
* @property {number} num the number of vertices so far.
* @property {number} lat the current latitude (degrees).
* @property {number} lon the current longitude (degrees).
* @summary Initialize a PolygonArea object.
* @classdesc Computes the area and perimeter of a geodesic polygon.
* This object is usually instantiated by
* {@link module:geodesic/Geodesic.Geodesic#Polygon Geodesic.Polygon}.
* @param {object} geod a {@link module:geodesic/Geodesic.Geodesic
* Geodesic} object.
* @param {bool} [polyline = false] if true the new PolygonArea object
* describes a polyline instead of a polygon.
*/
p.PolygonArea = function(geod, polyline) {
this._geod = geod;
this.a = this._geod.a;
this.f = this._geod.f;
this._area0 = 4 * Math.PI * geod._c2;
this.polyline = !polyline ? false : polyline;
this._mask = g.LATITUDE | g.LONGITUDE | g.DISTANCE |
(this.polyline ? g.NONE : g.AREA | g.LONG_UNROLL);
if (!this.polyline)
this._areasum = new a.Accumulator(0);
this._perimetersum = new a.Accumulator(0);
this.Clear();
};
/**
* @summary Clear the PolygonArea object, setting the number of vertices to
* 0.
*/
p.PolygonArea.prototype.Clear = function() {
this.num = 0;
this._crossings = 0;
if (!this.polyline)
this._areasum.Set(0);
this._perimetersum.Set(0);
this._lat0 = this._lon0 = this.lat = this.lon = NaN;
};
/**
* @summary Add the next vertex to the polygon.
* @param {number} lat the latitude of the point (degrees).
* @param {number} lon the longitude of the point (degrees).
* @description This adds an edge from the current vertex to the new vertex.
*/
p.PolygonArea.prototype.AddPoint = function(lat, lon) {
var t;
if (this.num === 0) {
this._lat0 = this.lat = lat;
this._lon0 = this.lon = lon;
} else {
t = this._geod.Inverse(this.lat, this.lon, lat, lon, this._mask);
this._perimetersum.Add(t.s12);
if (!this.polyline) {
this._areasum.Add(t.S12);
this._crossings += transit(this.lon, lon);
}
this.lat = lat;
this.lon = lon;
}
++this.num;
};
/**
* @summary Add the next edge to the polygon.
* @param {number} azi the azimuth at the current the point (degrees).
* @param {number} s the length of the edge (meters).
* @description This specifies the new vertex in terms of the edge from the
* current vertex.
*/
p.PolygonArea.prototype.AddEdge = function(azi, s) {
var t;
if (this.num) {
t = this._geod.Direct(this.lat, this.lon, azi, s, this._mask);
this._perimetersum.Add(s);
if (!this.polyline) {
this._areasum.Add(t.S12);
this._crossings += transitdirect(this.lon, t.lon2);
}
this.lat = t.lat2;
this.lon = t.lon2;
}
++this.num;
};
/**
* @summary Compute the perimeter and area of the polygon.
* @param {bool} reverse if true then clockwise (instead of
* counter-clockwise) traversal counts as a positive area.
* @param {bool} sign if true then return a signed result for the area if the
* polygon is traversed in the "wrong" direction instead of returning the
* area for the rest of the earth.
* @return {object} r where r.number is the number of vertices, r.perimeter
* is the perimeter (meters), and r.area (only returned if polyline is
* false) is the area (meters<sup>2</sup>).
* @description Arbitrarily complex polygons are allowed. In the case of
* self-intersecting polygons the area is accumulated "algebraically",
* e.g., the areas of the 2 loops in a figure-8 polygon will partially
* cancel. If the object is a polygon (and not a polyline), the perimeter
* includes the length of a final edge connecting the current point to the
* initial point. If the object is a polyline, then area is nan. More
* points can be added to the polygon after this call.
*/
p.PolygonArea.prototype.Compute = function(reverse, sign) {
var vals = {number: this.num}, t, tempsum;
if (this.num < 2) {
vals.perimeter = 0;
if (!this.polyline)
vals.area = 0;
return vals;
}
if (this.polyline) {
vals.perimeter = this._perimetersum.Sum();
return vals;
}
t = this._geod.Inverse(this.lat, this.lon, this._lat0, this._lon0,
this._mask);
vals.perimeter = this._perimetersum.Sum(t.s12);
tempsum = new a.Accumulator(this._areasum);
tempsum.Add(t.S12);
vals.area = AreaReduceA(tempsum, this._area0,
this._crossings + transit(this.lon, this._lon0),
reverse, sign);
return vals;
};
/**
* @summary Compute the perimeter and area of the polygon with a tentative
* new vertex.
* @param {number} lat the latitude of the point (degrees).
* @param {number} lon the longitude of the point (degrees).
* @param {bool} reverse if true then clockwise (instead of
* counter-clockwise) traversal counts as a positive area.
* @param {bool} sign if true then return a signed result for the area if the
* polygon is traversed in the "wrong" direction instead of returning the
* area for the rest of the earth.
* @return {object} r where r.number is the number of vertices, r.perimeter
* is the perimeter (meters), and r.area (only returned if polyline is
* false) is the area (meters<sup>2</sup>).
* @description A new vertex is *not* added to the polygon.
*/
p.PolygonArea.prototype.TestPoint = function(lat, lon, reverse, sign) {
var vals = {number: this.num + 1}, t, tempsum, crossings, i;
if (this.num === 0) {
vals.perimeter = 0;
if (!this.polyline)
vals.area = 0;
return vals;
}
vals.perimeter = this._perimetersum.Sum();
tempsum = this.polyline ? 0 : this._areasum.Sum();
crossings = this._crossings;
for (i = 0; i < (this.polyline ? 1 : 2); ++i) {
t = this._geod.Inverse(
i === 0 ? this.lat : lat, i === 0 ? this.lon : lon,
i !== 0 ? this._lat0 : lat, i !== 0 ? this._lon0 : lon,
this._mask);
vals.perimeter += t.s12;
if (!this.polyline) {
tempsum += t.S12;
crossings += transit(i === 0 ? this.lon : lon,
i !== 0 ? this._lon0 : lon);
}
}
if (this.polyline)
return vals;
vals.area = AreaReduceB(tempsum, this._area0, crossings, reverse, sign);
return vals;
};
/**
* @summary Compute the perimeter and area of the polygon with a tentative
* new edge.
* @param {number} azi the azimuth of the edge (degrees).
* @param {number} s the length of the edge (meters).
* @param {bool} reverse if true then clockwise (instead of
* counter-clockwise) traversal counts as a positive area.
* @param {bool} sign if true then return a signed result for the area if the
* polygon is traversed in the "wrong" direction instead of returning the
* area for the rest of the earth.
* @return {object} r where r.number is the number of vertices, r.perimeter
* is the perimeter (meters), and r.area (only returned if polyline is
* false) is the area (meters<sup>2</sup>).
* @description A new vertex is *not* added to the polygon.
*/
p.PolygonArea.prototype.TestEdge = function(azi, s, reverse, sign) {
var vals = {number: this.num ? this.num + 1 : 0}, t, tempsum, crossings;
if (this.num === 0)
return vals;
vals.perimeter = this._perimetersum.Sum() + s;
if (this.polyline)
return vals;
tempsum = this._areasum.Sum();
crossings = this._crossings;
t = this._geod.Direct(this.lat, this.lon, azi, s, this._mask);
tempsum += t.S12;
crossings += transitdirect(this.lon, t.lon2);
crossings += transit(t.lon2, this._lon0);
t = this._geod.Inverse(t.lat2, t.lon2, this._lat0, this._lon0, this._mask);
vals.perimeter += t.s12;
tempsum += t.S12;
vals.area = AreaReduceB(tempsum, this._area0, crossings, reverse, sign);
return vals;
};
})(geodesic.PolygonArea, geodesic.Geodesic,
geodesic.Math, geodesic.Accumulator);