From 6472d9d9925b448e8cf11405681ca3d573b50dba Mon Sep 17 00:00:00 2001 From: "S.Listl" <S.Listl@SLISTL.aditosoftware.local> Date: Mon, 19 Aug 2019 13:31:34 +0200 Subject: [PATCH] Location_lib with functions for area search --- process/Location_lib/Location_lib.aod | 9 ++ process/Location_lib/process.js | 166 ++++++++++++++++++++++++++ 2 files changed, 175 insertions(+) create mode 100644 process/Location_lib/Location_lib.aod create mode 100644 process/Location_lib/process.js diff --git a/process/Location_lib/Location_lib.aod b/process/Location_lib/Location_lib.aod new file mode 100644 index 0000000000..1a9a229f65 --- /dev/null +++ b/process/Location_lib/Location_lib.aod @@ -0,0 +1,9 @@ +<?xml version="1.0" encoding="UTF-8"?> +<process xmlns="http://www.adito.de/2018/ao/Model" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" VERSION="1.2.1" xsi:schemaLocation="http://www.adito.de/2018/ao/Model adito://models/xsd/process/1.2.1"> + <name>Location_lib</name> + <majorModelMode>DISTRIBUTED</majorModelMode> + <process>%aditoprj%/process/Location_lib/process.js</process> + <variants> + <element>LIBRARY</element> + </variants> +</process> diff --git a/process/Location_lib/process.js b/process/Location_lib/process.js new file mode 100644 index 0000000000..335f688606 --- /dev/null +++ b/process/Location_lib/process.js @@ -0,0 +1,166 @@ +import("Util_lib"); +import("system.eMath"); + +function GeoLocationUtils () {} + +/** + * average earth radius in kilometers + */ +GeoLocationUtils.EARTH_RADIUS_AVG = 6371.001; + +/** + * earth radius at the equator in kilometers + */ +GeoLocationUtils.EARTH_RADIUS_EQUATOR = 6378.137; + +/** + * earth radius at the poles in kilometers + */ +GeoLocationUtils.EARTH_RADIUS_POLE = 6356.752; + +/** + * Calculates the distance between two given locations using the haversine formula. + * This method is reasonably accurate for most applications (in most cases the error is under 0.3%). + * + * @param {String|Number} pLatA latitude of the first location + * @param {String|Number} pLonA longitude of the first location + * @param {String|Number} pLatB latitude of the second location + * @param {String|Number} pLonB longitude of the second location + * @return {Number} the distance in kilometers + */ +GeoLocationUtils.distanceHaversine = function (pLatA, pLonA, pLatB, pLonB) +{ + pLatA = GeoLocationUtils.degreeToRadians(pLatA); + pLonA = GeoLocationUtils.degreeToRadians(pLonA); + pLatB = GeoLocationUtils.degreeToRadians(pLatB); + pLonB = GeoLocationUtils.degreeToRadians(pLonB); + + var R = GeoLocationUtils.EARTH_RADIUS_AVG; //average Earth radius in km https://en.wikipedia.org/wiki/Earth_radius + + var deltaLat = eMath.subDec(pLatA, eMath.absDec(pLatB)); //lat can be negative + var deltaLon = eMath.subDec(pLonA, pLonB); + + var sinLat = GeoLocationUtils.haversine(deltaLat); + var sinLon = GeoLocationUtils.haversine(deltaLon); + + return R * 2 * Math.asin(Math.sqrt(Math.pow(sinLat, 2) + Math.cos(pLatA) * Math.cos(eMath.absDec(pLatB)) * Math.pow(sinLon, 2))); +} + +/** + * JavaScript function to calculate the geodetic distance between two points specified by latitude/longitude using the Vincenty inverse formula for ellipsoids. + * This is more accurate than the haversine formula but also more complex and thus theoretically not as fast. + * + * Taken from http://movable-type.co.uk/scripts/latlong-vincenty.html + * + * @param {Number} latA: first point in decimal degrees + * @param {Number} lonA: first point in decimal degrees + * @param {Number} latB: second point in decimal degrees + * @param {Number} lonB: second point in decimal degrees + * @returns {Number} distance in kilometres between points + */ +GeoLocationUtils.distanceVincenty = function (latA, lonA, latB, lonB) +{ + let phi1 = GeoLocationUtils.degreeToRadians(latA), lambda1 = GeoLocationUtils.degreeToRadians(lonA); + let phi2 = GeoLocationUtils.degreeToRadians(latB), lambda2 = GeoLocationUtils.degreeToRadians(lonB); + + let a_ = 6378137; + let b = 6356752.314245; + let f = 1/298.257223563; + let epsilon = Math.pow(2, -52); + + let L = lambda2 - lambda1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanphi. + let tanU1 = (1-f) * Math.tan(phi1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; + let tanU2 = (1-f) * Math.tan(phi2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2; + + let antipodal = Math.abs(L) > Math.PI/2 || Math.abs(phi2-phi1) > Math.PI/2; + + let maxIterations = 100; + + let lambda = L, sinlambda = null, coslambda = null; // lambda = difference in longitude on an auxiliary sphere + let sigma = antipodal ? Math.PI : 0, sinsigma = 0, cossigma = antipodal ? -1 : 1, sinSqsigma = null; // sigma = angular distance P₠P₂ on the sphere + let cos2sigmaM = 1; // sigmaM = angular distance on the sphere from the equator to the midpoint of the line + let sinalpha = null, cosSqalpha = 1; // alpha = azimuth of the geodesic at the equator + let C = null; + + let lambdaP = null, iterations = 0; + do { + sinlambda = Math.sin(lambda); + coslambda = Math.cos(lambda); + sinSqsigma = (cosU2*sinlambda) * (cosU2*sinlambda) + (cosU1*sinU2-sinU1*cosU2*coslambda) * (cosU1*sinU2-sinU1*cosU2*coslambda); + + if (Math.abs(sinSqsigma) < epsilon) + break; // co-incident/antipodal points (falls back on lambda/sigma = L) + + sinsigma = Math.sqrt(sinSqsigma); + cossigma = sinU1*sinU2 + cosU1*cosU2*coslambda; + sigma = Math.atan2(sinsigma, cossigma); + sinalpha = cosU1 * cosU2 * sinlambda / sinsigma; + cosSqalpha = 1 - sinalpha*sinalpha; + cos2sigmaM = (cosSqalpha != 0) ? (cossigma - 2*sinU1*sinU2/cosSqalpha) : 0; // on equatorial line cos²alpha = 0 (§6) + C = f/16*cosSqalpha*(4+f*(4-3*cosSqalpha)); + lambdaP = lambda; + lambda = L + (1-C) * f * sinalpha * (sigma + C*sinsigma*(cos2sigmaM+C*cossigma*(-1+2*cos2sigmaM*cos2sigmaM))); + + } while (Math.abs(lambda-lambdaP) > 1e-12 && ++iterations<maxIterations); + + if (iterations >= maxIterations) + return NaN; + + let uSq = cosSqalpha * (a_*a_ - b*b) / (b*b); + let A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); + let B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); + let deltasigma = B*sinsigma*(cos2sigmaM+B/4*(cossigma*(-1+2*cos2sigmaM*cos2sigmaM)- + B/6*cos2sigmaM*(-3+4*sinsigma*sinsigma)*(-3+4*cos2sigmaM*cos2sigmaM))); + + let s = b*A*(sigma-deltasigma); // s = length of the geodesic + + return s / 1000; //convert metres into kilometres +} + +/** + * haversine function + * + * @param {String|Number} pDelta req delta-value in radians + * @return {Number} hav(pDelta) + */ +GeoLocationUtils.haversine = function (pDelta) +{ + return Math.sin(eMath.divDec(pDelta, 2)); +} + +/** + * converts degrees to radians + * + * @param {String|Number} pAngle angle in degrees + * @return {Number} angle in radians + */ +GeoLocationUtils.degreeToRadians = function (pAngle) +{ + return pAngle * Math.PI / 180; +} + +function AreaSearchUtils () {} + +/** + * WIP + * @ignore + */ +AreaSearchUtils.findOrganisationsInArea = function (pLatitude, pLongitude, pMaxDistance) +{ + var locationData = []; + var locations = []; + + for (let i = 0, l = locationData.length; i < l; i++) + { + let distance = GeoLocationUtils.distanceVincenty(pLatitude, pLongitude, locationData[i][0], locationData[i][1]); + if (distace <= pMaxDistance) + locations.push(locationData[i]); + } + + return locations.sort(distanceSortFn); + + function distanceSortFn (a, b) + { + return a[3] - b[3]; + } +} \ No newline at end of file -- GitLab