/*
* niViz -- snow profiles visualization
* Copyright (C) 2015 WSL/SLF - Fluelastrasse 11 - 7260 Davos Dorf - Switzerland.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/*eslint yoda: 0, space-infix-ops: 0, complexity: 0*/
(function (niviz) {
'use strict';
// --- Module Dependencies ---
var properties = Object.defineProperties;
var Common = niviz.Common;
/** @module niviz */
/**
* Comprises the geo-coordinates of a station.
*
* @class Position
* @constructor
*/
function Position() {
/**
* @property uom
* @type String
*/
this.uom = 'm';
/**
* @property latitude
* @type Number
*/
/**
* @property longitude
* @type Number
*/
/**
* @property altitude
* @type Number
*/
/**
* @property x
* @type Number
*/
/**
* @property y
* @type Number
*/
/**
* Slope angle.
*
* @property angle
* @type Number
*/
/**
* Slope azimuth.
*
* @property azimuth
* @type Number
*/
}
/**
* @property directions
* @static
* @type {Array<String>}
*/
Position.directions = ['n/a', 'N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'];
Position.imagesrc = function (obj, zoom, size) {
zoom = zoom || 12;
size = size || '300x200';
if ((obj.latitude || obj.latitude === 0) &&
(obj.longitude || obj.longitude === 0)) {
var url = 'https://maps.googleapis.com/maps/api/staticmap?center='
+ obj.latitude + ',' + obj.longitude
+ '&zoom=' + zoom + '&size=' + size + '&sensor=false'
+ '&maptype=terrain'
+ '&key=AIzaSyCJs4S15Pd8Kyc4H8N63ro7rUKwWV6cEoM'
+ '&visual_refresh=true&markers=size:mid%7Ccolor:0xff0000%7C'
+ obj.latitude + ',' + obj.longitude;
return url;
}
return null;
};
properties(Position.prototype, {
/**
* Get the URL to a map service for this position.
* @property link
* @type {String}
*/
link: {
get: function () {
var service = Common.defaults.map_service;
if (service === 'maps.google.com') {
return 'http://maps.google.com/maps?z=18&t=m&q=loc:'
+ Number(this.latitude).toFixed(10) + '+' + Number(this.longitude).toFixed(10);
} else if (service === 'mapy.cz') {
return 'https://mapy.cz/zimni?'
+ 'x=' + Number(this.longitude).toFixed(10) + '&y='
+ Number(this.latitude).toFixed(10)
+ '&z=14&l=0&pano=1&source=base';
} else {
if (!this.x || !this.y)
return 'http://map.geo.admin.ch/?swisssearch='
+ Number(this.latitude).toFixed(10) + ',' + Number(this.longitude).toFixed(10)
+ '&zoom=8&crosshair=circle&bgLayer=ch.swisstopo.pixelkarte-farbe'
+ '&time_current=latest&lang=de&topic=ech&layers_opacity=0.35'
+ '&layers=ch.swisstopo.hangneigung-ueber_30';
return 'https://map.geo.admin.ch/?Y=' + this.x + '&X=' + this.y
+ '&zoom=8&crosshair=circle&bgLayer=ch.swisstopo.pixelkarte-farbe'
+ '&time_current=latest&lang=de&topic=ech&layers_opacity=0.35'
+ '&layers=ch.swisstopo.hangneigung-ueber_30';
}
}
},
/**
* Get imagesrc url for the current location.
* @property imagesrc
* @type {String}
*/
imagesrc: {
get: function () {
return Position.imagesrc(this);
}
},
/**
* Get the exposition of a location.
* @property direction
* @type {String}
*/
direction: {
get: function () {
if (this.aspect) {
if (this.aspect !== 'n/a') return this.aspect;
if (this.aspect === 'n/a' && this.angle === 0) return 'flat';
}
var direction = '';
if (this.azimuth !== undefined) {
var index = this.azimuth % 360;
if (index < 0 ) index = 360 + index;
index = Math.floor((index + 22.5) / 360 * 8) % 8;
direction = Position.directions[index];
}
return direction;
}
},
azimuth: {
get: function () {
return this._azimuth;
},
set: function (value) {
this._azimuth = value;
this._aspect = Position.angle_to_dir(value);
}
},
aspect: {
get: function () {
return this._aspect;
},
set: function (value) {
this._aspect = value;
this._azimuth = Position.dir_to_angle(value);
}
},
utm: {
get: function () {
var zone = null, x = null, y = null, conversion;
if ((this.latitude || this.latitude === 0) &&
(this.longitude || this.longitude === 0)) {
conversion = Position.getUTMZone(this.latitude, this.longitude);
zone = conversion.number + conversion.letter;
conversion = Position.WGS84_to_UTM(this.latitude, this.longitude);
x = conversion.east;
y = conversion.north;
}
return { zone: zone, east: x, north: y };
}
}
});
/**
* Convert direction to angle
*
* @method dir_to_angle
* @static
*
* @param {String} intercardinal or cardinal direction
*
* @return {Number} Angle [0; 360]
*/
Position.dir_to_angle = function (direction) {
if (direction === 'N') return 0;
if (direction === 'NE') return 45;
if (direction === 'E') return 90;
if (direction === 'SE') return 135;
if (direction === 'S') return 180;
if (direction === 'SW') return 225;
if (direction === 'W') return 270;
if (direction === 'NW') return 315;
return null;
};
/**
* Convert angle to direction
*
* @method angle_to_dir
* @static
*
* @param {Number} angle [0; 360]
*
* @return {String} intercardinal or cardinal direction
*/
Position.angle_to_dir = function (angle) {
if (!angle && angle !== 0) return null;
if (typeof angle !== 'number')
throw new Error('Error while converting angle to direction: angle not a number');
angle = angle % 360;
if (angle >= 337.5 || angle < 22.5) return 'N';
else if (angle >= 22.5 && angle < 67.5) return 'NE';
else if (angle >= 67.5 && angle < 112.5) return 'E';
else if (angle >= 112.5 && angle < 157.5) return 'SE';
else if (angle >= 157.5 && angle < 202.5) return 'S';
else if (angle >= 202.5 && angle < 247.5) return 'SW';
else if (angle >= 247.5 && angle < 292.5) return 'W';
else if (angle >= 292.5 && angle < 337.5) return 'NW';
return null;
};
var to_rad = Math.PI / 180.;
var to_deg = 180. / Math.PI;
function fmod (a, b) {
return Number((a - (Math.floor(a / b) * b)).toPrecision(8));
}
Position.ellipsoid = {
a: 6378137.,
b: 6356752.3142
};
/**
* returns the epsg code matching certain UTM zone
*
* @param {String} coordparam UTM zone
* @return epsg code
*/
Position.zone_to_epsg = function (zone) {
// UTM Zone information
var parsedZone = zone.match(/(\d{1,2})([A-Z])/);
var zoneNumber = parseInt(parsedZone[1]);
var zoneLetter = parsedZone[2];
if (zoneLetter >= 'N') {
return 32600 + zoneNumber;
} else { //southern hemisphere
return 32700 + zoneNumber;
}
};
/**
* returns the epsg code matching certain UTM zone
*
* @param {String} coordparam UTM zone
* @return epsg code
*/
Position.epsg_to_zone = function (epsg) {
var letter, number;
if ((epsg >= 32601) && (epsg <= 32660)) {
//northern hemisphere
number = epsg - 32600;
letter = 'N';
} else if ((epsg >= 32701) && (epsg <= 32760)) {
//southern hemisphere
number = epsg - 32700;
letter = 'M';
}
return {
letter: letter,
number: number
};
};
/**
* This routine determines the correct UTM letter designator for the given latitude␊
* UTM limits its coverage to [80S , 84N], outside of this, returns Y/Z/A/B for the zone
* Copied from MeteoIO: dataClasses/CoordsAlgorithms.cc
*
* @method getUTMZone
* @static
*
* @param {Number} i_latitude
* @param {Number} i_longitude
*
* @return {String} concat zone number + zone letter
*/
Position.getUTMZone = function (i_latitude, i_longitude) {
//computing zone number, assuming longitude in [-180. ; 180[
var ZoneNumber = Math.round((i_longitude + 180.) / 6.) + 1;
// Special zones for Scandinavia
if ( i_latitude >= 72.0 && i_latitude < 84.0 ) {
if ( i_longitude >= 0.0 && i_longitude < 9.0 ) ZoneNumber = 31;
else if ( i_longitude >= 9.0 && i_longitude < 21.0 ) ZoneNumber = 33;
else if ( i_longitude >= 21.0 && i_longitude < 33.0 ) ZoneNumber = 35;
else if ( i_longitude >= 33.0 && i_longitude < 42.0 ) ZoneNumber = 37;
}
if ( i_latitude >= 56.0 && i_latitude < 64.0 &&
i_longitude >= 3.0 && i_longitude < 12.0 ) {
ZoneNumber = 32;
}
//getting zone letter
var zoneLetter = 'Z';
if ((0 >= i_longitude) && (i_latitude > 84)) zoneLetter = 'Y';
else if ((0 < i_longitude) && (i_latitude > 84)) zoneLetter = 'Z';
else if ((84 >= i_latitude) && (i_latitude >= 72)) zoneLetter = 'X';
else if ((72 > i_latitude) && (i_latitude >= 64)) zoneLetter = 'W';
else if ((64 > i_latitude) && (i_latitude >= 56)) zoneLetter = 'V';
else if ((56 > i_latitude) && (i_latitude >= 48)) zoneLetter = 'U';
else if ((48 > i_latitude) && (i_latitude >= 40)) zoneLetter = 'T';
else if ((40 > i_latitude) && (i_latitude >= 32)) zoneLetter = 'S';
else if ((32 > i_latitude) && (i_latitude >= 24)) zoneLetter = 'R';
else if ((24 > i_latitude) && (i_latitude >= 16)) zoneLetter = 'Q';
else if ((16 > i_latitude) && (i_latitude >= 8)) zoneLetter = 'P';
else if (( 8 > i_latitude) && (i_latitude >= 0)) zoneLetter = 'N';
else if (( 0 > i_latitude) && (i_latitude >= -8)) zoneLetter = 'M';
else if ((-8 > i_latitude) && (i_latitude >= -16)) zoneLetter = 'L';
else if ((-16 > i_latitude) && (i_latitude >= -24)) zoneLetter = 'K';
else if ((-24 > i_latitude) && (i_latitude >= -32)) zoneLetter = 'J';
else if ((-32 > i_latitude) && (i_latitude >= -40)) zoneLetter = 'H';
else if ((-40 > i_latitude) && (i_latitude >= -48)) zoneLetter = 'G';
else if ((-48 > i_latitude) && (i_latitude >= -56)) zoneLetter = 'F';
else if ((-56 > i_latitude) && (i_latitude >= -64)) zoneLetter = 'E';
else if ((-64 > i_latitude) && (i_latitude >= -72)) zoneLetter = 'D';
else if ((-72 > i_latitude) && (i_latitude >= -80)) zoneLetter = 'C';
else if ((0 >= i_longitude) && (i_latitude <= -80)) zoneLetter = 'A';
else if ((0 < i_longitude) && (i_latitude <= -80)) zoneLetter = 'B';
return {
number: ZoneNumber,
letter: zoneLetter
};
};
/**
* @brief Coordinate conversion: from WGS84 Lat/Long to UTM grid.␊
* Copied from MeteoIO: dataClasses/CoordsAlgorithms.cc
*
* @param[in] lat_in Decimal Latitude
* @param[in] long_in Decimal Longitude
* @param[in] coordparam UTM zone to convert to
* @param[out] east_out easting coordinate (Swiss system)
* @param[out] north_out northing coordinate (Swiss system)
*/
Position.WGS84_to_UTM = function (lat_in, long_in) {
var east_out, north_out;
//Geometric constants
var a = Position.ellipsoid.a; //major ellipsoid semi-axis
var b = Position.ellipsoid.b;//minor ellipsoid semi-axis
var e2 = (a*a-b*b) / (a*a);//ellispoid eccentricity, squared (=(a²-b²)/a²)
var eP2 = e2 / (1.-e2);//second ellispoid eccentricity, squared (=(a²-b²)/b²)
var k0 = 0.9996;//scale factor for the projection
//getting position parameters
var zone;
long_in = fmod(long_in+360.+180., 360.) - 180.; //normalized to [-180. ; 180.[
var Long = long_in * to_rad;
var Lat = lat_in * to_rad;
var zoneNumber = Position.getUTMZone(lat_in, long_in, zone).number;
//+3 puts origin in middle of zone
var long0 = ((zoneNumber - 1)*6 - 180 + 3) * to_rad;
//Geometrical parameters
//radius of curvature of the earth perpendicular to the meridian plane
var nu = a / Math.sqrt(1.-e2*Math.pow(Math.sin(Lat), 2));
var p = (Long-long0);
//calculating first the coefficients of the series, then the Meridional Arc M itself
var n = (a-b)/(a+b);
var n2=n*n, n3=n*n*n, n4=n2*n2, n5=n2*n3, n6=n3*n3;
var A = a * (1. - n + 5./4.*(n2 - n3) + 81./64.*(n4 - n5));
var B = (3./2.*a) * (n - n2 + 7./8.*(n3 - n4) + 55./64.*(n5 - n6));
var C = (15./16.*a) * (n2 - n3 + 3./4.*(n4 - n5));
var D = (35./48.*a) * (n3 - n4 + 11./16.*(n5 - n6));
var E = (315./512.*a) * (n4 - n5); //correction of ~0.03mm
var M = A*Lat - B*Math.sin(2.*Lat) + C*Math.sin(4.*Lat)
- D*Math.sin(6.*Lat) + E*Math.sin(8.*Lat);
//calculating the coefficients for the series
var K1 = M*k0;
var K2 = 1./4.*k0*nu*Math.sin(2.*Lat);
var K3 = (k0*nu*Math.sin(Lat)*Math.pow(Math.cos(Lat), 3)*1./24.)
* (5. - Math.pow(Math.tan(Lat), 2) + 9.*eP2*Math.pow(Math.cos(Lat), 2)
+ 4.*eP2*eP2*Math.pow(Math.cos(Lat), 4));
var K4 = k0*nu*Math.cos(Lat);
var K5 = (k0*nu*Math.pow(Math.cos(Lat), 3)*1./6.)
* (1. - Math.pow(Math.tan(Lat), 2) + eP2*Math.pow(Math.cos(Lat), 2));
north_out = K1 + K2*p*p + K3*p*p*p*p;
east_out = K4*p + K5*p*p*p + 500000.0;
if (Lat < 0)
north_out += 10000000.0; //offset for southern hemisphere
return {
east: east_out,
north: north_out
};
};
/**
* Coordinate conversion: from UTM grid to WGS84 Lat/Long.
* Copied from MeteoIO: dataClasses/CoordsAlgorithms.cc
*
* @param[in] east_in easting coordinate (UTM)
* @param[in] north_in northing coordinate (UTM)
* @param[in] zone UTM zone of the easting/northing
* @param[out] lat_out Decimal Latitude
* @param[out] long_out Decimal Longitude
*/
Position.UTM_to_WGS84 = function (east_in, north_in, zone) {
var lat_out, long_out, sin = Math.sin;
// Geometric constants
var a = Position.ellipsoid.a; //major ellipsoid semi-axis
var b = Position.ellipsoid.b;//minor ellipsoid semi-axis
var e2 = (a*a-b*b) / (a*a);//ellispoid eccentricity, squared (=(a²-b²)/a²)
var eP2 = e2 / (1.-e2);//second ellispoid eccentricity, squared (=(a²-b²)/b²)
var k0 = 0.9996;//scale factor for the projection
// UTM Zone information
var parsedZone = zone.match(/(\d{1,2})([A-Z])/);
var zoneNumber = parseInt(parsedZone[1]);
var zoneLetter = parsedZone[2];
// set reference parameters: central meridian of the zone, true northing and easting
// please note that the special zones still use the reference meridian as given by
// their zone number (ie: even if it might not be central anymore)
// +3 puts origin in "middle" of zone as required for the projection
// meridian (might not be the middle for special zones)
var long0 = (zoneNumber - 1)*6 - 180 + 3;
if (zoneLetter <= 'M') {
north_in -= 10000000.0; //offset used for southern hemisphere
}
east_in -= 500000.0; //longitude offset: x coordinate is relative to central meridian
//calculating footprint latitude fp (it should be done using a few iterations)
var arc = north_in/k0; //Meridional arc
var mu = arc / (a*(1.-e2/4.-3.*e2*e2/64.-5.*e2*e2*e2/256.));
var e1 = (1.-b/a) / (1.+b/a); //simplification of [1 - (1 - e2)1/2]/[1 + (1 - e2)1/2]
var J1 = (3./2.*e1 - 27./32.*e1*e1*e1);
var J2 = (21./16.*e1*e1 - 55./32.*e1*e1*e1*e1);
var J3 = (151./96.*e1*e1*e1);
var J4 = (1097./512.*e1*e1*e1*e1);
var fp = mu + J1*sin(2.*mu) + J2*sin(4.*mu) + J3*sin(6.*mu) + J4*sin(8.*mu);
//calculating the parameters
var C1 = eP2 * Math.pow(Math.cos(fp), 2);
var T1 = Math.pow( Math.tan(fp), 2 );
var R1 = a*(1.-e2) / Math.pow((1.-e2*Math.pow(sin(fp), 2)), 1.5);
var N1 = a / Math.sqrt(1.-e2*Math.pow(sin(fp), 2));
var D = east_in / (N1*k0);
var D2=D*D, D3=D*D*D;
//calculating the coefficients of the series for latitude and longitude
var Q1 = N1*Math.tan(fp)/R1;
var Q2 = 0.5*D2;
var Q3 = (5. + 3.*T1 + 10.*C1 - 4.*C1*C1 - 9.*eP2) * 1./24.*D2*D2;
var Q4 = (61. + 90.*T1 + 298.*C1 + 45.*T1*T1 - 3.*C1*C1 - 252.*eP2) * 1./720.*D3*D3;
var Q5 = D;
var Q6 = (1. + 2.*T1 + C1) * 1./6.*D3;
var Q7 = (5. - 2.*C1 + 28.*T1 - 3.*C1*C1 + 8.*eP2 + 24.*T1*T1) * 1./120.*D2*D3;
lat_out = (fp - Q1 * (Q2 - Q3 + Q4))*to_deg;
long_out = long0 + ((Q5 - Q6 + Q7)/Math.cos(fp))*to_deg;
return {
latitude: lat_out,
longitude: long_out
};
};
Position.prototype.toString = function () {
};
/**
* Clone current object
*
* @method clone
* @return {Position}
*/
Position.prototype.clone = function () {
var pos = new Position();
var fields = [ 'latitude', 'longitude', 'altitude', 'x', 'y',
'angle', 'azimuth', 'aspect', 'description' ];
fields.forEach(function (prop) {
if (this[prop] !== undefined) pos[prop] = this[prop];
}, this);
return pos;
};
// --- Module Export ---
niviz.Position = Position;
}(niviz));
