Home > utilities > wgs2utm.m

wgs2utm

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

 Description:
    Convert WGS84 coordinates (Latitude, Longitude) into UTM coordinates
    (northing, easting) according to (optional) input UTM zone and
    hemisphere.

 Input:
    Lat: WGS84 Latitude scalar, vector or array in decimal degrees.  
    Lon: WGS84 Longitude scalar, vector or array in decimal degrees. 
    utmzone (optional): UTM longitudinal zone. Scalar or same size as Lat
       and Lon.
    utmhemi (optional): UTM hemisphere as a single character, 'N' or 'S',
       or array of 'N' or 'S' characters of same size as Lat and Lon.

 Output:
    x: UTM easting in meters.
    y: UTM northing in meters.
    utmzone: UTM longitudinal zone.
    utmhemi: UTM hemisphere as array of 'N' or 'S' characters. 

 Author notes:
    I downloaded and tried deg2utm.m from Rafael Palacios but found
    differences of up to 1m with my reference converters in southern
    hemisphere so I wrote my own code based on "Map Projections - A
    Working Manual" by J.P. Snyder (1987). Quick quality control performed
    only by comparing with LINZ converter
    (www.linz.govt.nz/apps/coordinateconversions/) and Chuck Taylor's 
    (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html) on a 
    few test points, so use results with caution. Equations not suitable
    for a latitude of +/- 90deg.

    UPDATE: Following requests, this new version allows forcing UTM zone
    in input.

 Examples:

    % set random latitude and longitude arrays
    Lat= 90.*(2.*rand(3)-1)
    Lon= 180.*(2.*rand(3)-1)

    % let the function find appropriate UTM zone and hemisphere from data 
    [x1,y1,utmzone1,utmhemi1] = wgs2utm(Lat,Lon)

    % forcing unique UTM zone and hemisphere for all data entries
    % note: resulting easting and northing are way off the usual values
    [x2,y2,utmzone2,utmhemi2] = wgs2utm(Lat,Lon,60,'S')

    % forcing different UTM zone and hemisphere for each data entry
    % note: resulting easting and northing are way off the usual values
    utmzone = floor(59.*rand(3))+1
    utmhemi = char(78 + 5.*round(rand(3)))
    [x3,y3,utmzone3,utmhemi3] = wgs2utm(Lat,Lon,utmzone,utmhemi)

 Author: 
   Alexandre Schimel
   MetOcean Solutions Ltd
   New Plymouth, New Zealand

 Version 2:
   February 2011
