Changeset 18478


Ignore:
Timestamp:
09/02/14 12:05:17 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: new version after my suggestion to the author

Location:
issm/trunk-jpl/src/m/coordsystems
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r18460 r18478  
    1 function [x,y,f]=ll2utm(lat,lon,datum)
     1function [x,y,f]=ll2utm(varargin)
    22%LL2UTM Lat/Lon to UTM coordinates precise conversion.
    3 %   [X,Y]=LL2UTM2(LAT,LON) or LL2UTM([LAT,LON]) converts coordinates
     3%       [X,Y]=LL2UTM2(LAT,LON) or LL2UTM([LAT,LON]) converts coordinates
    44%       LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84.
    55%
    … …  
    77%       have the same size as inputs.
    88%
    9 %       LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be
    10 %       a string in the following list:
     9%       LL2UTM(...,DATUM) uses specific DATUM for conversion. DATUM can be one
     10%       of the following char strings:
    1111%               'wgs84': World Geodetic System 1984 (default)
    1212%               'nad27': North American Datum 1927
    … …  
    1818%       meters) and F is flattening of the user-defined ellipsoid.
    1919%
    20 %       [X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value
    21 %       stands for southern hemisphere points).
    22 %
    23 %
    24 %       XY = LL2UTM(...) or without any output argument returns a 2-column
     20%       LL2UTM(...,ZONE) forces the UTM ZONE (scalar integer) instead of
     21%       automatic set.
     22%
     23%       [X,Y,ZONE]=LL2UTM(...) returns also the computed UTM ZONE (negative
     24%       value for southern hemisphere points).
     25%
     26%
     27%       XY=LL2UTM(...) or without any output argument returns a 2-column
    2528%       matrix [X,Y].
    2629%
    27 %       Notice:
     30%       Note:
    2831%               - LL2UTM does not perform cross-datum conversion.
    2932%               - precision is near a millimeter.
    … …  
    3437%                  Notes Techniques NT/G 76, janvier 1995.
    3538%
    36 %   Author: Francois Beauducel, <beauducel@ipgp.fr>
    37 %   Created: 2003-12-02
    38 %   Updated: 2014-04-20
     39%       Acknowledgments: Mathieu.
     40%
     41%
     42%       Author: Francois Beauducel, <beauducel@ipgp.fr>
     43%       Created: 2003-12-02
     44%       Updated: 2014-08-24
    3945
    4046
    … …  
    7480];
    7581
    76 if nargin == 1
    77         if size(lat,2) ~= 2
    78                 error('Single input argument must be a 2-column matrix [LAT,LON].')
    79         end
    80         lon = lat(:,2);
    81         lat = lat(:,1);
    82 end
    83        
     82% constants
     83D0 = 180/pi;    % conversion rad to deg
     84K0 = 0.9996;    % UTM scale factor
     85X0 = 500000;    % UTM false East (m)
     86
     87% defaults
     88datum = 'wgs84';
     89zone = [];
     90
    8491if nargin < 1
    8592        error('Not enough input arguments.')
    8693end
    8794
     95if isnumeric(varargin{1}) & size(varargin{1},2) == 2
     96        lat = varargin{1}(:,1);
     97        lon = varargin{1}(:,2);
     98        v = 1;
     99elseif nargin > 1 & isnumeric(varargin{2})
     100        lat = varargin{1};
     101        lon = varargin{2};
     102        v = 2;
     103else
     104        error('Single input argument must be a 2-column matrix [LAT,LON].')
     105end
     106
    88107if all([numel(lat),numel(lon)] > 1) && any(size(lat) ~= size(lon))
    89108        error('LAT and LON must be the same size or scalars.')
    90109end
    91110
    92 if nargin < 3
    93         datum = 'wgs84';
     111for n = (v+1):nargin
     112        % LL2UTM(...,DATUM)
     113        if ischar(varargin{n}) | (isnumeric(varargin{n}) & numel(varargin{n})==2)
     114                datum = varargin{n};
     115        % LL2UTM(...,ZONE)
     116        elseif isnumeric(varargin{n}) & isscalar(varargin{n})
     117                        zone = round(varargin{n});
     118        else
     119                error('Unknown argument #%d. See documentation.',n)
     120        end
    94121end
    95122
    96123if ischar(datum)
     124        % LL2UTM(...,DATUM) with DATUM as char
    97125        if ~any(strcmpi(datum,datums(:,1)))
    98126                error('Unkown DATUM name "%s"',datum);
    … …  
    102130        F1 = datums{k,3};       
    103131else
    104         if numel(datum) ~= 2
    105                 error('User defined DATUM must be a vector [A,F].');
    106         end
     132        % LL2UTM(...,DATUM) with DATUM as [A,F] user-defined
    107133        A1 = datum(1);
    108134        F1 = datum(2);
    109135end
    110136
    111 % constants
    112 D0 = 180/pi;    % conversion rad to deg
    113 
    114 K0 = 0.9996;    % UTM scale factor
    115 X0 = 500000;    % UTM false East (m)
     137p1 = lat/D0;                    % Phi = Latitude (rad)
     138l1 = lon/D0;                    % Lambda = Longitude (rad)
     139
     140% UTM zone automatic setting
     141if isempty(zone)
     142        F0 = round((l1*D0 + 183)/6);
     143else
     144        F0 = zone;
     145end
    116146
    117147B1 = A1*(1 - 1/F1);
    118148E1 = sqrt((A1*A1 - B1*B1)/(A1*A1));
    119 
    120 p1 = lat/D0;                                    % Phi = Latitude (rad)
    121 l1 = lon/D0;                                    % Lambda = Longitude (rad)
    122 F0 = round((l1*D0 + 183)/6);    % UTM zone
    123149P0 = 0/D0;
    124 L0 = (6*F0 - 183)/D0;                   % UTM origin longitude (rad)
    125 Y0 = 1e7*(p1 < 0);                              % UTM false northern (m)
     150L0 = (6*F0 - 183)/D0;           % UTM origin longitude (rad)
     151Y0 = 1e7*(p1 < 0);              % UTM false northern (m)
    126152N = K0*A1;
    127153

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

    r18460 r18478  
    11function [lat,lon]=utm2ll(x,y,f,datum,varargin)
    22%UTM2LL UTM to Lat/Lon coordinates precise conversion.
    3 %   [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters)
     3%       [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters)
    44%       defined in the UTM ZONE (integer) to latitude LAT and longitude LON
    55%       (in degrees). Default datum is WGS84.
    … …  
    3030%                  Notes Techniques NT/G 76, janvier 1995.
    3131%
    32 %   Author: Francois Beauducel, <beauducel@ipgp.fr>
    33 %   Created: 2001-08-23
    34 %   Updated: 2014-04-20
     32%       Author: Francois Beauducel, <beauducel@ipgp.fr>
     33%       Created: 2001-08-23
     34%       Updated: 2014-04-20
     35
    3536
    3637%       Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
    … …  
    105106eps = 1e-11;    % minimum residue for latitude computation
    106107
    107 K0 = 0.9996;                                                            % UTM scale factor
    108 X0 = 500000;                                                            % UTM false East (m)
    109 Y0 = 1e7*(f < 0);                                                       % UTM false North (m)
    110 P0 = 0;                                                                         % UTM origin latitude (rad)
    111 L0 = (6*abs(f) - 183)/D0;                                       % UTM origin longitude (rad)
     108K0 = 0.9996;                                    % UTM scale factor
     109X0 = 500000;                                    % UTM false East (m)
     110Y0 = 1e7*(f < 0);                               % UTM false North (m)
     111P0 = 0;                                         % UTM origin latitude (rad)
     112L0 = (6*abs(f) - 183)/D0;                       % UTM origin longitude (rad)
    112113E1 = sqrt((A1^2 - (A1*(1 - 1/F1))^2)/A1^2);     % ellpsoid excentricity
    113114N = K0*A1;
Note: See TracChangeset for help on using the changeset viewer.