javascriptgiscoordinate-systemscoordinate-transformation

Converting lat/long coordinates to UTM using Javascript - wrong easting

I am trying to convert these coordinates to UTM

const longitude = 143.5279742;
const latitude = 44.062027;


// Calculate the UTM coordinates
const zoneNumber = Math.floor((longitude + 180) / 6) + 1;

const k = 0.9996;
const e = 0.081819190842965;
const e2 = e * e;
const e4 = e2 * e2;
const e6 = e4 * e2;
const N = 6378137 / Math.sqrt(1 - e2 * Math.pow(Math.sin(latitude * Math.PI / 180), 2));
const T = Math.pow(Math.tan(latitude * Math.PI / 180), 2);
const C = e4 * Math.pow(Math.cos(latitude * Math.PI / 180), 2);
const A = (longitude - ((zoneNumber - 1) * 6 - 180 + 3)) * Math.PI / 180;
const M = 6378137 * ((1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256) * latitude * Math.PI / 180 
                   - (3 * e2 / 8 + 3 * e4 / 32 + 45 * e6 / 1024) * Math.sin(2 * latitude * Math.PI / 180)
                   + (15 * e4 / 256 + 45 * e6 / 1024) * Math.sin(4 * latitude * Math.PI / 180)
                   - (35 * e6 / 3072) * Math.sin(6 * latitude * Math.PI / 180));

const A_prime = A - (1 - T + C) * Math.pow(A, 3) / 6 + (5 - 18 * T + Math.pow(T, 2) + 72 * C - 58 * e2) * Math.pow(A, 5) / 120;

const x = k * N * (A_prime * 180 / Math.PI + 0.9996 * 500000 + 0.5);
const y = k * (M + N * Math.tan(latitude * Math.PI / 180) * (Math.pow(A_prime, 2) / 2 + (5 - T + 9 * C + 4 * C * C) * Math.pow(A_prime, 4) / 24 + (61 - 58 * T + Math.pow(T, 2) + 600 * C - 330 * e2) * Math.pow(A_prime, 6) / 720));

// Format the output as a string
const utmString = zoneNumber + " " + Math.round(x) + " " + Math.round(y);



text_UTM.setValue(utmString);

Blockquote

The northing produced and zone are correct however the easting is wrong. it should produce the easting 702469.05 but the number is completely wrong. Any advice?

Solution

  • I recommend to use this library https://www.npmjs.com/package/utm-latlng it has already the formulas for the convertion

    utm.convertLatLngToUtm(44.062027, 143.5279742, 1);

    Take a look to this fiddle, I have copied the entire library, sorry I haven't found a CDN of this, this is the file used: https://github.com/shahid28/utm-latlng/blob/master/UTMLatLngFront.js

    Result:

    {"Easting":702469,"Northing":4881869,"ZoneNumber":54,"ZoneLetter":"T"}


    /*********************/
/** start copied from https://github.com/shahid28/utm-latlng/blob/master/UTMLatLngFront.js **/
/*********************/
/*********************/
/*********************/

var datumName = "WGS 84";
var a;
var eccSquared;
var status = false;

function UTMLatLng (datumNameIn){

  if (datumNameIn !== undefined)
  {
      datumName = datumNameIn;
  }
  this.setEllipsoid(datumName);
}

