79 lines
3.4 KiB
JavaScript
79 lines
3.4 KiB
JavaScript
// Converts UTM coordinates to latitude and longitude
|
|
|
|
const a = 6378137.0; // WGS84 Major axis
|
|
const eccSquared = 0.00669438; // WGS84 eccentricity squared
|
|
const k0 = 0.9996; // Scale factor
|
|
|
|
function UTMtoLatLon(UTMEasting, UTMNorthing, UTMZoneNumber, Hemisphere) {
|
|
const e1 = (1 - Math.sqrt(1 - eccSquared)) / (1 + Math.sqrt(1 - eccSquared));
|
|
let x = UTMEasting - 500000.0;
|
|
let y = UTMNorthing;
|
|
|
|
if (Hemisphere === 'S') {
|
|
y -= 10000000.0;
|
|
}
|
|
|
|
const LongOrigin = (UTMZoneNumber - 1) * 6 - 180 + 3;
|
|
const eccPrimeSquared = eccSquared / (1 - eccSquared);
|
|
const M = y / k0;
|
|
const mu = M / (a * (1 - eccSquared / 4 - (3 * eccSquared ** 2) / 64 - (5 * eccSquared ** 3) / 256));
|
|
|
|
let phi1Rad = mu + ((3 * e1) / 2 - (27 * e1 ** 3) / 32) * Math.sin(2 * mu) + ((21 * e1 ** 2) / 16 - (55 * e1 ** 4) / 32) * Math.sin(4 * mu) + ((151 * e1 ** 3) / 96) * Math.sin(6 * mu);
|
|
|
|
const N1 = a / Math.sqrt(1 - eccSquared * Math.sin(phi1Rad) ** 2);
|
|
const T1 = Math.tan(phi1Rad) ** 2;
|
|
const C1 = eccPrimeSquared * Math.cos(phi1Rad) ** 2;
|
|
const R1 = (a * (1 - eccSquared)) / Math.pow(1 - eccSquared * Math.sin(phi1Rad) ** 2, 1.5);
|
|
const D = x / (N1 * k0);
|
|
|
|
let lat = phi1Rad - ((N1 * Math.tan(phi1Rad)) / R1) * (D ** 2 / 2 - ((5 + 3 * T1 + 10 * C1 - 4 * C1 ** 2 - 9 * eccPrimeSquared) * D ** 4) / 24 + ((61 + 90 * T1 + 298 * C1 + 45 * T1 ** 2 - 252 * eccPrimeSquared - 3 * C1 ** 2) * D ** 6) / 720);
|
|
lat *= 180 / Math.PI;
|
|
|
|
let lon = (D - ((1 + 2 * T1 + C1) * D ** 3) / 6 + ((5 - 2 * C1 + 28 * T1 - 3 * C1 ** 2 + 8 * eccPrimeSquared + 24 * T1 ** 2) * D ** 5) / 120) / Math.cos(phi1Rad);
|
|
lon = LongOrigin + (lon * 180) / Math.PI;
|
|
|
|
return { latitude: lat, longitude: lon };
|
|
}
|
|
|
|
// Converts latitude and longitude to UTM coordinates
|
|
function LatLontoUTM(lat, lon) {
|
|
const LatRad = (lat * Math.PI) / 180.0;
|
|
const LongRad = (lon * Math.PI) / 180.0;
|
|
|
|
const ZoneNumber = Math.floor((lon + 180) / 6) + 1;
|
|
const LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3;
|
|
const LongOriginRad = (LongOrigin * Math.PI) / 180.0;
|
|
|
|
const eccPrimeSquared = eccSquared / (1 - eccSquared);
|
|
const N = a / Math.sqrt(1 - eccSquared * Math.sin(LatRad) ** 2);
|
|
const T = Math.tan(LatRad) ** 2;
|
|
const C = eccPrimeSquared * Math.cos(LatRad) ** 2;
|
|
const A = Math.cos(LatRad) * (LongRad - LongOriginRad);
|
|
|
|
const M =
|
|
a *
|
|
((1 - eccSquared / 4 - (3 * eccSquared ** 2) / 64 - (5 * eccSquared ** 3) / 256) * LatRad -
|
|
((3 * eccSquared) / 8 + (3 * eccSquared ** 2) / 32 + (45 * eccSquared ** 3) / 1024) * Math.sin(2 * LatRad) +
|
|
((15 * eccSquared ** 2) / 256 + (45 * eccSquared ** 3) / 1024) * Math.sin(4 * LatRad) -
|
|
((35 * eccSquared ** 3) / 3072) * Math.sin(6 * LatRad));
|
|
|
|
let UTMEasting = k0 * N * (A + ((1 - T + C) * A ** 3) / 6 + ((5 - 18 * T + T ** 2 + 72 * C - 58 * eccPrimeSquared) * A ** 5) / 120) + 500000.0;
|
|
|
|
let UTMNorthing = k0 * (M + N * Math.tan(LatRad) * (A ** 2 / 2 + ((5 - T + 9 * C + 4 * C ** 2) * A ** 4) / 24 + ((61 - 58 * T + T ** 2 + 600 * C - 330 * eccPrimeSquared) * A ** 6) / 720));
|
|
|
|
const Hemisphere = lat < 0 ? 'S' : 'N';
|
|
if (lat < 0) {
|
|
UTMNorthing += 10000000.0;
|
|
}
|
|
|
|
return {
|
|
easting: UTMEasting,
|
|
northing: UTMNorthing,
|
|
zoneNumber: ZoneNumber,
|
|
hemisphere: Hemisphere
|
|
};
|
|
}
|
|
|
|
// Exporting functions as a module
|
|
export { UTMtoLatLon, LatLontoUTM };
|