327 lines
11 KiB
JavaScript
327 lines
11 KiB
JavaScript
const proj4 = require('proj4'),
|
||
utils = require('./utils'),
|
||
turf = require('@turf/turf'),
|
||
{ Errors } = require('./constants'),
|
||
{ AppInputError } = require('./app_error');
|
||
|
||
//'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]';
|
||
//+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
|
||
|
||
proj4.defs([
|
||
[
|
||
'EPSG:4326',
|
||
'+title=WGS 84 (long/lat) +proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees'
|
||
]
|
||
]);
|
||
|
||
function UTMtoLnLats(fromPrj, lnglats, reverse = false) {
|
||
if (!fromPrj || utils.isEmptyArray(lnglats)) return null;
|
||
|
||
return lnglats.map(co => {
|
||
const res = proj4(fromPrj, proj4('EPSG:4326'), reverse ? co.reverse() : co);
|
||
return [utils.truncR(res[0], 6), utils.truncR(res[1], 6)];
|
||
});
|
||
}
|
||
|
||
/**
|
||
* Return UTM zone number given a Central Meridian
|
||
* @param {number} cm the Central Meridian
|
||
* @returns zone number
|
||
*/
|
||
function zoneFromCM(cm) {
|
||
let zoneNo = 0;
|
||
if (cm >= 3 && cm <= 177)
|
||
zoneNo = 31 + Math.floor(Math.round(cm - 3) / 6);
|
||
else if (cm >= 183 && cm <= 357)
|
||
zoneNo = 1 + Math.floor(Math.round(cm - 183) / 6);
|
||
|
||
return zoneNo;
|
||
}
|
||
|
||
/**
|
||
* Return Central Meridian given a UTM zone
|
||
* @param {*} zone UTM zone number
|
||
* @returns Central Meridian
|
||
*/
|
||
function cmFromZone(zone) {
|
||
let cm = 0;
|
||
if (zone > 30)
|
||
cm = (zone * 6) - 183;
|
||
else
|
||
cm = (zone * 6) + 177;
|
||
return cm;
|
||
}
|
||
|
||
/**
|
||
* Convert to -180..180 range
|
||
* param value number in 0..359 range
|
||
*/
|
||
function to180Range(value) {
|
||
// let newVal = value;
|
||
// // convert heading to range of [-180,180] degrees
|
||
// while (newVal > 180)
|
||
// newVal -= 360;
|
||
// while (newVal < -180)
|
||
// newVal += 360;
|
||
// return newVal;
|
||
return utils.isNumber(value) ? ((value + 540) % 360 - 180) : 0.0;
|
||
}
|
||
|
||
/**
|
||
* Convert to 0..359 range
|
||
* param value number in -180..180 range
|
||
*/
|
||
function to360Range(value) {
|
||
// let newVal = value;
|
||
// while (newVal < 0)
|
||
// newVal += 360;
|
||
// while (newVal >= 360)
|
||
// newVal -= 360;
|
||
|
||
// return newVal;
|
||
return ((value + 360) % 360);
|
||
}
|
||
|
||
/**
|
||
* Check whether a point inside or on the edge of a polygon using Ray Casting algorithm
|
||
* When a point is given, then we virtually draw a line from a point far away outside from the polygon to the given point.
|
||
* Then, we calculate the number of intersections of the virtual line with the edges of the polygon.
|
||
* If the number of intersections is odd, then according to the Ray Casting theory we can conclude that the point is inside the polygon.
|
||
* Otherwise, the point is outside the polygon.
|
||
* @param {*} polyPoints vertexes list
|
||
* @param {*} x x coordinate
|
||
* @param {*} y y coordinate
|
||
*/
|
||
function pointInOrOnPolygon(polyPoints, x, y) {
|
||
let intersections = 0;
|
||
|
||
for (let i = 0, j = polyPoints.length - 1; i < polyPoints.length; j = i++) {
|
||
let xi = polyPoints[i].x, yi = polyPoints[i].y;
|
||
let xj = polyPoints[j].x, yj = polyPoints[j].y;
|
||
|
||
if (Math.abs(yj - yi) <= 1e-4 && Math.abs(yj - y) <= 1e-4 && x >= Math.min(xj, xi) && x <= Math.max(xj, xi)) { // Check if point is on an horizontal polygon boundary
|
||
return true;
|
||
}
|
||
if (y >= Math.min(yj, yi) && y <= Math.max(yj, yi) && x <= Math.max(xj, xi) && yj != yi) {
|
||
const ss = (y - yj) * (xi - xj) / (yi - yj) + xj;
|
||
if (Math.abs(ss - x) <= 1e-4) { // Check if point is on the polygon boundary (other than horizontal)
|
||
return true;
|
||
}
|
||
if (Math.abs(xj - xi) <= 1e-4 || x <= ss) {
|
||
intersections++;
|
||
}
|
||
}
|
||
}
|
||
// If the number of edges we passed through is odd, then it’s in the polygon.
|
||
if (intersections % 2 != 0) {
|
||
return true;
|
||
} else {
|
||
return false;
|
||
}
|
||
}
|
||
|
||
function isValidLL(lat, lon) {
|
||
if (undefined === lat || undefined === lon) return false;
|
||
return (-90 <= lat && lat <= 90) && (-180 <= lon && lon <= 180);
|
||
}
|
||
|
||
function isValidLLAr(latlonAr, isLonLat = false) {
|
||
if (utils.isEmptyArray(latlonAr) || latlonAr.length < 2) return false;
|
||
return isLonLat ? isValidLL(latlonAr[1], latlonAr[0]) : isValidLL(latlonAr[0], latlonAr[1]);
|
||
}
|
||
|
||
/**
|
||
* Returns Central Meridian given a longitude
|
||
* @param {number} longitude
|
||
* @returns {number} Central Meridian
|
||
*/
|
||
function computeCMfromLon(longitude) {
|
||
return Math.trunc((longitude < 0 ? 360 + longitude : longitude) / 6) * 6 + 3;
|
||
}
|
||
|
||
/**
|
||
* Returns UTM Zone number given a Central Meridian
|
||
* @param {number} CM Central Meridian
|
||
* @returns {number} Zone number
|
||
*/
|
||
function computeZoneNofromCM(CM) {
|
||
let zoneno = 0;
|
||
if (CM >= 3 && CM <= 177)
|
||
zoneno = 31 + Math.round(CM - 3) / 6;
|
||
else if (CM >= 183 && CM <= 357)
|
||
zoneno = 1 + Math.round(CM - 183) / 6;
|
||
|
||
if (zoneno > 60)
|
||
zoneno -= 60;
|
||
return zoneno;
|
||
}
|
||
|
||
/**
|
||
* Returns UTM zone number given a longitude
|
||
* @param {number} longitude
|
||
* @returns {number} zone number
|
||
*/
|
||
function zoneNoFromLon(longitude) {
|
||
if (!utils.isNumber(longitude))
|
||
AppInputError.throw(Errors.INVALID_LONGITUDE);
|
||
|
||
return Math.round((longitude + 180) / 6);
|
||
}
|
||
|
||
/**
|
||
* Distance between 2 latlong points in km using Haversine formular
|
||
* @param {number[]} latLng fisrt point [lat, lon]
|
||
* @param {number[]} latLng2 second point [lat, lon]
|
||
*/
|
||
function distance(latLng, latLng2) {
|
||
if (!(latLng && latLng.length) || !(latLng2 && latLng2.length))
|
||
return 0;
|
||
|
||
return turf.distance(latLng.slice().reverse(), latLng2.slice().reverse());
|
||
}
|
||
|
||
function getUTMProj(coorType, equo_crossing, zoneCM) {
|
||
const zone = coorType === 'U' ? zoneFromCM(zoneCM) : zoneCM;
|
||
const hemisphere = equo_crossing <= 1 ? 'N' : 'S';
|
||
let cmEW = coorType === 'U' ? zoneCM : cmFromZone(zone);
|
||
cmEW = to180Range(cmEW);
|
||
|
||
const falseNorthing = hemisphere === 'N' ? 0 : 10000000;
|
||
// AgNav always use WGS-84 geographic CR
|
||
const fromPrj = `PROJCS["WGS_1984_UTM_Zone_${zone}${hemisphere}",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",${cmEW}],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",${falseNorthing}],UNIT["Meter",1]]`;
|
||
return fromPrj;
|
||
}
|
||
|
||
/**
|
||
* Returns the destination point from ‘this’ point having travelled the given distance on the
|
||
* given initial bearing (bearing normally varies around path followed).
|
||
* Extracted from: https://github.com/chrisveness/geodesy/blob/master/latlon-spherical.js
|
||
* @param {number} distance - Distance travelled, in same units as earth radius (default: metres).
|
||
* @param {number} bearing - Initial bearing in degrees from north.
|
||
* @param {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
|
||
* @returns {number[]} Destination point in [Lat, Lon].
|
||
*
|
||
* @example
|
||
* const p1 = new LatLon(51.47788, -0.00147);
|
||
* const p2 = p1.destinationPoint(7794, 300.7); // 51.5136°N, 000.0983°W
|
||
*/
|
||
function destinationPoint(lat, lon, distance, bearing, radius = 6371e3) {
|
||
// sinφ2 = sinφ1⋅cosδ + cosφ1⋅sinδ⋅cosθ
|
||
// tanΔλ = sinθ⋅sinδ⋅cosφ1 / cosδ−sinφ1⋅sinφ2
|
||
// see mathforum.org/library/drmath/view/52049.html for derivation
|
||
|
||
const δ = distance / radius; // angular distance in radians
|
||
const θ = turf.degreesToRadians(bearing);
|
||
|
||
const φ1 = turf.degreesToRadians(lat), λ1 = turf.degreesToRadians(lon);
|
||
|
||
const sinφ2 = Math.sin(φ1) * Math.cos(δ) + Math.cos(φ1) * Math.sin(δ) * Math.cos(θ);
|
||
const φ2 = Math.asin(sinφ2);
|
||
const y = Math.sin(θ) * Math.sin(δ) * Math.cos(φ1);
|
||
const x = Math.cos(δ) - Math.sin(φ1) * sinφ2;
|
||
const λ2 = λ1 + Math.atan2(y, x);
|
||
|
||
const dlat = turf.radiansToDegrees(φ2);
|
||
const dlon = turf.radiansToDegrees(λ2);
|
||
|
||
return [dlat, dlon];
|
||
}
|
||
|
||
|
||
/**
|
||
* Returns the destination point from this point having travelled the given initial bearing, speed and time
|
||
* @param {*} lat Current Latitude
|
||
* @param {*} lon Current Longitude
|
||
* @param {*} grSpeed Ground speed in m/s
|
||
* @param {*} heading Initial bearing in degrees from north
|
||
* @param {*} timeAdv Projected time in seconds
|
||
* @returns {*} { lat, lon }
|
||
*/
|
||
function projectLatLong(lat, lon, grSpeed, heading, timeAdv) {
|
||
const destPoint = destinationPoint(lat, lon, (grSpeed * timeAdv), heading);
|
||
return ({ lat: destPoint[0], lon: destPoint[1] });
|
||
}
|
||
|
||
/**
|
||
* Return the Bearing/Heading [0..360] of the line of two given points and the true North;
|
||
* in Mecarter projection such as UTM
|
||
* @param {*} x1 easting of point 1
|
||
* @param {*} y1 northing of point 1
|
||
* @param {*} x2 easting of point 2
|
||
* @param {*} y2 northing of point 1
|
||
*/
|
||
function utmBearing(x1, y1, x2, y2) {
|
||
return to360Range(Math.atan2(x2 - x1, y2 - y1) * 180 / Math.PI);
|
||
}
|
||
|
||
/**
|
||
* Returns overall bounding box and mutates areas, if any, each with a bounding box for 'bbox' prop
|
||
* @param {*} areas job areas list
|
||
* @returns bbox extent in [minX, minY, maxX, maxY]
|
||
*/
|
||
function updateAreasBBoxLL(areas) {
|
||
if (utils.isEmptyArray(areas)) return areas;
|
||
|
||
const bbox = [Number.MAX_VALUE, Number.MAX_VALUE, (-1 * Number.MAX_VALUE), (-1 * Number.MAX_VALUE)]; // minX, minY, maxX, maxY
|
||
let areaBbox;
|
||
|
||
for (let i = 0; i < areas.length; i++) {
|
||
areaBbox = areas[i]['bbox'];
|
||
if (!areaBbox)
|
||
areaBbox = areas[i]['bbox'] = turf.bbox(areas[i].geometry);
|
||
|
||
if (areaBbox[0] < bbox[0]) bbox[0] = areaBbox[0];
|
||
if (areaBbox[1] < bbox[1]) bbox[1] = areaBbox[1];
|
||
if (areaBbox[2] > bbox[2]) bbox[2] = areaBbox[2];
|
||
if (areaBbox[3] > bbox[3]) bbox[3] = areaBbox[3];
|
||
}
|
||
|
||
return bbox;
|
||
}
|
||
|
||
/**
|
||
* Returns overall bounding box (areas with UTM coordinates)
|
||
* @param {*} areas job areas list
|
||
* @returns bbox extent in [minX, minY, maxX, maxY]
|
||
*/
|
||
function getAreasBBox(areas) {
|
||
if (utils.isEmptyArray(areas)) return areas;
|
||
|
||
const bbox = [Number.MAX_VALUE, Number.MAX_VALUE, (-1 * Number.MAX_VALUE), (-1 * Number.MIN_VALUE)]; // minX, minY, maxX, maxY
|
||
let areaBbox;
|
||
|
||
for (let i = 0; i < areas.length; i++) {
|
||
areaBbox = areas[i].bbox;
|
||
if (utils.isEmptyArray(areaBbox)) continue;
|
||
|
||
if (areaBbox[0] < bbox[0]) bbox[0] = areaBbox[0];
|
||
if (areaBbox[1] < bbox[1]) bbox[1] = areaBbox[1];
|
||
if (areaBbox[2] > bbox[2]) bbox[2] = areaBbox[2];
|
||
if (areaBbox[3] > bbox[3]) bbox[3] = areaBbox[3];
|
||
}
|
||
|
||
return bbox;
|
||
}
|
||
|
||
/**
|
||
* Returns the prefered UTM zone and the hemisphere given the LL bounding box
|
||
* @param {number[]} bbox bbox extent in [minX, minY, maxX, maxY]
|
||
* @returns prefered zone info { zone: [zone_number], hemisphere: ['N'|'S'] }
|
||
*/
|
||
function calcRefZonebyBbox(bbox) {
|
||
let refZone = { zone: undefined, hemisphere: undefined };
|
||
if (utils.isEmptyArray(bbox) || bbox.length < 4)
|
||
return refZone;
|
||
|
||
let centerX = (bbox[2] - bbox[0]) / 2, centerY = (bbox[3] - bbox[1]) / 2;
|
||
refZone.zone = zoneNoFromLon((bbox[0] + centerX));
|
||
refZone.hemisphere = (bbox[1] + centerY) > 0 ? 'N' : 'S';
|
||
|
||
return refZone;
|
||
}
|
||
|
||
module.exports = {
|
||
UTMtoLnLats, zoneFromCM, cmFromZone, to180Range, to360Range, pointInOrOnPolygon, isValidLL, isValidLLAr, computeCMfromLon, computeZoneNofromCM, zoneNoFromLon, distance,
|
||
getUTMProj, projectLatLong, utmBearing, updateAreasBBoxLL, calcRefZonebyBbox, getAreasBBox
|
||
}
|