UTMLatLng.prototype.convertLatLngToUtm = function (latitude, longitude,precision){
    if (this.status)
    {
        return 'No ecclipsoid data associated with unknown datum: ' + datumName;
    }


    if(!Number.isInteger(precision))
    {
        return 'Precision is not integer number.';
    }

    latitude = parseFloat(latitude);
    longitude = parseFloat(longitude);

    var LongTemp = longitude;
    var LatRad = this.toRadians(latitude);
    var LongRad = this.toRadians(LongTemp);

    if (LongTemp >= 8 && LongTemp <= 13 && latitude > 54.5 && latitude < 58) {
        ZoneNumber = 32;
    } else if (latitude >= 56.0 && latitude < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0) {
        ZoneNumber = 32;
    } else {
        ZoneNumber = ((LongTemp + 180) / 6) + 1;

        if (latitude >= 72.0 && latitude < 84.0) {
            if (LongTemp >= 0.0 && LongTemp < 9.0) {
                ZoneNumber = 31;
            } else if (LongTemp >= 9.0 && LongTemp < 21.0) {
                ZoneNumber = 33;
            } else if (LongTemp >= 21.0 && LongTemp < 33.0) {
                ZoneNumber = 35;
            } else if (LongTemp >= 33.0 && LongTemp < 42.0) {
                ZoneNumber = 37;
            }
        }
    }
    ZoneNumber = parseInt(ZoneNumber);

    var LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3;  //+3 puts origin in middle of zone
    var LongOriginRad = this.toRadians(LongOrigin);

    var UTMZone = this.getUtmLetterDesignator(latitude);

    var eccPrimeSquared = (this.eccSquared) / (1 - this.eccSquared);

    var N = this.a / Math.sqrt(1 - this.eccSquared * Math.sin(LatRad) * Math.sin(LatRad));
    var T = Math.tan(LatRad) * Math.tan(LatRad);
    var C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad);
    var A = Math.cos(LatRad) * (LongRad - LongOriginRad);

    var M = this.a * ((1 - this.eccSquared / 4 - 3 * this.eccSquared * this.eccSquared / 64 - 5 * this.eccSquared * this.eccSquared * this.eccSquared / 256) * LatRad
            - (3 * this.eccSquared / 8 + 3 * this.eccSquared * this.eccSquared / 32 + 45 * this.eccSquared * this.eccSquared * this.eccSquared / 1024) * Math.sin(2 * LatRad)
            + (15 * this.eccSquared * this.eccSquared / 256 + 45 * this.eccSquared * this.eccSquared * this.eccSquared / 1024) * Math.sin(4 * LatRad)
            - (35 * this.eccSquared * this.eccSquared * this.eccSquared / 3072) * Math.sin(6 * LatRad));

    var UTMEasting = parseFloat(0.9996 * N * (A + (1 - T + C) * A * A * A / 6
            + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
            + 500000.0);

    var UTMNorthing = parseFloat(0.9996 * (M + N * Math.tan(LatRad) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
            + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)));

    if (latitude < 0)
        UTMNorthing += 10000000.0;
    UTMNorthing = parseInt(Math.round(UTMNorthing));
    UTMEasting = parseInt(Math.round(UTMEasting));
    return {Easting: UTMEasting, Northing: UTMNorthing, ZoneNumber: parseInt(ZoneNumber), ZoneLetter: UTMZone};


};

UTMLatLng.prototype.convertUtmToLatLng = function (UTMEasting, UTMNorthing, UTMZoneNumber, UTMZoneLetter){
    var e1 = (1 - Math.sqrt(1 - this.eccSquared)) / (1 + Math.sqrt(1 - this.eccSquared));
    var x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
    var y = UTMNorthing;
    var ZoneNumber = UTMZoneNumber;
    var ZoneLetter = UTMZoneLetter;
    var NorthernHemisphere;
    if(UTMEasting===undefined)
    {
        return "Please pass the UTMEasting!";
    }
    if(UTMNorthing===undefined)
    {
        return "Please pass the UTMNorthing!";
    }
    if(UTMZoneNumber===undefined)
    {
        return "Please pass the UTMZoneNumber!";
    }
    if(UTMZoneLetter===undefined)
    {
        return "Please pass the UTMZoneLetter!";
    }

    if ('N' === ZoneLetter) {
        NorthernHemisphere = 1;
    } else {
        NorthernHemisphere = 0;
        y -= 10000000.0;
    }

    var LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3;

    var eccPrimeSquared = (this.eccSquared) / (1 - this.eccSquared);

    M = y / 0.9996;
    var mu = M / (this.a * (1 - this.eccSquared / 4 - 3 * this.eccSquared * this.eccSquared / 64 - 5 * this.eccSquared * this.eccSquared * this.eccSquared / 256));

    var phi1Rad = mu + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.sin(2 * mu)
            + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.sin(4 * mu)
            + (151 * e1 * e1 * e1 / 96) * Math.sin(6 * mu);
    var phi1 = this.toDegrees(phi1Rad);

    var N1 = this.a / Math.sqrt(1 - this.eccSquared * Math.sin(phi1Rad) * Math.sin(phi1Rad));
    var T1 = Math.tan(phi1Rad) * Math.tan(phi1Rad);
    var C1 = eccPrimeSquared * Math.cos(phi1Rad) * Math.cos(phi1Rad);
    var R1 = this.a * (1 - this.eccSquared) / Math.pow(1 - this.eccSquared * Math.sin(phi1Rad) * Math.sin(phi1Rad), 1.5);
    var D = x / (N1 * 0.9996);

    var Lat = phi1Rad - (N1 * Math.tan(phi1Rad) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24
            + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720);
    Lat = this.toDegrees(Lat);

    var Long = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1)
            * D * D * D * D * D / 120) / Math.cos(phi1Rad);
    Long = LongOrigin + this.toDegrees(Long);
    return {lat: Lat, lng: Long};
};

