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