deg2utm

PURPOSE ^

-------------------------------------------------------------------------

SYNOPSIS ^

function [x,y,utmzone] = deg2utm(Lat,Lon)

DESCRIPTION ^

 -------------------------------------------------------------------------
 [x,y,utmzone] = deg2utm(Lat,Lon)

 Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
 Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.

 Inputs:
    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84

 Outputs:
    x, y , utmzone.   See example

 Example 1:
    Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
    Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
    [x,y,utmzone] = deg2utm(Lat,Lon);
    fprintf('%7.0f ',x)
       458731  407653  239027  230253  343898  362850
    fprintf('%7.0f ',y)
      4462881 5126290 4163083 3171843 4302285 2772478
    utmzone =
       30 T
       32 T
       11 S
       28 R
       15 S
       51 R

 Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
    LatDMS=[40 18 55.56; 46 17 2.04];
    LonDMS=[-3 29  8.58;  7 48 4.44];
    Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
    Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
    [x,y,utmzone] = deg2utm(Lat,Lon)

 Author: 
   Rafael Palacios
   Universidad Pontificia Comillas
   Madrid, Spain
 Version: Apr/06, Jun/06, Aug/06, Aug/06
 Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern 
    hemisphere coordinates. 
 Aug/06: corrected m-Lint warnings
-------------------------------------------------------------------------

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function  [x,y,utmzone] = deg2utm(Lat,Lon)
0002 % -------------------------------------------------------------------------
0003 % [x,y,utmzone] = deg2utm(Lat,Lon)
0004 %
0005 % Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
0006 % Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.
0007 %
0008 % Inputs:
0009 %    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
0010 %    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
0011 %
0012 % Outputs:
0013 %    x, y , utmzone.   See example
0014 %
0015 % Example 1:
0016 %    Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
0017 %    Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
0018 %    [x,y,utmzone] = deg2utm(Lat,Lon);
0019 %    fprintf('%7.0f ',x)
0020 %       458731  407653  239027  230253  343898  362850
0021 %    fprintf('%7.0f ',y)
0022 %      4462881 5126290 4163083 3171843 4302285 2772478
0023 %    utmzone =
0024 %       30 T
0025 %       32 T
0026 %       11 S
0027 %       28 R
0028 %       15 S
0029 %       51 R
0030 %
0031 % Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
0032 %    LatDMS=[40 18 55.56; 46 17 2.04];
0033 %    LonDMS=[-3 29  8.58;  7 48 4.44];
0034 %    Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
0035 %    Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
0036 %    [x,y,utmzone] = deg2utm(Lat,Lon)
0037 %
0038 % Author:
0039 %   Rafael Palacios
0040 %   Universidad Pontificia Comillas
0041 %   Madrid, Spain
0042 % Version: Apr/06, Jun/06, Aug/06, Aug/06
0043 % Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern
0044 %    hemisphere coordinates.
0045 % Aug/06: corrected m-Lint warnings
0046 %-------------------------------------------------------------------------
0047 
0048 % Argument checking
0049 %
0050 error(nargchk(2, 2, nargin));  %2 arguments required
0051 n1=length(Lat);
0052 n2=length(Lon);
0053 if (n1~=n2)
0054    error('Lat and Lon vectors should have the same length');
0055 end
0056 
0057 
0058 % Memory pre-allocation
0059 %
0060 x=zeros(n1,1);
0061 y=zeros(n1,1);
0062 utmzone(n1,:)='60 X';
0063 
0064 % Main Loop
0065 %
0066 for i=1:n1
0067    la=Lat(i);
0068    lo=Lon(i);
0069 
0070    sa = 6378137.000000 ; sb = 6356752.314245;
0071          
0072    %e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
0073    e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
0074    e2cuadrada = e2 ^ 2;
0075    c = ( sa ^ 2 ) / sb;
0076    %alpha = ( sa - sb ) / sa;             %f
0077    %ablandamiento = 1 / alpha;   % 1/f
0078 
0079    lat = la * ( pi / 180 );
0080    lon = lo * ( pi / 180 );
0081 
0082    Huso = fix( ( lo / 6 ) + 31);
0083    S = ( ( Huso * 6 ) - 183 );
0084    deltaS = lon - ( S * ( pi / 180 ) );
0085 
0086    if (la<-72), Letra='C';
0087    elseif (la<-64), Letra='D';
0088    elseif (la<-56), Letra='E';
0089    elseif (la<-48), Letra='F';
0090    elseif (la<-40), Letra='G';
0091    elseif (la<-32), Letra='H';
0092    elseif (la<-24), Letra='J';
0093    elseif (la<-16), Letra='K';
0094    elseif (la<-8), Letra='L';
0095    elseif (la<0), Letra='M';
0096    elseif (la<8), Letra='N';
0097    elseif (la<16), Letra='P';
0098    elseif (la<24), Letra='Q';
0099    elseif (la<32), Letra='R';
0100    elseif (la<40), Letra='S';
0101    elseif (la<48), Letra='T';
0102    elseif (la<56), Letra='U';
0103    elseif (la<64), Letra='V';
0104    elseif (la<72), Letra='W';
0105    else Letra='X';
0106    end
0107 
0108    a = cos(lat) * sin(deltaS);
0109    epsilon = 0.5 * log( ( 1 +  a) / ( 1 - a ) );
0110    nu = atan( tan(lat) / cos(deltaS) ) - lat;
0111    v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
0112    ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2;
0113    a1 = sin( 2 * lat );
0114    a2 = a1 * ( cos(lat) ) ^ 2;
0115    j2 = lat + ( a1 / 2 );
0116    j4 = ( ( 3 * j2 ) + a2 ) / 4;
0117    j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
0118    alfa = ( 3 / 4 ) * e2cuadrada;
0119    beta = ( 5 / 3 ) * alfa ^ 2;
0120    gama = ( 35 / 27 ) * alfa ^ 3;
0121    Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
0122    xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;
0123    yy = nu * v * ( 1 + ta ) + Bm;
0124 
0125    if (yy<0)
0126        yy=9999999+yy;
0127    end
0128 
0129    x(i)=xx;
0130    y(i)=yy;
0131    if ~isnan(Lat(i))
0132        utmzone(i,:)=sprintf('%02d %c',Huso,Letra);
0133    else
0134        utmzone(i,:)='';
0135    end
0136 end
0137

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005