Hi Al,
Yes. Although I am only converting the function that converts latitude and longitude to UTM coordinates:
/**
* Converts latitude/longitude to UTM coordinate.
*
* Implements Karney’s method, using Krüger series to order n^6, giving results accurate to 5nm for
* distances up to 3900km from the central meridian.
*
* @returns {Utm} UTM coordinate.
* @throws {Error} If point not valid, if point outside latitude range.
*
* @example
* var latlong = new LatLon(48.8582, 2.2945);
* var utmCoord = latlong.toUtm(); // utmCoord.toString(): '31 N 448252 5411933'
*/
LatLon.prototype.toUtm = function() {
if (isNaN(this.lat) || isNaN(this.lon)) throw new Error('Invalid point');
if (!(-80<=this.lat && this.lat<=84)) throw new Error('Outside UTM limits');
var falseEasting = 500e3, falseNorthing = 10000e3;
var zone = Math.floor((this.lon+180)/6) + 1; // longitudinal zone
var λ0 = ((zone-1)*6 - 180 + 3).toRadians(); // longitude of central meridian
// ---- handle Norway/Svalbard exceptions
// grid zones are 8° tall; 0°N is offset 10 into latitude bands array
var mgrsLatBands = 'CDEFGHJKLMNPQRSTUVWXX'; // X is repeated for 80-84°N
var latBand = mgrsLatBands.charAt(Math.floor(this.lat/8+10));
// adjust zone & central meridian for Norway
if (zone==31 && latBand=='V' && this.lon>= 3) { zone++; λ0 += (6).toRadians(); }
// adjust zone & central meridian for Svalbard
if (zone==32 && latBand=='X' && this.lon< 9) { zone--; λ0 -= (6).toRadians(); }
if (zone==32 && latBand=='X' && this.lon>= 9) { zone++; λ0 += (6).toRadians(); }
if (zone==34 && latBand=='X' && this.lon< 21) { zone--; λ0 -= (6).toRadians(); }
if (zone==34 && latBand=='X' && this.lon>=21) { zone++; λ0 += (6).toRadians(); }
if (zone==36 && latBand=='X' && this.lon< 33) { zone--; λ0 -= (6).toRadians(); }
if (zone==36 && latBand=='X' && this.lon>=33) { zone++; λ0 += (6).toRadians(); }
var φ = this.lat.toRadians(); // latitude ± from equator
var λ = this.lon.toRadians() - λ0; // longitude ± from central meridian
var a = this.datum.ellipsoid.a, f = this.datum.ellipsoid.f;
// WGS 84: a = 6378137, b = 6356752.314245, f = 1/298.257223563;
var k0 = 0.9996; // UTM scale on the central meridian
// ---- easting, northing: Karney 2011 Eq 7-14, 29, 35:
var e = Math.sqrt(f*(2-f)); // eccentricity
var n = f / (2 - f); // 3rd flattening
var n2 = n*n, n3 = n*n2, n4 = n*n3, n5 = n*n4, n6 = n*n5; // TODO: compare Horner-form accuracy?
var cosλ = Math.cos(λ), sinλ = Math.sin(λ), tanλ = Math.tan(λ);
var τ = Math.tan(φ); // τ ≡ tanφ, τʹ ≡ tanφʹ; prime (ʹ) indicates angles on the conformal sphere
var σ = Math.sinh(e*Math.atanh(e*τ/Math.sqrt(1+τ*τ)));
var τʹ = τ*Math.sqrt(1+σ*σ) - σ*Math.sqrt(1+τ*τ);
var ξʹ = Math.atan2(τʹ, cosλ);
var ηʹ = Math.asinh(sinλ / Math.sqrt(τʹ*τʹ + cosλ*cosλ));
var A = a/(1+n) * (1 + 1/4*n2 + 1/64*n4 + 1/256*n6); // 2πA is the circumference of a meridian
var α = [ null, // note α is one-based array (6th order Krüger expressions)
1/2*n - 2/3*n2 + 5/16*n3 + 41/180*n4 - 127/288*n5 + 7891/37800*n6,
13/48*n2 - 3/5*n3 + 557/1440*n4 + 281/630*n5 - 1983433/1935360*n6,
61/240*n3 - 103/140*n4 + 15061/26880*n5 + 167603/181440*n6,
49561/161280*n4 - 179/168*n5 + 6601661/7257600*n6,
34729/80640*n5 - 3418889/1995840*n6,
212378941/319334400*n6 ];
var ξ = ξʹ;
for (var j=1; j<=6; j++) ξ += α[j] * Math.sin(2*j*ξʹ) * Math.cosh(2*j*ηʹ);
var η = ηʹ;
for (var j=1; j<=6; j++) η += α[j] * Math.cos(2*j*ξʹ) * Math.sinh(2*j*ηʹ);
var x = k0 * A * η;
var y = k0 * A * ξ;
// ---- convergence: Karney 2011 Eq 23, 24
var pʹ = 1;
for (var j=1; j<=6; j++) pʹ += 2*j*α[j] * Math.cos(2*j*ξʹ) * Math.cosh(2*j*ηʹ);
var qʹ = 0;
for (var j=1; j<=6; j++) qʹ += 2*j*α[j] * Math.sin(2*j*ξʹ) * Math.sinh(2*j*ηʹ);
var γʹ = Math.atan(τʹ / Math.sqrt(1+τʹ*τʹ)*tanλ);
var γʺ = Math.atan2(qʹ, pʹ);
var γ = γʹ + γʺ;
// ---- scale: Karney 2011 Eq 25
var sinφ = Math.sin(φ);
var kʹ = Math.sqrt(1 - e*e*sinφ*sinφ) * Math.sqrt(1 + τ*τ) / Math.sqrt(τʹ*τʹ + cosλ*cosλ);
var kʺ = A / a * Math.sqrt(pʹ*pʹ + qʹ*qʹ);
var k = k0 * kʹ * kʺ;
// ------------
// shift x/y to false origins
x = x + falseEasting; // make x relative to false easting
if (y < 0) y = y + falseNorthing; // make y in southern hemisphere relative to false northing
// round to reasonable precision
x = Number(x.toFixed(6)); // nm precision
y = Number(y.toFixed(6)); // nm precision
var convergence = Number(γ.toDegrees().toFixed(9));
var scale = Number(k.toFixed(12));
var h = this.lat>=0 ? 'N' : 'S'; // hemisphere
return new Utm(zone, h, x, y, this.datum, convergence, scale);
};
This is my unfinished code. Its just a case of working through the routine above and wondering how code can be so difficult to read. Its all untested and I'm not certain my hyperbolic functions are correct.
Code: Select all
Function ConvertToUTM pLat, pLng
# pLat pLon are numbers. negative values indicate
-- 0-180W in Latitude or
-- 0-90S in Longitude
if pLat is not a number then return "bad value - latitude"
if pLng is not a number then return "bad value - longitude"
if pLat <= -80 OR pLat >= 84 then return "Outside UTM Limits - Latitude"
if pLng < 180 OR pLat > 180 then return "Outside UTM Limits - Longitude"
put 500000 into tFalseEasting
put 10000000 into tFalseNorthing
put the floor of ((pLng+180)/6) + 1 into tZone-- longitude zone (sixty 6 degree zones)
put ((tZone-1)*6)-180 + 3 into tLngCM // Longitude of centeral meridan of zone
-- JavaScript code uses lambda0 for this variable
// ------ handle special cases Norway/Svalbard exceptions
// grid zones are 8 degrees tall; 0 deg North is offset 10 into latitude bands array
put "CDEFGHJKLMNPQRSTUVWXX" into tMgrsLatBands -- X is repeated for posn's 80-84 degrees North
put char (Floor((pLat/8) + 10 )) of tMgrsLatBands into tLatBand -- single letter - N is equator to 8 deg north
-- adjust tZone and central meridian for Norway special case
if (tZone=31 AND tLatBand="V" AND pLon>=3) then
add 1 to tzone
add 6 to tLngCM -- still in degrees
end if
if (tZone=32 AND tLatBand="X" AND pLon<9) then
subtract 1 from tzone
subtract 6 from tLngCM -- still in degrees
end if
if (tZone=32 AND tLatBand="X" AND pLon>=9) then
add 1 to tzone
add 6 to tLngCM -- still in degrees
end if
if (tZone=34 AND tLatBand="X" AND pLon<21) then
subtract 1 from tzone
subtract 6 from tLngCM -- still in degrees
end if
if (tZone=34 AND tLatBand="X" AND pLon>=21) then
add 1 to tzone
add 6 to tLngCM -- still in degrees
end if
if (tZone=36 AND tLatBand="X" AND pLon<33) then
subtract 1 from tzone
subtract 6 from tLngCM -- still in degrees
end if
if (tZone=36 AND tLatBand="X" AND pLon>=33) then
add 1 to tzone
add 6 to tLngCM -- still in degrees
end if
// *****************************************
put pLat*pi/180 into tPhi -- Latitude in radians
put (pLng-tLngCM)*pi/180 into tLambda -- longitude ref the longitude of central meriadian in radians
-- WGS84 datum conversion factors
put 6378137 into a
put 1/298.257223563 into f
put 0.9996 into kZero -- UTM scale on the central meridian
put sqrt(f*(2-f)) into e -- eccentricity
put f / (2-f) into n -- 3rd flattening
put n*n into n2
put n*n2 into n3
put n*n3 into n4
put n*n4 into n5
put n*n5 into n6
put cos(tLambda) into cosLambda
put sin(tLambda) into sinLambda
put tan(tLambda) into tanLambda
put tan(tPhi) into tTau
put e*tTau/(sqrt(1+tTau*tTau))
put sinh(e*atanh(e*tTau/sqrt(1+(tTau*tTau)))) into tDelta
end ConvertToUTM
Private Function aTanh pX
--returns the inverse hyperbolic tangent
put (ln((1+pX)/(1-pX)))/2 into tAnswer
return tAnswer
end aTanh
Private Function Sinh pX
put (exp(pX) - exp(-pX))/2 into tAnswer
return tAnswer
end Sinh
Simon