function [x,y,utmzone] = wgs2utm(Lat,Lon) % This function converts the vectors of Lat/Lon coordinates to those of UTM. % Inputs: % Lat (WGS84 Latitude vector) in decimal degrees % Lon (WGS84 Longitude vector) in decimal degrees % Outputs: % x - UTM easting in meters % y - UTM northing in meters % utmzone - UTM longitudinal zone % % Example: % Lat=[48.866667; 34.05; -36.85]; % Lon=[2.333056; -118.25; 174.783333]; % [x,y,utmzone]=wgs2utm(Lat,Lon) % returns % x = % 1.0e+005 * % 4.5109 % 3.8463 % 3.0237 % y = % 1.0e+006 * % 5.4128 % 3.7684 % 5.9195 % utmzone = % 31U % 11S % 60H % % Source: DMA Technical Manual 8358.2, Fairfax, VA %% Argument checking error(nargchk(2,2,nargin)); % 2 arguments are required n1=size(Lat); n2=size(Lon); if (n1~=n2) error('Lat and Lon should have same size'); return end %% Converting coordinates to radians lat = Lat*pi/180; lon = Lon*pi/180; %% WGS84 parameters a = 6378137; % semi-major axis of the Earth ellipsoid b = 6356752.314245; % semi-minor axis of the Earth ellipsoid e = sqrt(1-(b/a)^2); % first eccentricity %% UTM parameters Lon0 = floor(Lon/6)*6+3; % reference longitude in degrees lon0 = Lon0.*pi/180; % reference longitude in radians k0 = 0.9996; % scale on central meridian FE = 500000; % false easting FN = (Lat < 0).*10000000; % false northing %% Equations parameters eps = e^2/(1-e^2); % squared second eccentricity % N is the radius of curvature of the earth perpendicular to meridian plane % Also, distance from a point to polar axis N = a./sqrt(1-e^2*sin(lat).^2); T = tan(lat).^2; C = eps*(cos(lat)).^2; A = (lon-lon0).*cos(lat); % M: true distance along the central meridian from the equator to lat M = a*( (1 - e^2/4 - 3*e^4/64 - 5*e^6/256)* lat ... -(3*e^2/8 + 3*e^4/32 + 45*e^6/1024)* sin(2*lat) ... +(15*e^4/256 + 45*e^6/1024) * sin(4*lat) ... -(35*e^6/3072 ) * sin(6*lat)); %% Computing easting x = FE + k0*N.*( A ... +(1-T+C) .* A.^3/6 ... +(5-18*T+T.^2+72*C-58*eps) .* A.^5/120); %% Computing northing y = FN + k0*M + k0*N.*tan(lat).*( A.^2/2 ... +(5-T+9*C+4*C.^2) .* A.^4/24 ... +(61-58*T+T.^2+600*C-330*eps) .* A.^6/720); %% Determining the UTM zone lonzone = floor(Lon0./6)+31; latband = floor((Lat+80)./8)+1; LB='CDEFGHJKLMNPQRSTUVWX'; latband=LB(latband)'; utmzone = [num2str(lonzone,'%02g') latband];