------------------------------------------------------------------------- [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 -------------------------------------------------------------------------
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