Code covered by the BSD License  

Highlights from
Geodetic Transformations Toolbox

Geodetic Transformations Toolbox

by

 

19 Jan 2006 (Updated )

Set of tools for transformation used in geodesy, especially when using GPS or mapping

ELL=utm2ell(PRO,ell,UND,UseMGRSold,FileOut)
function ELL=utm2ell(PRO,ell,UND,UseMGRSold,FileOut)

% UTM2ELL performs transformation from a utm mapping projection to ellipsoidal coordinates
%
% ELL=utm2ell(PRO,ell,UND,UseMGRSold,FileOut)
%
% Also necessary:   Projections2.mat   Ellipsoids.mat   (see beneath)
%
% Inputs:  PRO  Cell array nx1 with one string per line containing coordinates
%               If for (some) points orthometric height is available, PRO may be a nx2-cell array
%               with doubles as height information
%               PRO may also be a file name with ASCII data to be processed. No point IDs, only UTM
%               coordinates with possible height separated by blanks.
%
%               Possible formats are:  
%                    grid    32U460000.123 5498765.123   arbitrary digits  (string) or
%                            32U 460000.123 5498765.123  arbitrary digits  (string)
%                   mgrsX    32UMV6000098765                    no digits  (string)
%                mgrsXold    32UMF6000098765                    no digits  (string)
%                            X is the number of MGRS integer places (0-5) 
%
%                The appropriate format is detected automatically, except mgrsXold.
%                This will only be used when UseMGRSold is set to a value other than 0.
%                The format in the input PRO must not vary!
%
%          ell  the underlying ellipsoid as string in lower case letters
%               Standard if omitted or set to [] is 'grs80'.
%               See Ellipsoids.m for details.                 
%
%          UND  undulation values from geoid model to calculate projected height from ellipsoidal height.
%               If omitted or set to [], no correction is done on the height in ELL.
%               Ellipsoidal height = Projected height + Undulation value
%
%   UseMGRSold  When MGRS is used in old format, this input parameter has te be set to a value
%               other than 0. Default if omitted or set to [] is 0.
%
%      FileOut  File to write the output to. If omitted, no output file is generated.
%
% Outputs: ELL  nx2- resp. nx3- matrix with longitude, latitude (and height) on the underlying
%               ellipsoid in [degree]
%               Southern hemisphere is signalled by negative latitude.
%
% utm2ell is meant as first step for geodetic transformations, when utm projection coordinates have to
% be transformed either to global coordinates or other projections.
%
% Note: tm2ell may also used with UTM projection parameters instead of utm2ell.
%       The difference is regarding the input data. tm2ell will work on double arrays like they are outputted
%       by ell2tm while utm2ell will only work on cell arrays like outputted by ell2utm.
%       For limitation differences regarding poles and non-standard zones see ell2tm / ell2utm.

% Author:
% Peter Wasmeier, Technical University of Munich
% p.wasmeier@bv.tum.de
% Jan 18, 2006

%% Do some input checking

% Load input file if specified
if ischar(PRO)
    fid=fopen(PRO,'r');
    i=1;
    while 1 % Pre-allocation counter
        line=fgetl(fid);
        if ~ischar(line)
            break
        end
        i=i+1;
    end
    Coo=cell(i,1);
    PRO=cell(i,1);
    i=1;
    while 1
        line=fgetl(fid);
        if ~ischar(line)
            break
        end
        Coo{i,1}=strtrim(line);
        i=i+1;
    end
    fclose(fid);
    for i=1:length(Coo)
        letters=find(and(double(Coo{i})>=65,double(Coo{i})<=91));
        if length(letters)==1
            tpos=2;
        elseif length(letters)==3
            tpos=1;
        else
            error('Format description of input file is not valid (1 or 3 letters needed)!')
        end
        temp=strtrim(Coo{i}(letters(end)+1:end));
        spac=findstr(temp,' ');
        if length(spac)>=tpos,
            PRO{i,1}=[Coo{i}(1:letters(end)) strtrim(temp(1:spac(tpos)))];
            PRO{i,2}=str2num(temp(spac(tpos)+1:end));
        else
            PRO{i,1}=Coo{i};
        end
    end
end

% Check input size
if ~any(ismember(size(PRO),[1 2]))
    error('Coordinate list PRO must be nx1 or nx2 when being a cell array!')
elseif (ismember(size(PRO,1),[1 2]))&&(~ismember(size(PRO,2),[1 2]))
    PRO=PRO';
end

