How Compute & Use UTM Scale Factor?

June 8, 2011

Recently i do some work on computing point scale factor in utm coordinate system. The input of system is Easting,Northing,height and zone number in UTM (WGS84) and output is point scale factor.

First of all i should mention that the thing that we use as scale factor most of the time is Grid factor  (Combined Scale Factor that is the product of UTM Scale factor and Elevation Factor). I explain more.

Consider three coordinate systems.

System1:The system that we do observation in the surface of the earth with traditional surveying instrument (like total station and theodolites but not GPS)

System2: WGS84 (for example) or any other ellipsoid that made a coordinate system (we make a coordinate system using an ellipsoid by projecting the point into the ellipsoid plane and measure geodetic fi and landa angles and alos height of the point)

System3: is a projection surface(UTM) . this is what made (fi,landa) to (x,y) and make it possible to use paper maps to represent the earth curved surface (Earth is a curved surface in all points)

for converting a length from system1 to susyem3 we should use two kind of scale factors.

for example you measure a length on earth with GPS (E,N,h) and want to convert it to local length . what should we do?

first you should compute a scale factor that bring the point of observation from elevation h  into the surface of ellipsoid (h=0). it is because all of computations is on the surface. the name of this kind of scale factor is "Elevation Factor" (EF). it's formula is a very simple relation.

EF=(Re/Re+h) where Re= mean earth radius =6371000m and "h" is ellipsoid height of point.

second scale factor is what that transform between projection coordinate system(UTM) and WGS84 coordinate system.it's formula (with negligible truncation) is:

k=k0(1+(landa-landa0)^2*cos(fi)^2/2) where:

k=scale factor

k0=central meridian scale factor =0.9996 (in UTM Coordinate system)

(fi,landa) are the position of the point in WGS84

landa0=longitude of central meridian in that zone in utm Coordinate System

and at last for computing "grid scale factor" or combined scale factor,  you should multiply k*EF=Grid Factor

Grid factor is used to convert length from surface of earth into projection system (UTM) and if you want convert UTM length into earth surface length you should divide length into the Grid Factor.

you are done here but if you want go further, in measuring a length in fact we have 2 points, Start point and End point. it is better that you compute point scale factor in Start and End (or Start Middle and End) and then use average (or Weighted average) Scale Factor. in short distances it is not problem to use ONE SF but in long distances you need to average.

I wrote a C# program that takes E.N,h and zone and compute (fi,landa) and then Scal factor,Elevation Factor and finally Grid factor (combined Factor). here is a source code.

For rebuilding project you should just make a form and put some textbox controls and give desired input:

 

using System;

using System.Collections.Generic;

using System.ComponentModel;

using System.Data;

using System.Drawing;

using System.Linq;

using System.Text;

using System.Windows.Forms;

namespace UTM

{

    public partial class Form1 : Form

    {

        public double eccSquared = 0.00669438; // eccentricity (0.081819191 ^ 2) WGS84

        public double dEquatorialRadius = 6378137.0; // WGS84 (note above: varies from 6,356.750 km to 6,378.135 km)

        public double dMeanRadius = 6371000;//Earth Mean Radius

        public double dScaleFactor = 0.9996; // scale factor, used as k0

        public double dDenominatorOfFlatteningRatio = 298.257223563;

        public double

                      dCvtLiterHa2GalAc = 0.10691,   // l/ha -> gal/ac (dCvtLiter2Gal / dCvtHa2Ac)

                      dCvtLiterHa2OzAc = 13.68416,   // l/ha -> fl oz/ac (dCvtLiter2Oz / dCvtHa2Ac)

                      dCvtDeg2Rad = Math.PI / 180.0, // 0.0174532925199432957 ... remember: pi radians = 180 deg

                      dCvtRad2Deg = 180.0 / Math.PI; // 57.2957795130823208767 ...

 

        public Form1()

        {

            InitializeComponent();

        }

 

        public int iUTM2LatLon(double UTMNorthing,

              double UTMEasting,

              string sUTMZone,  // expected format "12N"

              ref double dLat,

              ref double dLon)

        {

            // if deviation from WGS84 is desired, do this (after loading the array, duh):

            //  cEllipsoid[] ellipsoidRg = EllipsoidLoad();

            //  dEquatorialRadius = ellipsoid[index-of-desired-reference-ellipsoid].EquatorialRadius;

            //  eccSquared = ellipsoid[index-of-desired-reference-ellipsoid].eccentricitySquared;

 

            // populate North/South

            char cZoneLetter = sUTMZone[sUTMZone.Length - 1];

            bool bNorthernHemisphere = (cZoneLetter >= 'N');

 

            string sZoneNum = sUTMZone.Substring(0, sUTMZone.Length - 1);

            int iZoneNumber = Convert.ToInt32(sZoneNum);

 

            double x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude

            double y = UTMNorthing;

            if (!bNorthernHemisphere) // point is in southern hemisphere

                y -= 10000000.0; // remove 10,000,000 meter offset used for southern hemisphere

 

            double dLongOrigin = (iZoneNumber - 1) * 6 - 180 + 3; // +3 puts origin in middle of zone

 

            double eccPrimeSquared = (eccSquared) / (1 - eccSquared);

 

            double M = y / dScaleFactor;

            double mu = M / (dEquatorialRadius * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256));

 

