Index: /issm/trunk-jpl/src/m/coordsystems/ll2utm.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/ll2utm.m (revision 18460) +++ /issm/trunk-jpl/src/m/coordsystems/ll2utm.m (revision 18460) @@ -0,0 +1,200 @@ +function [x,y,f]=ll2utm(lat,lon,datum) +%LL2UTM Lat/Lon to UTM coordinates precise conversion. +% [X,Y]=LL2UTM2(LAT,LON) or LL2UTM([LAT,LON]) converts coordinates +% LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84. +% +% LAT and LON can be scalars, vectors or matrix. Outputs X and Y will +% have the same size as inputs. +% +% LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be +% a string in the following list: +% 'wgs84': World Geodetic System 1984 (default) +% 'nad27': North American Datum 1927 +% 'clk66': Clarke 1866 +% 'nad83': North American Datum 1983 +% 'grs80': Geodetic Reference System 1980 +% 'int24': International 1924 / Hayford 1909 +% or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in +% meters) and F is flattening of the user-defined ellipsoid. +% +% [X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value +% stands for southern hemisphere points). +% +% +% XY = LL2UTM(...) or without any output argument returns a 2-column +% matrix [X,Y]. +% +% Notice: +% - LL2UTM does not perform cross-datum conversion. +% - precision is near a millimeter. +% +% +% Reference: +% I.G.N., Projection cartographique Mercator Transverse: Algorithmes, +% Notes Techniques NT/G 76, janvier 1995. +% +% Author: Francois Beauducel, +% Created: 2003-12-02 +% Updated: 2014-04-20 + + +% Copyright (c) 2001-2014, François Beauducel, covered by BSD License. +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are +% met: +% +% * Redistributions of source code must retain the above copyright +% notice, this list of conditions and the following disclaimer. +% * Redistributions in binary form must reproduce the above copyright +% notice, this list of conditions and the following disclaimer in +% the documentation and/or other materials provided with the distribution +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +% Available datums +datums = [ ... + { 'wgs84', 6378137.0, 298.257223563 }; + { 'nad83', 6378137.0, 298.257222101 }; + { 'grs80', 6378137.0, 298.257222101 }; + { 'nad27', 6378206.4, 294.978698214 }; + { 'int24', 6378388.0, 297.000000000 }; + { 'clk66', 6378206.4, 294.978698214 }; +]; + +if nargin == 1 + if size(lat,2) ~= 2 + error('Single input argument must be a 2-column matrix [LAT,LON].') + end + lon = lat(:,2); + lat = lat(:,1); +end + +if nargin < 1 + error('Not enough input arguments.') +end + +if all([numel(lat),numel(lon)] > 1) && any(size(lat) ~= size(lon)) + error('LAT and LON must be the same size or scalars.') +end + +if nargin < 3 + datum = 'wgs84'; +end + +if ischar(datum) + if ~any(strcmpi(datum,datums(:,1))) + error('Unkown DATUM name "%s"',datum); + end + k = find(strcmpi(datum,datums(:,1))); + A1 = datums{k,2}; + F1 = datums{k,3}; +else + if numel(datum) ~= 2 + error('User defined DATUM must be a vector [A,F].'); + end + A1 = datum(1); + F1 = datum(2); +end + +% constants +D0 = 180/pi; % conversion rad to deg + +K0 = 0.9996; % UTM scale factor +X0 = 500000; % UTM false East (m) + +B1 = A1*(1 - 1/F1); +E1 = sqrt((A1*A1 - B1*B1)/(A1*A1)); + +p1 = lat/D0; % Phi = Latitude (rad) +l1 = lon/D0; % Lambda = Longitude (rad) +F0 = round((l1*D0 + 183)/6); % UTM zone +P0 = 0/D0; +L0 = (6*F0 - 183)/D0; % UTM origin longitude (rad) +Y0 = 1e7*(p1 < 0); % UTM false northern (m) +N = K0*A1; + +C = coef(E1,0); +B = C(1)*P0 + C(2)*sin(2*P0) + C(3)*sin(4*P0) + C(4)*sin(6*P0) + C(5)*sin(8*P0); +YS = Y0 - N*B; + +C = coef(E1,2); +L = log(tan(pi/4 + p1/2).*(((1 - E1*sin(p1))./(1 + E1*sin(p1))).^(E1/2))); +z = complex(atan(sinh(L)./cos(l1 - L0)),log(tan(pi/4 + asin(sin(l1 - L0)./cosh(L))/2))); +Z = N.*C(1).*z + N.*(C(2)*sin(2*z) + C(3)*sin(4*z) + C(4)*sin(6*z) + C(5)*sin(8*z)); +xs = imag(Z) + X0; +ys = real(Z) + YS; + +% outputs zone if needed: scalar value if unique, or vector/matrix of the +% same size as x/y in case of crossed zones +if nargout > 2 + fu = unique(F0.*sign(lat)); + if isscalar(fu) + f = fu; + else + f = F0; + end +end + +if nargout < 2 + x = [xs(:),ys(:)]; +else + x = xs; + y = ys; +end + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function c = coef(e,m) +%COEF Projection coefficients +% COEF(E,M) returns a vector of 5 coefficients from: +% E = first ellipsoid excentricity +% M = 0 for transverse mercator +% M = 1 for transverse mercator reverse coefficients +% M = 2 for merdian arc + + +if nargin < 2 + m = 0; +end + +switch m + case 0 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + -105/4096, 0, -45/1024, 0, -3/32, 0, -3/8, 0, 0; + 525/16384, 0, 45/1024, 0, 15/256, 0, 0, 0, 0; + -175/12288, 0, -35/3072, 0, 0, 0, 0, 0, 0; + 315/131072, 0, 0, 0, 0, 0, 0, 0, 0]; + + case 1 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + 1/61440, 0, 7/2048, 0, 1/48, 0, 1/8, 0, 0; + 559/368640, 0, 3/1280, 0, 1/768, 0, 0, 0, 0; + 283/430080, 0, 17/30720, 0, 0, 0, 0, 0, 0; + 4397/41287680, 0, 0, 0, 0, 0, 0, 0, 0]; + + case 2 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + -901/184320, 0, -9/1024, 0, -1/96, 0, 1/8, 0, 0; + -311/737280, 0, 17/5120, 0, 13/768, 0, 0, 0, 0; + 899/430080, 0, 61/15360, 0, 0, 0, 0, 0, 0; + 49561/41287680, 0, 0, 0, 0, 0, 0, 0, 0]; + +end +c = zeros(size(c0,1),1); + +for i = 1:size(c0,1) + c(i) = polyval(c0(i,:),e); +end Index: /issm/trunk-jpl/src/m/coordsystems/utm2ll.m =================================================================== --- /issm/trunk-jpl/src/m/coordsystems/utm2ll.m (revision 18459) +++ /issm/trunk-jpl/src/m/coordsystems/utm2ll.m (revision 18460) @@ -1,117 +1,190 @@ -function [Lat,Lon] = utm2ll(xx,yy,utmzone) -% ------------------------------------------------------------------------- -% [Lat,Lon] = utm2ll(x,y,utmzone) +function [lat,lon]=utm2ll(x,y,f,datum,varargin) +%UTM2LL UTM to Lat/Lon coordinates precise conversion. +% [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters) +% defined in the UTM ZONE (integer) to latitude LAT and longitude LON +% (in degrees). Default datum is WGS84. % -% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84). -% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez. +% X and Y can be scalars, vectors or matrix. Outputs LAT and LON will +% have the same size as inputs. % -% Inputs:on) -% -3.485713 7.801235 -119.955246 -17.759537 -94.799019 121.640266 +% For southern hemisphere points, use negative zone -ZONE. % -% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds -% [Lat, Lon]=utm2ll(x,y,utmzone); -% LatDMS=dms2mat(deg2dms(Lat)) -%LatDMS = -% 40.00 18.00 55.55 -% 46.00 17.00 2.01 -% 37.00 34.00 40.17 -% 28.00 38.00 44.33 -% 38.00 51.00 19.96 -% 25.00 3.00 42.41 -% LonDMS=dms2mat(deg2dms(Lon)) -%LonDMS = -% -3.00 29.00 8.61 -% 7.00 48.00 4.40 -% -119.00 57.00 18.93 -% -17.00 45.00 34.33 -% -94.00 47.00 56.47 -% 121.00 38.00 24.96 +% UTM2LL(X,Y,ZONE,DATUM) uses specific DATUM for conversion. DATUM can be +% a string in the following list: +% 'wgs84': World Geodetic System 1984 (default) +% 'nad27': North American Datum 1927 +% 'clk66': Clarke 1866 +% 'nad83': North American Datum 1983 +% 'grs80': Geodetic Reference System 1980 +% 'int24': International 1924 / Hayford 1909 +% or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in +% meters) and F is flattening of the user-defined ellipsoid. % -% Author: -% Rafael Palacios -% Universidad Pontificia Comillas -% Madrid, Spain -% Version: Apr/06, Jun/06, Aug/06 -% Aug/06: corrected m-Lint warnings -%------------------------------------------------------------------------- +% Notice: +% - UTM2LL does not perform cross-datum conversion. +% - precision is near a millimeter. +% +% +% Reference: +% I.G.N., Projection cartographique Mercator Transverse: Algorithmes, +% Notes Techniques NT/G 76, janvier 1995. +% +% Author: Francois Beauducel, +% Created: 2001-08-23 +% Updated: 2014-04-20 -% Argument checking +% Copyright (c) 2001-2014, François Beauducel, covered by BSD License. +% All rights reserved. % -error(nargchk(3, 3, nargin)); %3 arguments required -n1=length(xx); -n2=length(yy); -n3=size(utmzone,1); -if (n1~=n2 || n1~=n3) - error('x,y and utmzone vectors should have the same number or rows'); -end -c=size(utmzone,2); -if (c~=4) - error('utmzone should be a vector of strings like "30 T"'); +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are +% met: +% +% * Redistributions of source code must retain the above copyright +% notice, this list of conditions and the following disclaimer. +% * Redistributions in binary form must reproduce the above copyright +% notice, this list of conditions and the following disclaimer in +% the documentation and/or other materials provided with the distribution +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. + +% Available datums +datums = [ ... + { 'wgs84', 6378137.0, 298.257223563 }; + { 'nad83', 6378137.0, 298.257222101 }; + { 'grs80', 6378137.0, 298.257222101 }; + { 'nad27', 6378206.4, 294.978698214 }; + { 'int24', 6378388.0, 297.000000000 }; + { 'clk66', 6378206.4, 294.978698214 }; +]; + +if nargin < 3 + error('Not enough input arguments.') end -% Memory pre-allocation -% -Lat=zeros(n1,1); -Lon=zeros(n1,1); +if all([numel(x),numel(y)] > 1) && any(size(x) ~= size(y)) + error('X and Y must be the same size or scalars.') +end -% Main Loop -% -for i=1:n1 - if (utmzone(i,4)>'X' || utmzone(i,4)<'C') - fprintf('utm2ll: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n'); +if ~isnumeric(f) || ~isscalar(f) || f ~= round(f) + error('ZONE must be a scalar integer.') +end + +if nargin < 4 + datum = 'wgs84'; +end + +if ischar(datum) + if ~any(strcmpi(datum,datums(:,1))) + error('Unkown DATUM name "%s"',datum); end - if (utmzone(i,4)>'M') - hemis='N'; % Northern hemisphere - else - hemis='S'; + k = find(strcmpi(datum,datums(:,1))); + A1 = datums{k,2}; + F1 = datums{k,3}; +else + if numel(datum) ~= 2 + error('User defined DATUM must be a vector [A,F].'); end + A1 = datum(1); + F1 = datum(2); +end - x=xx(i); - y=yy(i); - zone=str2double(utmzone(i,1:2)); +% constants +D0 = 180/pi; % conversion rad to deg +maxiter = 100; % maximum iteration for latitude computation +eps = 1e-11; % minimum residue for latitude computation - sa = 6378137.000000 ; sb = 6356752.314245; +K0 = 0.9996; % UTM scale factor +X0 = 500000; % UTM false East (m) +Y0 = 1e7*(f < 0); % UTM false North (m) +P0 = 0; % UTM origin latitude (rad) +L0 = (6*abs(f) - 183)/D0; % UTM origin longitude (rad) +E1 = sqrt((A1^2 - (A1*(1 - 1/F1))^2)/A1^2); % ellpsoid excentricity +N = K0*A1; - % e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa; - e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb; - e2cuadrada = e2 ^ 2; - c = ( sa ^ 2 ) / sb; - % alpha = ( sa - sb ) / sa; %f - % ablandamiento = 1 / alpha; % 1/f +% computing parameters for Mercator Transverse projection +C = coef(E1,0); +YS = Y0 - N*(C(1)*P0 + C(2)*sin(2*P0) + C(3)*sin(4*P0) + C(4)*sin(6*P0) + C(5)*sin(8*P0)); - X = x - 500000; +C = coef(E1,1); +zt = complex((y - YS)/N/C(1),(x - X0)/N/C(1)); +z = zt - C(2)*sin(2*zt) - C(3)*sin(4*zt) - C(4)*sin(6*zt) - C(5)*sin(8*zt); +L = real(z); +LS = imag(z); - if hemis == 'S' || hemis == 's' - Y = y - 10000000; - else - Y = y; - end +l = L0 + atan(sinh(LS)./cos(L)); +p = asin(sin(L)./cosh(LS)); - S = ( ( zone * 6 ) - 183 ); - lat = Y / ( 6366197.724 * 0.9996 ); - v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996; - a = X / v; - a1 = sin( 2 * lat ); - a2 = a1 * ( cos(lat) ) ^ 2; - j2 = lat + ( a1 / 2 ); - j4 = ( ( 3 * j2 ) + a2 ) / 4; - j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3; - alfa = ( 3 / 4 ) * e2cuadrada; - beta = ( 5 / 3 ) * alfa ^ 2; - gama = ( 35 / 27 ) * alfa ^ 3; - Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 ); - b = ( Y - Bm ) / v; - Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2; - Eps = a * ( 1 - ( Epsi / 3 ) ); - nab = ( b * ( 1 - Epsi ) ) + lat; - senoheps = ( exp(Eps) - exp(-Eps) ) / 2; - Delt = atan(senoheps / (cos(nab) ) ); - TaO = atan(cos(Delt) * tan(nab)); - longitude = (Delt *(180 / pi ) ) + S; - latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ... - (180 / pi); +L = log(tan(pi/4 + p/2)); - Lat(i)=latitude; - Lon(i)=longitude; +% calculates latitude from the isometric latitude +p = 2*atan(exp(L)) - pi/2; +p0 = NaN; +n = 0; +while any(isnan(p0) | abs(p - p0) > eps) & n < maxiter + p0 = p; + es = E1*sin(p0); + p = 2*atan(((1 + es)./(1 - es)).^(E1/2).*exp(L)) - pi/2; + n = n + 1; +end +if nargout < 2 + lat = D0*[p(:),l(:)]; +else + lat = p*D0; + lon = l*D0; end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function c = coef(e,m) +%COEF Projection coefficients +% COEF(E,M) returns a vector of 5 coefficients from: +% E = first ellipsoid excentricity +% M = 0 for transverse mercator +% M = 1 for transverse mercator reverse coefficients +% M = 2 for merdian arc + + +if nargin < 2 + m = 0; +end + +switch m + case 0 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + -105/4096, 0, -45/1024, 0, -3/32, 0, -3/8, 0, 0; + 525/16384, 0, 45/1024, 0, 15/256, 0, 0, 0, 0; + -175/12288, 0, -35/3072, 0, 0, 0, 0, 0, 0; + 315/131072, 0, 0, 0, 0, 0, 0, 0, 0]; + + case 1 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + 1/61440, 0, 7/2048, 0, 1/48, 0, 1/8, 0, 0; + 559/368640, 0, 3/1280, 0, 1/768, 0, 0, 0, 0; + 283/430080, 0, 17/30720, 0, 0, 0, 0, 0, 0; + 4397/41287680, 0, 0, 0, 0, 0, 0, 0, 0]; + + case 2 + c0 = [-175/16384, 0, -5/256, 0, -3/64, 0, -1/4, 0, 1; + -901/184320, 0, -9/1024, 0, -1/96, 0, 1/8, 0, 0; + -311/737280, 0, 17/5120, 0, 13/768, 0, 0, 0, 0; + 899/430080, 0, 61/15360, 0, 0, 0, 0, 0, 0; + 49561/41287680, 0, 0, 0, 0, 0, 0, 0, 0]; + +end +c = zeros(size(c0,1),1); + +for i = 1:size(c0,1) + c(i) = polyval(c0(i,:),e); +end