Changeset 27627

Timestamp:

03/01/23 12:03:43 (2 years ago)

Author:

Mathieu Morlighem

Message:

CHG: latest version

Location:

issm/trunk-jpl/src/m/coordsystems

Files:

• ll2utm.m (modified) (7 diffs)
• utm2ll.m (modified) (6 diffs)

issm/trunk-jpl/src/m/coordsystems/ll2utm.m ¶

@@ -18 +18 @@
 %       meters) and F is flattening of the user-defined ellipsoid.
 %
-%       LL2UTM(...,ZONE) forces the UTM ZONE (scalar integer) instead of
-%       automatic set.
+%       LL2UTM(...,ZONE) forces the UTM ZONE (scalar integer or same size as
+%   LAT and LON) instead of automatic set.
 %
 %       [X,Y,ZONE]=LL2UTM(...) returns also the computed UTM ZONE (negative
@@ -37 +37 @@
 %                  Notes Techniques NT/G 76, janvier 1995.
 %
-%       Acknowledgments: Mathieu.
+%       Acknowledgments: Mathieu, Frederic Christen.
 %
 %
 %       Author: Francois Beauducel, <beauducel@ipgp.fr>
 %       Created: 2003-12-02
-%       Updated: 2014-08-24
-
-
-%       Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%       Updated: 2019-05-29
+
+
+%       Copyright (c) 2001-2019, François Beauducel, covered by BSD License.
 %       All rights reserved.
 %
@@ -93 +93 @@
 end
 
-if isnumeric(varargin{1}) & size(varargin{1},2) == 2
+if nargin > 1 && isnumeric(varargin{1}) && isnumeric(varargin{2}) ...
+                && (all(size(varargin{1})==size(varargin{2})) ...
+                || isscalar(varargin(1)) || isscalar(varargin{2}))
+        lat = varargin{1};
+        lon = varargin{2};
+        v = 2;
+elseif isnumeric(varargin{1}) && size(varargin{1},2) == 2
         lat = varargin{1}(:,1);
         lon = varargin{1}(:,2);
         v = 1;
-elseif nargin > 1 & isnumeric(varargin{2})
-        lat = varargin{1};
-        lon = varargin{2};
-        v = 2;
 else
         error('Single input argument must be a 2-column matrix [LAT,LON].')
@@ -109 +111 @@
 end
 
+if any(abs(lat)>90)
+        error('LAT absolute values must be lower than 90.')
+end
+
+% checks for DATUM and/or ZONE syntax
+% NOTE: the following strategy works in any case except if ZONE argument
+% has a size of 1x2 (in that case it will be interpreted as a DATUM). To
+% force the ZONE syntax with 2 elements, just use ZONE(:) to make a colum
+% vector of 2x1.
 for n = (v+1):nargin
         % LL2UTM(...,DATUM)
-        if ischar(varargin{n}) | (isnumeric(varargin{n}) & numel(varargin{n})==2)
+        if ischar(varargin{n}) || (isnumeric(varargin{n}) ...
+                        && all(size(varargin{n})==[1,2]))
                 datum = varargin{n};
         % LL2UTM(...,ZONE)
-        elseif isnumeric(varargin{n}) & isscalar(varargin{n})
-                        zone = round(varargin{n});
+        elseif isnumeric(varargin{n}) && (isscalar(varargin{n}) ...
+                        || (isscalar(lat) || all(size(varargin{n})==size(lat))) ...
+                        && (isscalar(lon) || all(size(varargin{n})==size(lon))))
+                zone = round(varargin{n});
         else
                 error('Unknown argument #%d. See documentation.',n)
@@ -142 +156 @@
         F0 = round((l1*D0 + 183)/6);
 else
-        F0 = zone;
+        F0 = abs(zone);
 end
 
@@ -148 +162 @@
 E1 = sqrt((A1*A1 - B1*B1)/(A1*A1));
 P0 = 0/D0;
-L0 = (6*F0 - 183)/D0;           % UTM origin longitude (rad)
+L0 = (6*F0 - 183)/D0;   % UTM origin longitude (rad)
 Y0 = 1e7*(p1 < 0);              % UTM false northern (m)
 N = K0*A1;
@@ -166 +180 @@
 % same size as x/y in case of crossed zones
 if nargout > 2
-        fu = unique(F0.*sign(lat));
+        f = F0.*sign(lat);
+        fu = unique(f);
         if isscalar(fu)
                 f = fu;
-        else
-                f = F0;
         end
 end

issm/trunk-jpl/src/m/coordsystems/utm2ll.m ¶

@@ -5 +5 @@
 %       (in degrees). Default datum is WGS84.
 %
-%       X and Y can be scalars, vectors or matrix. Outputs LAT and LON will
+%       X, Y and F can be scalars, vectors or matrix. Outputs LAT and LON will
 %       have the same size as inputs.
 %
@@ -32 +32 @@
 %       Author: Francois Beauducel, <beauducel@ipgp.fr>
 %       Created: 2001-08-23
-%       Updated: 2014-04-20
-
-
-%       Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%       Updated: 2019-05-29
+
+%       Revision history:
+%
+%       [2019-05-29]
+%         - fix an issue when X or Y are matrices.
+
+%       Copyright (c) 2001-2019, François Beauducel, covered by BSD License.
 %       All rights reserved.
 %
@@ -74 +78 @@
 end
 
-if all([numel(x),numel(y)] > 1) && any(size(x) ~= size(y))
-        error('X and Y must be the same size or scalars.')
-end
-
-if ~isnumeric(f) || ~isscalar(f) || f ~= round(f)
-        error('ZONE must be a scalar integer.')
+% checks if input arguments have compatible sizes using unique(complex)
+sz = [size(x);size(y);size(f)];
+sz = complex(sz(:,1),sz(:,2));
+if length(unique(sz(sz~=complex(1,1)))) > 1
+        error('X, Y and ZONE must be scalar or vector/matrix of the same size.')
+end
+
+if ~isnumeric(f) || any(f ~= round(f))
+        error('ZONE must be integer value.')
 end
 
@@ -101 +108 @@
 end
 
+% calculations are made on column vectors
+x = x(:);
+y = y(:);
+f = f(:);
+
 % constants
 D0 = 180/pi;    % conversion rad to deg
@@ -133 +145 @@
 p0 = NaN;
 n = 0;
-while any(isnan(p0) | abs(p - p0) > eps) & n < maxiter
+while any(isnan(p0) | abs(p - p0) > eps) && n < maxiter
         p0 = p;
         es = E1*sin(p0);
@@ -141 +153 @@
 
 if nargout < 2
-        lat = D0*[p(:),l(:)];
+        lat = D0*[p,l];
 else
-        lat = p*D0;
-        lon = l*D0;
+        % reshapes vectors to x/y/f original size
+        sz = max([size(x);size(y);size(f)]);
+        lat = reshape(p*D0,sz);
+        lon = reshape(l*D0,sz);
 end