            double e1 = (1 - Math.Sqrt(1 - eccSquared)) / (1 + Math.Sqrt(1 - eccSquared));

            // phi in radians

            double 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);

            // convert to degrees

            double phi1 = phi1Rad * dCvtRad2Deg;

 

            double N1 = dEquatorialRadius / Math.Sqrt(1 - eccSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad));

            double T1 = Math.Tan(phi1Rad) * Math.Tan(phi1Rad);

            double C1 = eccPrimeSquared * Math.Cos(phi1Rad) * Math.Cos(phi1Rad);

            double R1 = dEquatorialRadius * (1 - eccSquared) / Math.Pow(1 - eccSquared * Math.Sin(phi1Rad) * Math.Sin(phi1Rad), 1.5);

            double D = x / (N1 * dScaleFactor);

 

            // phi in radians

            dLat = 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);

            // convert to degrees

            dLat = dLat * dCvtRad2Deg;

 

            // lon in radians

            dLon = (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);

            // convert to degrees

            dLon = dLongOrigin + dLon * dCvtRad2Deg;

            return (0);

        }

        private void button1_Click(object sender, EventArgs e)

        {

            try

            {

                //cBaseCmnGIS icBaceCmnGIS = new cBaseCmnGIS();

                double UTMNorthing = System.Convert.ToDouble(txtN.Text);

                double UTMEasting = System.Convert.ToDouble(txtE.Text);

                string sUTMZone = txtZone.Text;  // expected format "12N"

                int iZoneNum = Convert.ToInt32(sUTMZone.Substring(0, sUTMZone.Length - 1));

                double dLat = 0;

                double dLon = 0;

                int iRes = iUTM2LatLon(UTMNorthing, UTMEasting, sUTMZone, ref dLat, ref dLon);

 

                double d,m,s;

                d=Math.Floor(dLat);

                m=Math .Floor ((dLat -d)*60);

                s = Math.Round (((dLat - d) * 60-m)*60,4);

                lblfi.Text = "fi= " + d.ToString() + "°  " + m.ToString() +"'  "+s.ToString() +"\"";

                d = Math.Floor(dLon);

                m = Math.Floor((dLon - d) * 60);

                s = Math.Round(((dLon - d) * 60 - m) * 60, 4);

                lblLanda.Text = "Landa: " + d.ToString() + "°  " + m.ToString() + "'  " + s.ToString() + "\"";

                double k1, k2, k;//ScaleFactor

                double dCentralMeridian = (183.0 - (6.0 * (double)iZoneNum)) * -1.0;

                double d1 = dCvtDeg2Rad * (dCentralMeridian - dLon);

                double d2 = Math.Cos(dLat * dCvtDeg2Rad);

                k1 = dScaleFactor * (1 + d1 * d1 * d2 * d2 / 2);

                k2 = dMeanRadius / (dMeanRadius + Convert.ToDouble(txth.Text));

                k = k1 * k2;

                txtProjectionSF.Text = k1.ToString();

                txtEF.Text = k2.ToString();

                txtCombinedSF.Text = k.ToString();

            }

            catch (Exception Err)

            {

                Console.WriteLine("{0} Exception caught.", Err);

            }

        }

    }

}

you can download it here

The precision is tested with Leica Geo Office tools (LGO) and i don't see any problem.

In this little application you enter E,N,h in UTM coordinate system and also your zone. at the end of zone number you  should add N or S for North hemisphere or South hemisphere. then press compute.

The default location that are fill in fields (this is the location of my home ;-)) may be a good guide.

I thank many internet source code and articles. 

If you have any question about similar problems i'll be so happy if i can help.

Blog post

Thank you.

Comment

I have done a survey using DGPS by collecting latlong and leaving behind at site temporary bench mack.

I then converted the lat long (WGS84) coordinates to UTM zone 43R coordinates using scale factor of 0.9996.

These utm cord were then exported to autocad for map making.

The above utm map with benchmark coordinates was supplied to another surveyor using a total station.

He says angle and distances dont mach.

plz suggest what to do....

Comment Comment

How do I calculate utm scale factor between the Central Merdian false easting of 500,000m and the edge of the UTM Zone that has scale factor 1.00000 and also going beyond where the scale factor is 1.0000.

Is there a way of calcuating the UTM scale fator by just giving how far the area of survey work is away from the CM. Also please tell me at which easting value the scale factor is 1.

I am referrig to secant map projection.

Yours gratefully,

Ian Waller

Comment Comment

good

Comment

thank you very much

Comment


Comment


This Free Website was created with the Jigsy Website Builder     Build your own Website for Free!