wgs2utm.html 15.4 KB
Newer Older
Geoffrey Cowles's avatar
Geoffrey Cowles committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of wgs2utm</title>
  <meta name="keywords" content="wgs2utm">
  <meta name="description" content="-------------------------------------------------------------------------">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">utilities</a> &gt; wgs2utm.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for utilities&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>wgs2utm
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>-------------------------------------------------------------------------</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function  [x,y,utmzone,utmhemi] = wgs2utm(Lat,Lon,utmzone,utmhemi) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment"> -------------------------------------------------------------------------
 [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 &quot;Map Projections - A
    Working Manual&quot; 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
-------------------------------------------------------------------------</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
101
<li><a href="get_polcoms_timeseriesv1.html" class="code" title="function [polcoms]=get_polcoms_timeseriesv1(rootfname,ts_controlfile,Mobj,inputConf,tseries_dir,mm,polcoms)">get_polcoms_timeseriesv1</a>	ts_controlfile='/users/modellers/rito/Models/MEDINA/tseries.MEDI29'</li></ul>
Geoffrey Cowles's avatar
Geoffrey Cowles committed
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
<!-- crossreference -->



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