m_ll2ll.m 5.34 KB
Newer Older
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
101
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
function [X,Y]=m_ll2ll(lon1,lat1, zone, hemisphere,ellipsoid);
% M_ll2LL Converts LON,LAT to long,lat coordinates using the current projection
% and a provided one. I really only want to do UTM to LATLON so keep it
% simple (this is provided to use m_map instead of proj.exe )
global MAP_PROJECTION MAP_VAR_LIST

% define a structure of various ellipsoids. each has a name, and
% a vector consisting of equatorial radius and flattening. the first
% two are somewhat special cases.

MAP_ELLIP = struct ( 'normal', [1.0, 0], ...
    'sphere', [6370997.0, 0], ...
    'grs80' , [6378137.0, 1/298.257], ...
    'grs67' , [6378160.0, 1/247.247], ...
    'wgs84' , [6378137.0, 1/298.257], ...
    'wgs72' , [6378135.0, 1/298.260], ...
    'wgs66' , [6378145.0, 1/298.250], ...
    'wgs60' , [6378165.0, 1/298.300], ...
    'clrk66', [6378206.4, 1/294.980], ...
    'clrk80', [6378249.1, 1/293.466], ...
    'intl24', [6378388.0, 1/297.000], ...
    'intl67', [6378157.5, 1/298.250]);

MAP_VAR_LIST.zone
MAP_VAR_LIST.hemisphere
MAP_VAR_LIST.ellipsoid
getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid)
%     [X,Y] = mu_ll2utm(lat1,lon1,MAP_VAR_LIST.zone,MAP_VAR_LIST.hemisphere, ...
% 	getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid));
    [Y,X] = mu_utm2ll(lon1, lat1, MAP_VAR_LIST.zone, ...
	MAP_VAR_LIST.hemisphere, getfield(MAP_ELLIP,MAP_VAR_LIST.ellipsoid));

return
function [x,y] = mu_ll2utm (lat,lon, zone, hemisphere,ellipsoid)
%mu_ll2utm		Convert geodetic lat,lon to X/Y UTM coordinates
%
%	[x,y] = mu_ll2utm (lat, lon, zone, hemisphere,ellipsoid)
%
%	input is latitude and longitude vectors, zone number, 
%		hemisphere(N=0,S=1), ellipsoid info [eq-rad, flat]
%	output is X/Y vectors
%
%	see also	mu_utm2ll, utmzone


% some general constants

DEG2RADS    = 0.01745329252;
RADIUS      = ellipsoid(1);
FLAT        = ellipsoid(2);
K_NOT       = 0.9996;
FALSE_EAST  = 500000;
FALSE_NORTH = 10000000;

% check for valid numbers

if (max(abs(lat)) > 90)
  error('latitude values exceed 89 degree');
  return;
end

if ((zone < 1) | (zone > 60))
  error ('utm zones only valid from 1 to 60');
  return;
end

% compute some geodetic parameters

lambda_not  = ((-180 + zone*6) - 3) * DEG2RADS;

e2  = 2*FLAT - FLAT*FLAT;
e4  = e2 * e2;
e6  = e4 * e2;
ep2 = e2/(1-e2);

% some other constants, vectors

lat = lat * DEG2RADS;
lon = lon * DEG2RADS;

sinL = sin(lat);
tanL = tan(lat);
cosL = cos(lat);

T = tanL.*tanL;
C = ep2 * (cosL.*cosL);
A = (lon - lambda_not).*cosL;
A2 = A.*A;
A4 = A2.*A2;
S = sinL.*sinL;

% solve for N

N = RADIUS ./ (sqrt (1-e2*S));

% solve for M

M0 = 1 - e2*0.25 - e4*0.046875 - e6*0.01953125;
M1 = e2*0.375 + e4*0.09375 + e6*0.043945313;
M2 = e4*0.05859375 + e6*0.043945313;
M3 = e6*0.011393229;
M = RADIUS.*(M0.*lat - M1.*sin(2*lat) + M2.*sin(4*lat) - M3.*sin(6*lat));

