-
Notifications
You must be signed in to change notification settings - Fork 3
/
polarstereo_inv.m
149 lines (133 loc) · 5.56 KB
/
polarstereo_inv.m
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
function [phi,lambda]=polarstereo_inv(x,y);
a = 6378137;
e = 0.08181919;
phi_c = 70;
lambda_0 = -45;
%POLARSTEREO_INV transforms map coordinates to lat/lon data for a polar stereographic system
% [LAT,LON]=POLARSTEREO_INV(X,Y,EARTHRADIUS,ECCENTRICITY,LAT_TRUE,LON_POSY)
% LAT and LON are in decimal degrees with negative numbers (-) for S and W.
% X and Y are the map coordinates (scalars, vectors, or matrices of equal size).
% EARTHRADIUS is the radius of the earth defined in the projection
% (default is 6378137.0 m, WGS84)
% ECCENTRICITY is the earth's misshapenness
% (default is 0.08181919)
% LAT_TRUE is the latitude of true scale in degrees, aka standard parallel
% (default is -70). Note that some NSIDC data use 70 and some use 71.
% LON_POSY is the meridian in degrees along the positive Y axis of the map
% (default is 0)
%
% The National Snow and Ice Data Center (NSIDC) and Scientific Committee
% on Antarctic Research (SCAR) use a version of the polar stereographic
% projection that Matlab does not have. This file does transformations to
% geographic coordinates from map coordinates to facilitate
% comparisons with other datasets.
%
% Equations from: Map Projections - A Working manual - by J.P. Snyder. 1987
% http://kartoweb.itc.nl/geometrics/Publications/Map%20Projections%20-%20A%20Working%20manual%20-%20by%20J.P.%20Snyder.pdf
% See the section on Polar Stereographic, with a south polar aspect and
% known phi_c not at the pole.
%
% See also: MINVTRAN, MFWDTRAN.
%
% Written by Andy Bliss, 9/12/2011
% Changes since version 01:
% 1. Split into two functions and vectorized code.
%%%%%%%%%%%%%
%some standard info
%%%%%%%%%%%%%
%WGS84 - radius: 6378137.0 eccentricity: 0.08181919
% in Matlab: axes2ecc(6378137.0, 6356752.3142)
%Hughes ellipsoid - radius: 6378.273 km eccentricity: 0.081816153
% Used for SSM/I http://nsidc.org/data/polar_stereo/ps_grids.html
%International ellipsoid (following Snyder) - radius: 6378388.0 eccentricity: 0.0819919
%{
%check the code using Snyder's example. Should be about -75 (S) and 150 (E).
x=-1540033.6;
y=-560526.4;
[phi,lambda]=polarstereo_inv(x,y,6378388.0,0.0819919,-71,-100);
phi,lambda
%%%%%%%%%%%%
%check with AntDEM
%%%%%%%%%%%%
%http://nsidc.org/data/docs/daac/nsidc0304_0305_glas_dems.gd.html
% Center Point of Corner Grid Cell
%x y Latitude Longitude
test=[-2812000.0 2299500.0 -57.3452815 -50.7255753
2863500.0 2299500.0 -57.0043684 51.2342036
-2812000.0 -2384000.0 -56.8847122 -130.2911169
2863500.0 -2384000.0 -56.5495152 129.7789915];
[phi,lambda]=polarstereo_inv(test(:,1),test(:,2),6378137.0,axes2ecc(6378137.0, 6356752.3),-70,0);
figure,hold on,plot(test(:,4),test(:,3),'.'),plot(lambda,phi,'r+')
[test(:,3) test(:,3)-phi],[test(:,4) test(:,4)-lambda]
%error is less than half a meter (probably just round-off error).
%%%%%%%%%%%%
%check with Greenland
%%%%%%%%%%%%
%projected from the WGS 84 Ellipsoid, with 70° N as the latitude of true scale and a rotation of 45.
test=[-890000.0 -629000.0 79.9641229 -99.7495626 %center point of cell
1720000.0 -629000.0 73.2101234 24.9126514
-890000.0 -3410000.0 58.2706251 -59.6277136
1720000.0 -3410000.0 55.7592932 -18.2336765];
[phi,lambda]=polarstereo_inv(test(:,1),test(:,2),6378273,0.081816153,70,-45); %slightly off
[phi2,lambda2]=polarstereo_inv(test(:,1),test(:,2),6378137.0,0.08181919,70,-45); %correct, FINALLY!
figure,hold on,plot(test(:,4),test(:,3),'.'),plot(lambda,phi,'r+'),plot(lambda2,phi2,'gx')
[test(:,3) test(:,3)-phi test(:,3)-phi2],[test(:,4) test(:,4)-lambda test(:,4)-lambda2]
%error is less than half a meter (probably just round-off error).
%}
%%%%%%%%%%%%
%input checking
%%%%%%%%%%%%
% if nargin < 3 || isempty(a)
% a=6378137.0; %radius of ellipsoid, WGS84
% end
% if nargin < 4 || isempty(e)
% e=0.08181919; %eccentricity, WGS84
% end
% if nargin < 5 || isempty(phi_c)
% phi_c=-70; %standard parallel, latitude of true scale
% end
% if nargin < 6 || isempty(lambda_0)
% lambda_0=0; %meridian along positive Y axis
% end
%convert to radians
phi_c=deg2rad(phi_c);
lambda_0=deg2rad(lambda_0);
%if the standard parallel is in S.Hemi., switch signs.
if phi_c < 0
pm=-1; %plus or minus, north lat. or south
phi_c=-phi_c;
lambda_0=-lambda_0;
x=-x;
y=-y;
else
pm=1;
end
%this is not commented very well. See Snyder for details.
t_c=tan(pi/4 - phi_c/2)/((1-e*sin(phi_c))/(1+e*sin(phi_c)))^(e/2);
m_c=cos(phi_c)/sqrt(1-e^2*(sin(phi_c))^2);
rho=sqrt(x.^2+y.^2);
t=rho*t_c/(a*m_c);
%iterate to find phi, with a threshold of pi*1e-8
% phi_alt=pi/2 - 2 * atan(t); %guess for phi
% phiold=phi_alt+10; %+10 to make sure it executes the while loop
% while any(abs(phi_alt(:)-phiold(:)) > pi*1e-8)
% phiold=phi_alt;
% phi_alt=pi/2 - 2*atan(t.*((1-e*sin(phi_alt))./(1+e*sin(phi_alt))).^(e/2));
% %add a break in case it doesn't converge?
% end
%find phi with a series instead of iterating.
chi=pi/2 - 2 * atan(t);
phi=chi+(e^2/2 + 5*e^4/24 + e^6/12 + 13*e^8/360)*sin(2*chi)...
+ (7*e^4/48 + 29*e^6/240 + 811*e^8/11520)*sin(4*chi)...
+ (7*e^6/120+81*e^8/1120)*sin(6*chi)...
+ (4279*e^8/161280)*sin(8*chi);
% phicheck=phi-phi_alt %difference is usually < 1e-10
lambda=lambda_0 + atan2(x,-y);
%correct the signs and phasing
phi=pm*phi;
% phi_alt=pm*phi_alt;
lambda=pm*lambda;
lambda=mod(lambda+pi,2*pi)-pi; %want longitude in the range -pi to pi
%convert back to degrees
phi=rad2deg(phi);
lambda=rad2deg(lambda);