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

ntv2trafo.m
function [OUT,UsedGridName]=ntv2trafo(IN,gridlat,gridlong,header,direction,FileOut)

% NTV2TRAFO performs a NTV2 grid transformation
%
% [OUT, UsedGridName] = ntv2trafo(IN, gridlat, gridlong, header, FileOut)
%
% Inputs:        IN  nx2 - matrix of points to be transformed [Long Lat] []
%                    2xn-matrices are allowed, but be careful with 2x2-matrices.
%                    IN may also be a file name with ASCII data to be processed. No point IDs, only
%                    coordinates as if it was a matrix.
%
%           gridlat  cell array of grids with NTv2 shift values (lat)  ["]
%          gridlong  cell array of grid with NTv2 shift values (long) ["]
%            header  cell array of structures containing grid header information
%                    NTV2TRAFO works on inverted longitude sign convention as produced by READNTV2
%                    i.e. western longitudes have negative sign (standardized NTv2 transformation
%                    grids show positive sign for western longitudes).
%                    Please use READNTV2 to get the proper inputs for NTV2TRAFO.
%                    With only one grid entered, the input variables may also be double/struct
%                    matrices instead of cell arrays.
%
%         direction  0 if transformation is in standard direction. (Default if omitted or set to [])
%                    1 if reverse transformation is to be performed.
%                    Per definition, shift values are defined in the source system. So they are meant
%                    to transform Sys1 -> Sys2, but not in reverse direction. If reverse direction is
%                    required, this will be done by this function using an iterative algorithm. See
%                    comments below.
%                   
%           FileOut  File to write the output OUT to. If omitted, no output file is generated.
%
% Outputs:      OUT  nx2 - matrix with the transformed output variables [Long Lat] []
%
%      UsedGridName  nx1 - cell array with the name of the sub_grid which was used for each
%                    transformed point. The names are taken from SUB_NAME field in the header.
%                    If the point is out of specified grids, '* NOT SPECIFIED *' is returned.
%
% For each point it is individually searched for the proper grid. If more than one grid is convenient,
% the algorithm takes the one with the lowest size of mesh (mesh diagonal).
% If a point is outside of the specified grids, a warning is thrown and the point is returned
% untransformed. Use 'warning off NTv2:NoGridAvailable' to suppress the warnings.
%
% Note: NTv2 is an gridwise affine transformation with linear interpolation between knot points.
%       This inherently means that there will be interpolation unsteadiness (a kink) at the knot
%       positions. As a result, different coordinates very close to a knot point, which lie in
%       different meshes may be mapped to identical target coordinates (depending on the mesh width
%       of the used subgrid).
%       As converse argument, when performing a backward transformation ([Sys1a->]Sys2->Sys1b), it may 
%       occur that other coordinates are found for Sys1b than the starting ones of Sys1a in such cases.
%       The error may be up to 0.1" which would lead to about 3 meters.

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

%% Do input checking, set defaults and adjust vectors

if nargin<6
    FileOut=[];
end
if nargin<5 || isempty(direction)
    direction = 0;
end
if nargin<4
    error('Not all parameters have been specified.');
end
if iscell(gridlat)
    if  any(diff([length(gridlat) length(gridlong) length(header)]))
        error('Input cell arrays must be of equal length.')
    end
    for i=1:length(gridlat)
        if any(size(gridlat{i})~=size(gridlong{i}))
            error('GridLat and GridLong dimensions do not match.');
        end
    end
else
    if any(size(gridlat)~=size(gridlong))
        error('GridLat and GridLong dimensions do not match.');
    end
    gridlat={gridlat};gridlong={gridlong};header={header};
end

% Load input file if specified
if ischar(IN) % Load file
    IN=load(IN);
end

if     (size(IN,1)~=2)&&(size(IN,2)~=2)
    error('Coordinate list IN must be a nx2-matrix!')
elseif (size(IN,1)==2)&&(size(IN,2)~=2)
    IN=IN';
end

IN=IN*3600;
OUT=zeros(size(IN));
UsedGridName=cell(size(IN,1),1);

if (direction)
    IN0=IN;
end
dIN=1;

%% Calculate the transformation
while (any(abs(dIN)>1e-10))
    % Find the correct grid to work in and do the calculations
    for k=1:size(IN,1);
        g_inc=inf;
        UsedGrid=0;
        for i=1:length(gridlat)
            if (header{i}.W_LONG <= IN(k,1))&&(IN(k,1) <= header{i}.E_LONG) &&...
                    (header{i}.S_LAT <= IN(k,2))&&(IN(k,2) <= header{i}.N_LAT)
                width=header{i}.LAT_INC^2+header{i}.LONG_INC^2;
                if (width)<g_inc
                    g_inc=width;
                    UsedGrid=i;
                end
            end
        end
        if UsedGrid==0
            warning('NTv2:NoGridAvailable',['Transformation point ',num2str(k),' is out of the specified grids.']);
            OUT(k,:)=IN(k,:);
            UsedGridName{k}='* NOT SPECIFIED *';
            continue;
        end
        UsedGridName{k}=header{UsedGrid}.SUB_NAME;
        fcol=(IN(k,1)-header{UsedGrid}.W_LONG)/header{UsedGrid}.LONG_INC;
        frow=(IN(k,2)-header{UsedGrid}.S_LAT)/header{UsedGrid}.LAT_INC;
        col=floor(fcol);
        row=floor(frow);
        if col==size(gridlat{UsedGrid},2)-1
            col=col-1;
        end
        if row==size(gridlat{UsedGrid},1)-1
            row=row-1;
        end
        dx=fcol-col;
        dy=frow-row;
        svlat=(1-dx)*(1-dy)*gridlat{UsedGrid}(size(gridlat{UsedGrid},1)-row,col+1)+...
            dx *(1-dy)*gridlat{UsedGrid}(size(gridlat{UsedGrid},1)-row,col+2)+...
            dx *   dy *gridlat{UsedGrid}(size(gridlat{UsedGrid},1)-row-1,col+2)+...
            (1-dx)*   dy *gridlat{UsedGrid}(size(gridlat{UsedGrid},1)-row-1,col+1);
        svlong=(1-dx)*(1-dy)*gridlong{UsedGrid}(size(gridlat{UsedGrid},1)-row,col+1)+...
            dx *(1-dy)*gridlong{UsedGrid}(size(gridlat{UsedGrid},1)-row,col+2)+...
            dx *   dy *gridlong{UsedGrid}(size(gridlat{UsedGrid},1)-row-1,col+2)+...
            (1-dx)*   dy *gridlong{UsedGrid}(size(gridlat{UsedGrid},1)-row-1,col+1);
        OUT(k,:)=(IN(k,:)+[svlong svlat]);
    end
    if (~direction)
        break
    end
    dIN=OUT-IN0;
    IN=IN-dIN;
end
if (direction)
    OUT=IN;
end
OUT=OUT/3600;

%% Write output to file if specified

if ~isempty(FileOut)
    fid=fopen(FileOut,'w+');
    fprintf(fid,'%12.10f  %12.10f\n',OUT');
    fclose(fid);
end

Contact us