% solve for x

X0 = A4.*A/120;
X1 = 5 - 18*T + T.*T + 72*C - 58*ep2;
X2 = A2.*A/6;
X3 = 1 - T + C;
x = N.*(A + X3.*X2 + X1.* X0);

% solve for y

Y0 = 61 - 58*T + T.*T + 600*C - 330*ep2;
Y1 = 5 - T + 9*C + 4*C.*C;

y = M + N.*tanL.*(A2/2 + Y1.*A4/24 + Y0.*A4.*A2/720);


% finally, do the scaling and false thing. if using a unit-normal radius,
% we don't bother.

x = x*K_NOT + (RADIUS>1) * FALSE_EAST;

y = y*K_NOT;
if (hemisphere)
  y = y + (RADIUS>1) * FALSE_NORTH;
end

return



%-------------------------------------------------------------------

function [lat,lon] = mu_utm2ll (x,y, zone, hemisphere,ellipsoid)
%mu_utm2ll		Convert X/Y UTM coordinates to geodetic lat,lon 
%
%	[lat,lon] = mu_utm2ll (x,y, zone, hemisphere,ellipsoid)
%
%	input is X/Y vectors, zone number, hemisphere(N=0,S=1),
%		ellipsoid info [eq-rad, flat]
%	output is lat/lon vectors
%
%	see also	mu_ll2utm, utmzone


% some general constants

DEG2RADS    = 0.01745329252;
RADIUS      = ellipsoid(1);
FLAT        = ellipsoid(2);
K_NOT       = 0.9996;
FALSE_EAST  = 500000;
FALSE_NORTH = 10000000;

if ((zone < 1) | (zone > 60))
  error ('utm zones only valid from 1 to 60');
  return;
end

% compute some geodetic parameters

e2  = 2*FLAT - FLAT*FLAT;
e4  = e2 * e2;
e6  = e4 * e2;
eps = e2 / (1-e2);
em1 = sqrt(1-e2);
e1  = (1-em1)/(1+em1);
e12 = e1*e1;

lambda_not  = ((-180 + zone*6) - 3) * DEG2RADS;

% remove the false things

x = x - (RADIUS>1)*FALSE_EAST;
if (hemisphere)
  y = y - (RADIUS>1)*FALSE_NORTH;
end

% compute the footpoint latitude

M = y/K_NOT;
mu = M/(RADIUS * (1 - 0.25*e2 - 0.046875*e4 - 0.01953125*e6));
foot = mu + (1.5*e1 - 0.84375*e12*e1)*sin(2*mu) ...
    + (1.3125*e12 - 1.71875*e12*e12)*sin(4*mu) ...
    + (1.57291666667*e12*e1)*sin(6*mu) ...
    + (2.142578125*e12*e12)*sin(8*mu);

% some other terms

sinF = sin(foot);
cosF = cos(foot);
tanF = tan(foot);

N = RADIUS ./ sqrt(1-e2*(sinF.*sinF));
T = tanF.*tanF;
T2 = T.*T;
C = eps * cosF.*cosF;
C2 = C.*C;
denom = sqrt(1-e2*(sinF.*sinF));
R = RADIUS * em1*em1 ./ (denom.*denom.*denom);
D = x./(N*K_NOT);
D2 = D.*D;
D4 = D2.*D2;

% can now compute the lat and lon

lat = foot - (N.*tanF./R) .* (0.5*D2 - (5 + 3*T + 10*C - 4*C2 - 9*eps).*D4/24 ...
    + (61 + 90*T + 298*C + 45*T2 - 252*eps - 3*C2) .* D4 .* D2/720);

lon = lambda_not + (D - (1 + 2*T +C).*D2.*D/6 + ...
    (5 - 2*C + 28*T - 3*C2 + 8*eps + 24*T2).*D4.*D./120)./cosF;


% convert back to degrees;

lat=lat/DEG2RADS;
lon=lon/DEG2RADS;

return