% Check input format and search for letters in UTM strings
letters=and(double(PRO{1,1})>=65,double(PRO{1,1})<=91);
woletters=find(letters);
if sum(letters)==1
    format='grid';
elseif sum(letters)==3
    format='mgrs';
    places=(length(PRO{1,1})-woletters(end))/2;
    format=[format,num2str(places)];
else
    error('Format description of cell array is not valid (1 or 3 letters needed)!')
end

% Allocate variables
nct=size(PRO,1);  % Number of coordinates to transform
if size(PRO,2)==1
    P=zeros(size(PRO,1),2);
else
    P=zeros(size(PRO,1),3);
end
ID=zeros(size(PRO,1),2);
L0=zeros(size(PRO,1),1);

% Default settings
if nargin<5
    FileOut=[];
end
if nargin<4 || isempty(UseMGRSold)
    UseMGRSold=0;
end
if nargin<3 || isempty(UND)
    UND=zeros(nct,1);
elseif (numel(UND)~=nct)||(~isvector(UND))
    error('Parameter ''UND'' must be a vector with the length of PRO!')
else
    UND=UND(:)';
end
if nargin<2 || isempty(ell)
    ell='grs80';
end

%% Load ellipsoids
load Ellipsoids;
if ~exist(ell,'var')
    error(['Ellipsoid ',ell,' is not defined in Ellipsoids.mat - check your definitions!.'])
end
eval(['ell=',ell,';']);

%% Undo the formatting
switch format(1:4)
    case 'grid'
        for i=1:nct
            letters=find(and(double(PRO{i,1})>=65,double(PRO{i,1})<=91));
            if ~any(ismember(letters,[1 2 3]))
                error(['Input value #',num2str(i),'is not in grid format.']);
            end
            ThisID=str2num(PRO{i,1}(1:letters(1)-1));
            if (isempty(ThisID))||(ThisID==0)
                ID(i,1)=0;    % Polar region
                L0=0;
            else
                ID(i,1)=ThisID;
                L0(i)=ID(i,1)*6-183;
            end
            ID(i,2)=double(PRO{i,1}(letters(1)));
            P(i,1:2)=str2num(PRO{i,1}(letters(1)+1:end));
            if (ID(i,1)~=0) && (ID(i,2)<78) % southern hemisphere
                P(i,2)=P(i,2)-1e7;
            end
            if size(PRO,2)==2
                P(i,3)=PRO{i,2};
            end
        end
    case 'mgrs'
        easting_table=[65:72; 74:78 80:82; 83:90; -4e5:1e5:3e5];
        pole_easting_table=[65:67 70:72 74:76 80:85 88:90; 0:1e5:17e5];
        northing_zone=[67:72 74:78 80:88];
        pole_northing_table=[65:72 74:78 80:90;-7e5:1e5:16e5;-12e5:1e5:11e5];
        if (UseMGRSold==0)
            northing_table=[65:72 74:78 80:86];
        else
            northing_table=[76:78 80:86 65:72 74:75];
        end
        for i=1:nct
            letters=find(and(double(PRO{i,1})>=65,double(PRO{i,1})<=90));
            if ~any(ismember(letters(1),[1 2 3]))||length(letters)~=3
                error(['Input value #',num2str(i),'is not in mgrs format.']);
            end
            ThisID=str2num(PRO{i,1}(1:letters(1)-1));
            if (isempty(ThisID))||(ThisID==0)
                ID(i,1)=0;    % Polar region
                L0=0;
                ID(i,2)=double(PRO{i,1}(letters(1)));
                Pshort=[str2num(PRO{i,1}(letters(3)+1:letters(3)+places)) str2num(PRO{i,1}(letters(3)+1+places:end))];
                [woi,woj]=find(pole_easting_table==double(PRO{i,1}(letters(2))));
                if any(double(PRO{i,1}(letters(1)))==[65 89])  % western
                    P(i,1)=2e6+Pshort(1)*10^(5-places)+pole_easting_table(2,woj)-18e5;
                else
                    P(i,1)=2e6+Pshort(1)*10^(5-places)+pole_easting_table(2,woj);
                end
                [woi,woj]=find(pole_northing_table==double(PRO{i,1}(letters(3))));
                if any(double(PRO{i,1}(letters(1)))==[65 66])  % south pole
                    P(i,2)=2e6+Pshort(2)*10^(5-places)+pole_northing_table(3,woj);
                else
                    P(i,2)=2e6+Pshort(2)*10^(5-places)+pole_northing_table(2,woj);
                end
            else
                ID(i,1)=ThisID;
                L0(i)=ID(i,1)*6-183;
                ID(i,2)=double(PRO{i,1}(letters(1)));
                Pshort=[str2num(PRO{i,1}(letters(3)+1:letters(3)+places)) str2num(PRO{i,1}(letters(3)+1+places:end))];
                [woi,woj]=find(easting_table==double(PRO{i,1}(letters(2))));
                P(i,1)=5e5+easting_table(4,woj)+Pshort(1)*10^(5-places);
                approx_north=885000*(find(northing_zone==ID(i,2))-10.5);
                shortened_north=(find(northing_table==double(PRO{i,1}(letters(3))))-1)*1e5;
                if mod(ID(i,1),2)==0    % even zone - starting with 'F'
                    shortened_north=shortened_north-5e5;
                end
                diffn=approx_north-shortened_north;
                multn=round(diffn/2e6)*2e6;
                P(i,2)=multn+shortened_north+Pshort(2)*10^(5-places);
            end
            if size(PRO,2)==2
                P(i,3)=PRO{i,2};
            end
        end