UTMLatLng.prototype.getUtmLetterDesignator = function (latitude) {
    latitude = parseFloat(latitude);
    if ((84 >= latitude) && (latitude >= 72))
        return 'X';
    else if ((72 > latitude) && (latitude >= 64))
        return 'W';
    else if ((64 > latitude) && (latitude >= 56))
        return 'V';
    else if ((56 > latitude) && (latitude >= 48))
        return 'U';
    else if ((48 > latitude) && (latitude >= 40))
        return 'T';
    else if ((40 > latitude) && (latitude >= 32))
        return 'S';
    else if ((32 > latitude) && (latitude >= 24))
        return 'R';
    else if ((24 > latitude) && (latitude >= 16))
        return 'Q';
    else if ((16 > latitude) && (latitude >= 8))
        return 'P';
    else if ((8 > latitude) && (latitude >= 0))
        return 'N';
    else if ((0 > latitude) && (latitude >= -8))
        return 'M';
    else if ((-8 > latitude) && (latitude >= -16))
        return 'L';
    else if ((-16 > latitude) && (latitude >= -24))
        return 'K';
    else if ((-24 > latitude) && (latitude >= -32))
        return 'J';
    else if ((-32 > latitude) && (latitude >= -40))
        return 'H';
    else if ((-40 > latitude) && (latitude >= -48))
        return 'G';
    else if ((-48 > latitude) && (latitude >= -56))
        return 'F';
    else if ((-56 > latitude) && (latitude >= -64))
        return 'E';
    else if ((-64 > latitude) && (latitude >= -72))
        return 'D';
    else if ((-72 > latitude) && (latitude >= -80))
        return 'C';
    else
        return 'Z';

};

UTMLatLng.prototype.setEllipsoid = function (name){
    switch (name) {
        case 'Airy':
            this.a = 6377563;
            this.eccSquared = 0.00667054;
            break;
        case 'Australian National':
            this.a = 6378160;
            this.eccSquared = 0.006694542;
            break;
        case 'Bessel 1841':
            this.a = 6377397;
            this.eccSquared = 0.006674372;
            break;
        case 'Bessel 1841 Nambia':
            this.a = 6377484;
            this.eccSquared = 0.006674372;
            break;
        case 'Clarke 1866':
            this.a = 6378206;
            this.eccSquared = 0.006768658;
            break;
        case 'Clarke 1880':
            this.a = 6378249;
            this.eccSquared = 0.006803511;
            break;
        case 'Everest':
            this.a = 6377276;
            this.eccSquared = 0.006637847;
            break;
        case 'Fischer 1960 Mercury':
            this.a = 6378166;
            this.eccSquared = 0.006693422;
            break;
        case 'Fischer 1968':
            this.a = 6378150;
            this.eccSquared = 0.006693422;
            break;
        case 'GRS 1967':
            this.a = 6378160;
            this.eccSquared = 0.006694605;
            break;
        case 'GRS 1980':
            this.a = 6378137;
            this.eccSquared = 0.00669438;
            break;
        case 'Helmert 1906':
            this.a = 6378200;
            this.eccSquared = 0.006693422;
            break;
        case 'Hough':
            this.a = 6378270;
            this.eccSquared = 0.00672267;
            break;
        case 'International':
            this.a = 6378388;
            this.eccSquared = 0.00672267;
            break;
        case 'Krassovsky':
            this.a = 6378245;
            this.eccSquared = 0.006693422;
            break;
        case 'Modified Airy':
            this.a = 6377340;
            this.eccSquared = 0.00667054;
            break;
        case 'Modified Everest':
            this.a = 6377304;
            this.eccSquared = 0.006637847;
            break;
        case 'Modified Fischer 1960':
            this.a = 6378155;
            this.eccSquared = 0.006693422;
            break;
        case 'South American 1969':
            this.a = 6378160;
            this.eccSquared = 0.006694542;
            break;
        case 'WGS 60':
            this.a = 6378165;
            this.eccSquared = 0.006693422;
            break;
        case 'WGS 66':
            this.a = 6378145;
            this.eccSquared = 0.006694542;
            break;
        case 'WGS 72':
            this.a = 6378135;
            this.eccSquared = 0.006694318;
            break;
        case 'ED50':
            this.a = 6378388;
            this.eccSquared = 0.00672267;
            break; // International Ellipsoid
        case 'WGS 84':
        case 'EUREF89': // Max deviation from WGS 84 is 40 cm/km see http://ocq.dk/euref89 (in danish)
        case 'ETRS89': // Same as EUREF89
            this.a = 6378137;
            this.eccSquared = 0.00669438;
            break;
        default:
            this.status = true;
            //   new Error('No ecclipsoid data associated with unknown datum: '.name);

    }
};

UTMLatLng.prototype.toDegrees = function (rad) {

    return rad / Math.PI * 180;
};

UTMLatLng.prototype.toRadians = function (deg) {

    return deg * Math.PI / 180;
};

function precisionRound(number, precision) {
      var factor = Math.pow(10, precision);
      return Math.round(number * factor) / factor;
  }
  
/*********************/
/** end copied from https://github.com/shahid28/utm-latlng/blob/master/UTMLatLngFront.js **/
/*********************/
/*********************/
/*********************/
  

const utm = new UTMLatLng();
const result = utm.convertLatLngToUtm(44.062027, 143.5279742, 1);
console.log(result);
document.querySelector(".result").innerHTML = JSON.stringify(result);
    <p class="result"></p>