Changeset 6341

Timestamp:

10/19/10 11:12:16 (14 years ago)

Author:

Mathieu Morlighem

Message:

Fixed projection routines thanks to Michael

Location:

issm/trunk/src/m/utils/LatLong

Files:

• latlongfromxy.m (deleted)
• ll2xy.m (modified) (1 diff)
• mapll.m (deleted)
• mapxy.m (deleted)
• utm2ll.m (modified) (2 diffs)
• xy2ll.m (added)
• xyscale_south.m (deleted)

issm/trunk/src/m/utils/LatLong/ll2xy.m

@@ -1 @@
-function [x,y]=ll2xy(lat,lon);
-%LL2XY - convert latitude and longitude coordinates to x and y
+function [x,y] = ll2xy(lat,lon,sgn)  
+%LL2XY - converts lat long to polar stereographic
+%
+%   Converts from geodetic latitude and longitude to Polar 
+%   Stereographic (X,Y) coordinates for the polar regions.
+%   Author: Michael P. Schodlok, December 2003 (map2ll)
 %
 %   Usage:
-%      [x,y]=ll2xy(lat,lon)
-%
-%   See also MAPLL, MAPXY
+%      [x,y] = map2ll(lat,lon,sgn)
+%      - sgn = Sign of latitude +1 : north latitude
+%                               -1 : south latitude
 
-lon=lon+360; % to have 0<lon<360
+if nargin ~= 3
+        help map2ll
+        error('bad usage');
+end
 
-re = 6378137.0;         % WGS84
-e2 = 0.00669437999015;  % WGS84
-sn=-1.0;        % because it's southern hemisphere
+% Conversion constant from degrees to radians
+cde  = 57.29577951;
+% Standard latitude for the SSM/I grid with no distorsion
+slat = 70.;
+% Radius of the earth in meters
+re   = 6378.273*10^3;
+% Eccentricity of the Hughes ellipsoid squared
+ex2   = .006693883;
+% Eccentricity of the Hughes ellipsoid
+ex    =  sqrt(ex2);
 
-a=re;
-e=sqrt(e2);
+if sgn == 1 
+        delta = 45.;
+else
+        delta = 0.0;
+end
 
-phi  = sn*lat*pi/180;
-lambda  = lon*pi/180; 
+latitude  = abs(lat) * pi/180.;
+longitude = (lon + delta) * pi/180.;
 
-qp = 1 - (1-e2)/2/e*log((1-e)/(1+e));
-%m=cos(phi)/sqrt(1-e2.*sin(phi)*sin(phi));
-q=(1-e2)*(sin(phi)./(1-e2*sin(phi).^2)-0.5/e*log((1-e*sin(phi))./(1+e*sin(phi))));
-rho = a*sqrt(qp-q);
-x = rho.*sin(lambda);
-y = - sn*rho.*cos(lambda);
+% compute X and Y in grid coordinates.
+T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2);
 
+if (90 - slat) <  1.e-5 
+        rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex));
+else
+        sl  = slat*pi/180.;
+        tc  = tan(pi/4.-sl/2.)/((1.-ex*sin(sl))/(1.+ex*sin(sl)))^(ex/2.);
+        mc  = cos(sl)/sqrt(1.0-ex2*(sin(sl)^2));
+        rho = re*mc*T/tc;
+end
 
-%lon=lon-360; % to recover the initial longitude
+y = -rho .* sgn .* cos(sgn.*longitude);
+x =  rho .* sgn .* sin(sgn.*longitude);
 
+[cnt1,cnt2] = find(latitude >= pi / 2.);
+
+if cnt1
+        x(cnt1,1) = 0.0
+        y(cnt1,1) = 0.0
+end

issm/trunk/src/m/utils/LatLong/utm2ll.m

@@ -50 +50 @@
 end
 
-
-
 % Memory pre-allocation
 %
 Lat=zeros(n1,1);
 Lon=zeros(n1,1);
-
 
 % Main Loop
@@ -119 +116 @@
 
 end
-