parameters and functions for UTMLat/Long conversions

Posted: September 3, 2011 in Java, Mapping
/*---------------------------------------------------------*/ /* parameters and functions for UTM<->Lat/Long conversions */ /* UTM conversion is based on WGS84 ellipsoid parameters */ /*---------------------------------------------------------*/

RADIUS = 6378137.0;
FLATTENING = 0.00335281068; /* GRS80 or WGS84 */
K_NOT = 0.9996;     /* UTM scale factor */
DEGREES_TO_RADIANS = 0.01745329252;
FALSE_EASTING = 500000.0;
FALSE_NORTHING = 10000000.0;

/*--------------------------------------------------------*/
/* These Java functions were converted from original      */
/* C fucntions provided by Dana Yoerger, Steve Gegg,      */
/* and Louis Whitcomb at WHOI                             */
/*--------------------------------------------------------*/

function geo_utm(lat, lon, zone)
{
   with(Math)
   {
      /* first compute the necessary geodetic parameters and constants */

      lambda_not = ((-180.0 + zone*6.0) -3.0)/RADIANS ;
      e_squared = 2.0 * FLATTENING - FLATTENING*FLATTENING;
      e_fourth = e_squared * e_squared;
      e_sixth = e_fourth * e_squared;
      e_prime_sq = e_squared/(1.0 - e_squared);
      sin_phi = sin(lat);
      tan_phi = tan(lat);
      cos_phi = cos(lat);
      N = RADIUS/sqrt(1.0 - e_squared*sin_phi*sin_phi);
      T = tan_phi*tan_phi;
      C = e_prime_sq*cos_phi*cos_phi;
      M = RADIUS*((1.0 - e_squared*0.25 -0.046875*e_fourth  -0.01953125*e_sixth)*lat-
	      (0.375*e_squared + 0.09375*e_fourth +
				 0.043945313*e_sixth)*sin(2.0*lat) +
	      (0.05859375*e_fourth + 0.043945313*e_sixth)*sin(4.0*lat) -
	      (0.011393229 * e_sixth)*sin(6.0*lat));
      A = (lon - lambda_not)*cos_phi;
      A_sq = A*A;
      A_fourth =  A_sq*A_sq;

      /* now go ahead and compute X and Y */

      x_utm = K_NOT*N*(A + (1.0 - T + C)*A_sq*A/6.0 +
		   (5.0 - 18.0*T + T*T + 72.0*C -
		    58.0*e_prime_sq)*A_fourth*A/120.0);

      /* note:  specific to UTM, vice general trasverse mercator.
         since the origin is at the equator, M0, the M at phi_0,
         always equals zero, and I won't compute it   */                                            

       y_utm = K_NOT*(M + N*tan_phi*(A_sq/2.0 +
			    (5.0 - T + 9.0*C + 4.0*C*C)*A_fourth/24.0 +
			    (61.0 - 58.0*T + T*T + 600.0*C -
			     330.0*e_prime_sq)*A_fourth*A_sq/720.0));

       /* now correct for false easting and northing */

       if( lat < 0)
       {
          y_utm +=10000000.0;
       }
       x_utm +=500000;
    }
    return true;
}

function utm_geo(x_utm, y_utm, zone)
{
   with(Math)
   {
      /* first, subtract the false easting */
      x_utm = x_utm - FALSE_EASTING;

      /* compute the necessary geodetic parameters and constants */

      e_squared = 2.0 * FLATTENING - FLATTENING*FLATTENING;
      e_fourth = e_squared * e_squared;
      e_sixth = e_fourth * e_squared;
      oneminuse = sqrt(1.0-e_squared);  

      /* compute the footpoint latitude */

      M = y_utm/K_NOT;
      mu =M/(RADIUS*(1.0 - 0.25*e_squared -
                  0.046875*e_fourth - 0.01953125*e_sixth));
      e1 = (1.0 - oneminuse)/(1.0 + oneminuse);
      e1sq =e1*e1;
      footpoint = mu + (1.5*e1 - 0.84375*e1sq*e1)*sin(2.0*mu) +
              (1.3125*e1sq - 1.71875*e1sq*e1sq)*sin(4.0*mu) +
              (1.57291666667*e1sq*e1)*sin(6.0*mu) +
              (2.142578125*e1sq*e1sq)*sin(8.0*mu);

      /* compute the other necessary terms */

      e_prime_sq =  e_squared/(1.0 -  e_squared);
      sin_phi = sin(footpoint);
      tan_phi = tan(footpoint);
      cos_phi = cos(footpoint);
      N = RADIUS/sqrt(1.0 - e_squared*sin_phi*sin_phi);
      T = tan_phi*tan_phi;
      Tsquared = T*T;
      C = e_prime_sq*cos_phi*cos_phi;
      Csquared = C*C;
      denom = sqrt(1.0-e_squared*sin_phi*sin_phi);
      R = RADIUS*oneminuse*oneminuse/(denom*denom*denom);
      D = x_utm/(N*K_NOT);
      Dsquared = D*D;
      Dfourth = Dsquared*Dsquared;

      lambda_not = ((-180.0 + zone*6.0) -3.0) * DEGREES_TO_RADIANS;

      /* now, use the footpoint to compute the real latitude and longitude */

      var lat = footpoint - (N*tan_phi/R)*(0.5*Dsquared - (5.0 + 3.0*T + 10.0*C -
                           4.0*Csquared - 9.0*e_prime_sq)*Dfourth/24.0 +
                           (61.0 + 90.0*T + 298.0*C + 45.0*Tsquared -
                            252.0*e_prime_sq -
                            3.0*Csquared)*Dfourth*Dsquared/720.0);
      var lon = lambda_not + (D - (1.0 + 2.0*T + C)*Dsquared*D/6.0 +
                         (5.0 - 2.0*C + 28.0*T - 3.0*Csquared + 8.0*e_prime_sq +
                          24.0*Tsquared)*Dfourth*D/120.0)/cos_phi;

   }
   return true;
}

Reference:
http://www.whoi.edu/marine/ndsf/utility/NDSFref.txt
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s