As the first step in the decommissioning of sasCommunity.org the site has been converted to read-only mode.


Here are some tips for How to share your SAS knowledge with your professional network.


Latitude/longitude to UTM conversion (and vice-versa)

From sasCommunity
Jump to: navigation, search

Geographic coordinates are most often reported in latitude/longitude or in Universal Transverse Mercator (UTM) coordinates. The following code converts between the two. It has been validated using Steven Dutch's web page [1] which offers both Microsoft Excel and javascript conversion tools.

The mathematics of the conversion is complicated. Latitude and longitude coordinates are in relation to the surface of an (approximate) sphere while UTM flattens the earth into 120 distinct zones with minimal distortion in each zone - two very different ways of representing location. For precise details see:

Army, Department of, 1973. Universal Transverse Mercator Grid, U. S. Army Technical Manual TM 5-241-8, 64 p.

NIMA Technical Report 8350.2, "Department of Defense World Geodetic System 1984, Its Definition and Relationships with Local Geodetic Systems," Second Edition, 1 September 1991 and its supplements.

In addition, several SAS conference presentations in the 1980s tackled this same problem. [2][3][4]


***************************************************************************
Convert utm to lat/lon
***************************************************************************;
 
data utm2lat;
input place & : $20. hemisphere $ zone northing easting;
datalines;
Eiffel Tower  N 31 5411949 448231   
Sydney Opera House  S 56 6252309 334897 
Buenos Aires Obelisk  S 21 6170000 373315 
Statue of Liberty  N 18 4504682 580720  
;
run;
%let a=6378137; *equatorial radius;
%let b=6356752.314; *polar radius (values of a and b are for WGS84 datum - for other datums see http://www.uwgb.edu/dutchs/usefuldata/utmformulas.htm );
%let k0=0.9996; *scale along central meridian;
 
%let e=%sysevalf((1-(&b**2/&a**2))**0.5); *eccentricity;
%let e1=%sysevalf((1-(1-&e**2)**0.5)/(1+(1-&e**2)**0.5));
%let ep2=%sysevalf(&e**2/(1-&e**2)); *e prime squared;
%let c1=%sysevalf((3*&e1/2)-(27*(&e1**3)/32));
%let c2=%sysevalf((21*(&e1**2)/16)-(55*(&e1**4)/32));
%let c3=%sysevalf((151*(&e1**3)/96));
%let c4=%sysevalf((1097*(&e1**4)/512));
 
 
*Intermediate variables below follow the use by Steven Dutch. CN = corrected northing, EP = east prime, AL = arc length, 
mu and phi are greek letters, others have no specific reference;
 
data utm2lat (drop=CN EP AL mu phi C1 T1 N1 R1 D Q1-Q7 DL);
set utm2lat;
if hemisphere="N" then CN=northing;
else CN=10000000-northing;
EP=500000-easting;
AL=CN/&k0;
mu=AL/(&a*((1-&e**2/4)-(3*&e**4/64)-(5*&e**6/256)));
phi=mu+&c1*sin(2*mu)+&c2*sin(4*mu)+&c3*sin(6*mu)+&c4*sin(8*mu);
C1=&ep2*cos(phi)**2;
T1=tan(phi)**2;
N1=&a/(1-(&e*sin(phi))**2)**0.5;
R1=&a*(1-&e*&e)/(1-(&e*sin(phi))**2)**1.5;
D=EP/(N1*&k0);
Q1=N1*tan(phi)/R1;
Q2=D**2/2;
Q3=(5+(3*T1)+(10*C1)-(4*C1**2)-(9*&ep2))*(D**4)/24;
Q4=(61+(90*T1)+(298*C1)+(45*T1**2)-(252*&ep2)-(3*C1**2))*(D**6)/720;
Q5=D;
Q6=(1+(2*T1)+C1)*(D**3)/6;
Q7=(5-(2*C1)+(28*T1)-(3*C1**2)+(8*&ep2)+(24*T1**2))*(D**5)/120;
DL=(Q5-Q6+Q7)/cos(phi);
latitude=180*(phi-Q1*(Q2+Q3+Q4))/constant('pi');
if hemisphere="S" then latitude=latitude*-1;; 
longitude=((6*zone)-183)-DL*180/constant('pi');
run;
 
***************************************************************************
Convert lat/lon to utm
***************************************************************************;
 
data lat2utm;
input place & : $20. lat lon;
datalines;
Eiffel Tower  48.8583 2.2942   
Sydney Opera House  -33.8566 151.2153 
Buenos Aires Obelisk  -34.6040 -58.3816 
Statue of Liberty  40.6890 -74.0447 
;
run;
%let a=6378137; *equatorial radius (default is for WGS84 datum);
%let b=6356752.314; *polar radius (default is for WGS84 datum);
%let k0=0.9996; *scale along central meridian;
%let e=%sysevalf((1-(&b**2/&a**2))**0.5); *eccentricity;
%let ep2=%sysevalf(&e**2/(1-&e**2)); *e prime squared;
%let n=%sysevalf((&a-&b)/(&a+&b));
%let A0=%sysevalf(&a*((1-&n)+((1-&n)*(5*&n**2)/4)+((1-&n)*(81/64*&n**4))));
%let B0=%sysevalf((3/2*&a*&n)*((1-&n)-((1-&n)*7/8*&n**2)+(55/64*&n**4)));
%let C0=%sysevalf((15/16*&a*&n**2)*((1-&n)+((1-&n)*(3/4*&n**2))));
%let D0=%sysevalf((35/48*&a*&n**3)*((1-&n)+(11/16*&n**2)));
%let E0=%sysevalf((315/51*&a*&n**4)*(1-&n));
 
*LM = longitude central meridian, MA = meridianal arc, RN = raw northing
*rho, nu are greek letters;
*all others are intermediate unnamed constants;;
 
data lat2utm (drop=LM deltalon latr lonr rho nu MA Ki Kii Kiii Kiv Kv A6 RN);
set lat2utm;
if lat lt 0 then hemisphere = "S";
else hemisphere = "N";
zone=31+floor(lon/6);
LM=6*zone-183;
deltalon=(lon-LM)*constant('pi')/180;
latr = lat*constant('pi')/180;
lonr = lon*constant('pi')/180;
rho = &a*(1-&e**2)/((1-(&e*sin(latr))**2)**1.5);
nu = &a/((1-(&e*sin(latr))**2)**0.5);
MA = (&A0*latr)-(&B0*sin(2*latr))+(&C0*sin(4*latr))-(&D0*sin(6*latr))+(&E0*sin(8*latr));
Ki = MA*&k0;
Kii = nu*sin(latr)*cos(latr)/2;
Kiii = ((nu*sin(latr)*cos(latr)**3)/24)*(5-tan(latr)**2+9*&ep2*cos(latr)**2+4*&ep2**2*cos(latr)**4)*&k0;
Kiv = nu*cos(latr)*&k0;
Kv = cos(latr)**3*nu/6*(1-tan(latr)**2+&ep2*cos(latr)**2)*&k0;
A6 = ((deltalon)**6*nu*sin(latr)*cos(latr)**5/720)*(61-58*tan(latr)**2+tan(latr)**4+270*&ep2*cos(latr)**2-330*&ep2*sin(latr)**2)*&k0;
RN =(Ki+Kii*deltalon**2+Kiii*deltalon**4);
if RN < 0 then northing=10000000+RN;
else northing=RN;
easting = 500000+(Kiv*deltalon+Kv*deltalon**3);
run;