end

%% Do the calculations

ELL=zeros(size(P));

% 1. eccentricity
e2=(ell.a^2-ell.b^2)/ell.a^2;
% 2. eccentricity
es2=(ell.a^2-ell.b^2)/ell.b^2;

% Transverse mercator - all except poles
if any(ID(:,1)>0)
    y=P(ID(:,1)>0,1)-5e5;
    x=P(ID(:,1)>0,2);
    m0=0.9996;

    n=(ell.a-ell.b)/(ell.a+ell.b);
    es2=(ell.a^2-ell.b^2)/ell.b^2;

    B=(1+n)/ell.a/(1+n^2/4+n^4/64)*x/m0;

    b1=3/2*n*(1-9/16*n^2)*sin(2*B);
    b2=n^2/16*(21-55/2*n^2)*sin(4*B);
    b3=151/96*n^3*sin(6*B);
    b4=1097/512*n^4*sin(8*B);

    Bf=B+b1+b2+b3+b4;

    % shortened Longitude:
    Vf=sqrt(1+es2*cos(Bf).^2);
    etaf=sqrt(es2*cos(Bf).^2);

    ys=y*ell.b/m0/ell.a^2;
    l=atan(Vf./cos(Bf).*sinh(ys).*(1-etaf.^4.*ys.^2/6-es2.*ys.^4/10));
    ELL(ID(:,1)>0,1)=l*180/pi-183+ID(ID(:,1)>0,1)*6;

    % Latitude:
    ELL(ID(:,1)>0,2)=atan(tan(Bf).*cos(Vf.*l).*(1-etaf.^2/6.*l.^4))*180/pi;
end

if any(ID(:,1)==0)    % now handle the poles - ups
    m0=0.994;
    C0=2*ell.a/sqrt(1-e2)*((1-sqrt(e2))/(1+sqrt(e2)))^(sqrt(e2)/2);
    As=e2/2+5*e2^2/24+e2^3/12+13*e2^4/360;
    Bs=7*e2^2/48+29*e2^3/240+811*e2^4/11520;
    Cs=7*e2^3/120+81*e2^4/1120;
    Ds=4279*e2^4/161280;
    
    PolePoints=find(ID(:,1)==0);
    for i=1:length(PolePoints)
        y=P(PolePoints(i),1)-2e6;
        x=P(PolePoints(i),2)-2e6;
        if (ID(PolePoints(i),2)>88) % north pole
            L=atan2(y,-x);
        else
            L=atan2(y,x);
        end
        if (y==0)&&(x==0)       % exactly the pole
            B=pi/2;
        else
            if (y==0)
                R=abs(x);
            elseif (x==0)
                R=abs(y);
            else
                R=abs(y/sin(L));
            end
            xi=pi/2-2*atan(R/m0/C0);
            B=xi+As*sin(2*xi)+Bs*sin(4*xi)+Cs*sin(6*xi)+Ds*sin(8*xi);
        end
        if (ID(PolePoints(i),2)<88)
            B=-B;
        end
        ELL(PolePoints(i),1:2)=[L B]*180/pi;
    end
end

% Height:
if (size(P,2)==3)
    ELL(:,3)=P(:,3)+UND;
end

%% Write output to file if specified

if ~isempty(FileOut)
    fid=fopen(FileOut,'w+');
    if (size(ELL,2)==3)
        fprintf(fid,'%12.10f  %12.10f  %12.10f\n',ELL');
    else
        fprintf(fid,'%12.10f  %12.10f\n',ELL');
    end
    fclose(fid);
end

Contact us