utm2deg

PURPOSE ^

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

SYNOPSIS ^

function [Lat,Lon] = utm2deg(xx,yy,utmzone)

DESCRIPTION ^

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

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

 Inputs:
    x, y , utmzone.

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

 Example 1:
 x=[ 458731;  407653;  239027;  230253;  343898;  362850];
 y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478];
 utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R'];
 [Lat, Lon]=utm2deg(x,y,utmzone);
 fprintf('%11.6f ',lat)
    40.315430   46.283902   37.577834   28.645647   38.855552   25.061780
 fprintf('%11.6f ',lon)
    -3.485713    7.801235 -119.955246  -17.759537  -94.799019  121.640266

 Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds
 [Lat, Lon]=utm2deg(x,y,utmzone);
 LatDMS=dms2mat(deg2dms(Lat))
LatDMS =
    40.00         18.00         55.55
    46.00         17.00          2.01
    37.00         34.00         40.17
    28.00         38.00         44.33
    38.00         51.00         19.96
    25.00          3.00         42.41
 LonDMS=dms2mat(deg2dms(Lon))
LonDMS =
    -3.00         29.00          8.61
     7.00         48.00          4.40
  -119.00         57.00         18.93
   -17.00         45.00         34.33
   -94.00         47.00         56.47
   121.00         38.00         24.96

 Author: 
   Rafael Palacios
   Universidad Pontificia Comillas
   Madrid, Spain
 Version: Apr/06, Jun/06, Aug/06
 Aug/06: corrected m-Lint warnings
-------------------------------------------------------------------------

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function  [Lat,Lon] = utm2deg(xx,yy,utmzone)
0002 % -------------------------------------------------------------------------
0003 % [Lat,Lon] = utm2deg(x,y,utmzone)
0004 %
0005 % Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84).
0006 % Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez.
0007 %
0008 % Inputs:
0009 %    x, y , utmzone.
0010 %
0011 % Outputs:
0012 %    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
0013 %    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
0014 %
0015 % Example 1:
0016 % x=[ 458731;  407653;  239027;  230253;  343898;  362850];
0017 % y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478];
0018 % utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R'];
0019 % [Lat, Lon]=utm2deg(x,y,utmzone);
0020 % fprintf('%11.6f ',lat)
0021 %    40.315430   46.283902   37.577834   28.645647   38.855552   25.061780
0022 % fprintf('%11.6f ',lon)
0023 %    -3.485713    7.801235 -119.955246  -17.759537  -94.799019  121.640266
0024 %
0025 % Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds
0026 % [Lat, Lon]=utm2deg(x,y,utmzone);
0027 % LatDMS=dms2mat(deg2dms(Lat))
0028 %LatDMS =
0029 %    40.00         18.00         55.55
0030 %    46.00         17.00          2.01
0031 %    37.00         34.00         40.17
0032 %    28.00         38.00         44.33
0033 %    38.00         51.00         19.96
0034 %    25.00          3.00         42.41
0035 % LonDMS=dms2mat(deg2dms(Lon))
0036 %LonDMS =
0037 %    -3.00         29.00          8.61
0038 %     7.00         48.00          4.40
0039 %  -119.00         57.00         18.93
0040 %   -17.00         45.00         34.33
0041 %   -94.00         47.00         56.47
0042 %   121.00         38.00         24.96
0043 %
0044 % Author:
0045 %   Rafael Palacios
0046 %   Universidad Pontificia Comillas
0047 %   Madrid, Spain
0048 % Version: Apr/06, Jun/06, Aug/06
0049 % Aug/06: corrected m-Lint warnings
0050 %-------------------------------------------------------------------------
0051 format long
0052 % Argument checking
0053 %
0054 error(nargchk(3, 3, nargin)); %3 arguments required
0055 n1=length(xx);
0056 n2=length(yy);
0057 n3=size(utmzone,1);
0058 if (n1~=n2 || n1~=n3)
0059    error('x,y and utmzone vectors should have the same number or rows');
0060 end
0061 c=size(utmzone,2);
0062 if (c~=4)
0063    error('utmzone should be a vector of strings like "30 T"');
0064 end
0065 
0066    
0067  
0068 % Memory pre-allocation
0069 %
0070 Lat=zeros(n1,1);
0071 Lon=zeros(n1,1);
0072 
0073 
0074 % Main Loop
0075 %
0076 for i=1:n1
0077    if (utmzone(i,4)>'X' || utmzone(i,4)<'C')
0078       fprintf('utm2deg: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n');
0079    end
0080    if (utmzone(i,4)>'M')
0081       hemis='N';   % Northern hemisphere
0082    else
0083       hemis='S';
0084    end
0085 
0086    x=xx(i);
0087    y=yy(i);
0088    zone=str2double(utmzone(i,1:2));
0089 
0090    sa = 6378137.000000 ; sb = 6356752.314245;
0091   
0092 %   e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
0093    e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
0094    e2cuadrada = e2 ^ 2;
0095    c = ( sa ^ 2 ) / sb;
0096 %   alpha = ( sa - sb ) / sa;             %f
0097 %   ablandamiento = 1 / alpha;   % 1/f
0098 
0099    X = x - 500000;
0100    
0101    if hemis == 'S' || hemis == 's'
0102        Y = y - 10000000;
0103    else
0104        Y = y;
0105    end
0106     
0107    S = ( ( zone * 6 ) - 183 ); 
0108    lat =  Y / ( 6366197.724 * 0.9996 );                                    
0109    v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
0110    a = X / v;
0111    a1 = sin( 2 * lat );
0112    a2 = a1 * ( cos(lat) ) ^ 2;
0113    j2 = lat + ( a1 / 2 );
0114    j4 = ( ( 3 * j2 ) + a2 ) / 4;
0115    j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
0116    alfa = ( 3 / 4 ) * e2cuadrada;
0117    beta = ( 5 / 3 ) * alfa ^ 2;
0118    gama = ( 35 / 27 ) * alfa ^ 3;
0119    Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
0120    b = ( Y - Bm ) / v;
0121    Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2;
0122    Eps = a * ( 1 - ( Epsi / 3 ) );
0123    nab = ( b * ( 1 - Epsi ) ) + lat;
0124    senoheps = ( exp(Eps) - exp(-Eps) ) / 2;
0125    Delt = atan(senoheps / (cos(nab) ) );
0126    TaO = atan(cos(Delt) * tan(nab));
0127    longitude = (Delt *(180 / pi ) ) + S;
0128    latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ...
0129                     (180 / pi);
0130    
0131    Lat(i)=latitude;
0132    Lon(i)=longitude;
0133    
0134 end

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