-------------------------------------------------------------------------

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function  [x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,utmhemi)
0002 % -------------------------------------------------------------------------
0003 % [x,y,utmzone] = wgs2utm(Lat,Lon,Zone)
0004 %
0005 % Description:
0006 %    Convert WGS84 coordinates (Latitude, Longitude) into UTM coordinates
0007 %    (northing, easting) according to (optional) input UTM zone and
0008 %    hemisphere.
0009 %
0010 % Input:
0011 %    Lat: WGS84 Latitude scalar, vector or array in decimal degrees.
0012 %    Lon: WGS84 Longitude scalar, vector or array in decimal degrees.
0013 %    utmzone (optional): UTM longitudinal zone. Scalar or same size as Lat
0014 %       and Lon.
0015 %    utmhemi (optional): UTM hemisphere as a single character, 'N' or 'S',
0016 %       or array of 'N' or 'S' characters of same size as Lat and Lon.
0017 %
0018 % Output:
0019 %    x: UTM easting in meters.
0020 %    y: UTM northing in meters.
0021 %    utmzone: UTM longitudinal zone.
0022 %    utmhemi: UTM hemisphere as array of 'N' or 'S' characters.
0023 %
0024 % Author notes:
0025 %    I downloaded and tried deg2utm.m from Rafael Palacios but found
0026 %    differences of up to 1m with my reference converters in southern
0027 %    hemisphere so I wrote my own code based on "Map Projections - A
0028 %    Working Manual" by J.P. Snyder (1987). Quick quality control performed
0029 %    only by comparing with LINZ converter
0030 %    (www.linz.govt.nz/apps/coordinateconversions/) and Chuck Taylor's
0031 %    (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html) on a
0032 %    few test points, so use results with caution. Equations not suitable
0033 %    for a latitude of +/- 90deg.
0034 %
0035 %    UPDATE: Following requests, this new version allows forcing UTM zone
0036 %    in input.
0037 %
0038 % Examples:
0039 %
0040 %    % set random latitude and longitude arrays
0041 %    Lat= 90.*(2.*rand(3)-1)
0042 %    Lon= 180.*(2.*rand(3)-1)
0043 %
0044 %    % let the function find appropriate UTM zone and hemisphere from data
0045 %    [x1,y1,utmzone1,utmhemi1] = wgs2utm(Lat,Lon)
0046 %
0047 %    % forcing unique UTM zone and hemisphere for all data entries
0048 %    % note: resulting easting and northing are way off the usual values
0049 %    [x2,y2,utmzone2,utmhemi2] = wgs2utm(Lat,Lon,60,'S')
0050 %
0051 %    % forcing different UTM zone and hemisphere for each data entry
0052 %    % note: resulting easting and northing are way off the usual values
0053 %    utmzone = floor(59.*rand(3))+1
0054 %    utmhemi = char(78 + 5.*round(rand(3)))
0055 %    [x3,y3,utmzone3,utmhemi3] = wgs2utm(Lat,Lon,utmzone,utmhemi)
0056 %
0057 % Author:
0058 %   Alexandre Schimel
0059 %   MetOcean Solutions Ltd
0060 %   New Plymouth, New Zealand
0061 %
0062 % Version 2:
0063 %   February 2011
0064 %-------------------------------------------------------------------------
0065 
0066 %% Argument checking
0067 if ~sum(double(nargin==[2,4]))
0068     error('Wrong number of input arguments');return
0069 end
0070 n1=size(Lat);
0071 n2=size(Lon);
0072 if (n1~=n2)
0073     error('Lat and Lon should have same size');return
0074 end
0075 if exist('utmzone','var') && exist('utmhemi','var')
0076     n3=size(utmzone);
0077     n4=size(utmhemi);
0078     if (sort(n3)~=sort(n4))
0079         error('utmzone and utmhemi should have same size');return
0080     end
0081     if max(n3)~=1 && max(n3)~=max(n1)
0082         error('utmzone should have either same size as Lat and Long, or size=1');return
0083     end
0084 end
0085 
0086 % expand utmzone and utmhemi if needed
0087 if exist('utmzone','var') && exist('utmhemi','var')
0088     n3=size(utmzone);
0089     n4=size(utmhemi);
0090     if n3==[1 1]
0091         utmzone = utmzone.*ones(size(Lat));
0092         utmhemi = char(utmhemi.*ones(size(Lat)));
0093     end
0094 end
0095     
0096 %% coordinates in radians
0097 lat = Lat.*pi./180;
0098 lon = Lon.*pi./180;
0099 
0100 %% WGS84 parameters
0101 a = 6378137;           %semi-major axis
0102 b = 6356752.314245;    %semi-minor axis
0103 % b = 6356752.314140;  %GRS80 value, originally used for WGS84 before refinements
0104 e = sqrt(1-(b./a).^2); % eccentricity
0105 
0106 %% UTM parameters
0107 % lat0 = 0;                % reference latitude, not used here
0108 if exist('utmzone','var')
0109     Lon0 = 6.*utmzone-183; % reference longitude in degrees
0110 else
0111     Lon0 = floor(Lon./6).*6+3; % reference longitude in degrees
0112 end
0113 lon0 = Lon0.*pi./180;      % in radians
0114 k0 = 0.9996;               % scale on central meridian
0115 
0116 FE = 500000;              % false easting
0117 if exist('utmhemi','var')
0118     FN = double(utmhemi=='S').*10000000;
0119 else
0120     FN = (Lat < 0).*10000000; % false northing
0121 end
0122 
0123 %% Equations parameters
0124 eps = e.^2./(1-e.^2);  % e prime square
0125 % N: radius of curvature of the earth perpendicular to meridian plane
0126 % Also, distance from point to polar axis
0127 N = a./sqrt(1-e.^2.*sin(lat).^2); 
0128 T = tan(lat).^2;                
0129 C = ((e.^2)./(1-e.^2)).*(cos(lat)).^2;
0130 A = (lon-lon0).*cos(lat);                            
0131 % M: true distance along the central meridian from the equator to lat
0132 M = a.*(  ( 1 - e.^2./4 - 3.*e.^4./64 - 5.*e.^6./256 )  .* lat         ...
0133          -( 3.*e.^2./8 + 3.*e.^4./32 + 45.*e.^6./1024 ) .* sin(2.*lat) ...
0134          +( 15.*e.^4./256 + 45.*e.^6./1024 )            .* sin(4.*lat) ...
0135          -(35.*e.^6./3072 )                             .* sin(6.*lat) );
0136 
0137 %% easting
0138 x = FE + k0.*N.*(                                  A       ...
0139                  + (1-T+C)                      .* A.^3./6 ...
0140                  + (5-18.*T+T.^2+72.*C-58.*eps) .* A.^5./120 );
0141                
0142 %% northing
0143 % M(lat0) = 0 so not used in following formula
0144 y = FN + k0.*M + k0.*N.*tan(lat).*(                                     A.^2./2  ...
0145                                    + (5-T+9.*C+4.*C.^2)              .* A.^4./24 ...
0146                                    + (61-58.*T+T.^2+600.*C-330.*eps) .* A.^6./720 );
0147                                  
0148 %% UTM zone
0149 if exist('utmzone','var') && exist('utmhemi','var')
0150     utmzone = utmzone;
0151     utmhemi = utmhemi;
0152 else 
0153    utmzone = floor(Lon0./6)+31;
0154    utmhemi = char( 83.* (Lat < 0) + 78.* (Lat >= 0) );
0155 end

Generated on Wed 20-Feb-2019 16:06:01 by m2html © 2005