Home > atmlab > math > gridinterp.m

gridinterp

PURPOSE ^

GRIDINTERP Change of rectangular grid in 1D to 5D

SYNOPSIS ^

function B = gridinterp(agrids,A,newgrids,varargin)

DESCRIPTION ^

 GRIDINTERP  Change of rectangular grid in 1D to 5D

    Atmlabs function for interpolation that can be seen as a change the
    underlaying grid(s). That is, conversion of data between two
    multi-dimensional rectangular grids. This is basically an interface to
    to the standard interp1 and interpn functions.

    For "point interpolation", see accompanying function *pointinterp*.

    The grids of the data to be interpolated (*A*) are packed into the vector
    array *agrids*. Grids for interpolated data are packed in the same in
    *newgrids*. The length of these grids determine the size of *B*.
    That is, the number of input arguments is not changed by the data
    dimension.

    All grids must be sorted in ascending order (demand of interpn).
    The Atmlab convention of that 1D objects (including grids) are column
    vectors is here strictly followed.

    The function allows a special treatment of extrapolation by the optional
    argument *extrap*. If set to false, standard matlab functionality is
    obtained. If set to true, then data are treated to be constant outside end
    points. That is, the data values at end points are assumed to be valid to
    -Inf andf Inf, respectively. This can also be seen as a "nearest"
    interpolation for positions outside the covered range. This feature is
    also applied on singleton dimensions (resulting in a constant value for
    that dimension between -Inf to Inf).

    A note: For "unsmooth" data, this function should be avoided when going
    from a finer to a coarser grid. For such cases, consider the function
    *resample_geodata*.

 FORMAT   B = gridinterp(agrids,A,newgrids[,iopt])

 OUT   B          Interpolated data.
 IN    agrids     Grids of A, as an array of vectors.
       A          Data to be interpolated.
       newgrids   Grids for interpolation, as an array of vectors.
 OPT   iopt       Interpolation option. See *interpn*. Default is 'linear'.
       extrap     Special treatment of extrapolation. See above.
                  Default is false.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

gridinterp.m

SOURCE CODE ^

0001 % GRIDINTERP  Change of rectangular grid in 1D to 5D
0002 %
0003 %    Atmlabs function for interpolation that can be seen as a change the
0004 %    underlaying grid(s). That is, conversion of data between two
0005 %    multi-dimensional rectangular grids. This is basically an interface to
0006 %    to the standard interp1 and interpn functions.
0007 %
0008 %    For "point interpolation", see accompanying function *pointinterp*.
0009 %
0010 %    The grids of the data to be interpolated (*A*) are packed into the vector
0011 %    array *agrids*. Grids for interpolated data are packed in the same in
0012 %    *newgrids*. The length of these grids determine the size of *B*.
0013 %    That is, the number of input arguments is not changed by the data
0014 %    dimension.
0015 %
0016 %    All grids must be sorted in ascending order (demand of interpn).
0017 %    The Atmlab convention of that 1D objects (including grids) are column
0018 %    vectors is here strictly followed.
0019 %
0020 %    The function allows a special treatment of extrapolation by the optional
0021 %    argument *extrap*. If set to false, standard matlab functionality is
0022 %    obtained. If set to true, then data are treated to be constant outside end
0023 %    points. That is, the data values at end points are assumed to be valid to
0024 %    -Inf andf Inf, respectively. This can also be seen as a "nearest"
0025 %    interpolation for positions outside the covered range. This feature is
0026 %    also applied on singleton dimensions (resulting in a constant value for
0027 %    that dimension between -Inf to Inf).
0028 %
0029 %    A note: For "unsmooth" data, this function should be avoided when going
0030 %    from a finer to a coarser grid. For such cases, consider the function
0031 %    *resample_geodata*.
0032 %
0033 % FORMAT   B = gridinterp(agrids,A,newgrids[,iopt])
0034 %
0035 % OUT   B          Interpolated data.
0036 % IN    agrids     Grids of A, as an array of vectors.
0037 %       A          Data to be interpolated.
0038 %       newgrids   Grids for interpolation, as an array of vectors.
0039 % OPT   iopt       Interpolation option. See *interpn*. Default is 'linear'.
0040 %       extrap     Special treatment of extrapolation. See above.
0041 %                  Default is false.
0042 %
0043 
0044 % 2006-08-22   Created by Patrick Eriksson.
0045 
0046 
0047 function B = gridinterp(agrids,A,newgrids,varargin)
0048 %
0049 [iopt,extrap] = optargs( varargin, { 'linear', false } );
0050 errid= ['atmlab', mfilename, 'badInput'];
0051 if atmlab('STRICT_ASSERT')
0052     %&%
0053     %&%
0054     %- Check input                                                            %&%
0055     %                                                                         %&%
0056     rqre_nargin( 3, nargin )                                                  %&%
0057     rqre_datatype( agrids, @iscell );                                         %&%
0058     rqre_datatype( A, @isnumeric );                                           %&%
0059     rqre_datatype( newgrids, @iscell );                                       %&%
0060     indim  = length(agrids);
0061     outdim = length(newgrids);
0062     if dimens(A) > length(agrids)                                             %&%
0063         error(errid, 'Dimensionality of A is higher than number of given grids.' );   %&%
0064     end                                                                       %&%
0065     if outdim < indim                                                         %&%
0066         error(errid,'There can not be fewer grids in *newgrids* than in *agrids*.');  %&%
0067     end                                                                       %&%
0068     for d = 1 : indim                                                         %&%
0069         rqre_datatype( agrids{d}, @istensor1, 'Grids in *agrids*' );            %&%
0070         rqre_gridmatch( A, d, agrids{d}, true );                                %&%
0071         rqre_datatype( newgrids{d}, @istensor1, 'Grids in *newgrids*' );        %&%
0072     end                                                                       %&%
0073     rqre_datatype( extrap, @isboolean );                                      %&%
0074     
0075     
0076     %- Check a or determine a's effective dimensions
0077     %
0078     if extrap    % Ignore singleton dimensions
0079         dims = find( size( A ) > 1 );
0080         A    = getdims( A, dims );
0081     else
0082         
0083         if outdim ~= indim                                                      %&%
0084             error( errid, ['With *extrap* false, the number of grids in *newgrids* ',... %&%
0085                 'and *agrids* must be the same.'] );  %&%
0086         end                                                                     %&%
0087         asize = size( A );                                                      %&%
0088         if length(asize) < outdim  ||  any( asize(1:outdim)==1 )                %&%
0089             error( errid,['With *extrap* false there can not be any singleton ',...     %&%
0090                 'dimmensions in *A*.'] );   %&%
0091         end                                                                     %&%
0092         dims = 1:outdim;
0093     end
0094     
0095     dim = length( dims );
0096     
0097     if dim == 0
0098         %
0099         B = A;
0100         
0101     elseif dim == 1
0102         %
0103         if extrap
0104             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0105         else
0106             xi = newgrids{dims(1)};
0107         end
0108         %                                                                     %&%
0109         if isempty(xi)                                                        %&%
0110             error(errid,'Empty grid found in *newgrid* for non-singleton dimension.');%&%
0111         end                                                                   %&%
0112         %
0113         B = interp1( agrids{dims(1)}, A, xi, iopt );
0114         
0115     elseif dim == 2
0116         %
0117         if extrap
0118             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0119             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0120         else
0121             xi = newgrids{dims(1)};
0122             yi = newgrids{dims(2)};
0123         end
0124         %                                                                     %&%
0125         if isempty(xi)||isempty(yi)                                            %&%
0126             error( errid,'Empty grid found in *newgrid* for non-singleton dimension.');%&%
0127         end                                                                   %&%
0128         %
0129         [xi,yi] = ndgrid( xi, yi );
0130         %
0131         B = interpn( agrids{dims(1)}, agrids{dims(2)}, A, xi, yi, iopt );
0132         
0133     elseif dim == 3
0134         %
0135         if extrap
0136             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0137             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0138             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0139         else
0140             xi = newgrids{dims(1)};
0141             yi = newgrids{dims(2)};
0142             zi = newgrids{dims(3)};
0143         end
0144         %                                                                     %&%
0145         if isempty(xi)||isempty(yi)||isempty(zi)                                %&%
0146             error(errid,'Empty grid found in *newgrid* for non-singleton dimension.');%&%
0147         end                                                                   %&%
0148         %
0149         [xi,yi,zi] = ndgrid( xi, yi, zi );
0150         %
0151         B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0152             A, xi, yi, zi, iopt );
0153         
0154     elseif dim == 4
0155         %
0156         if extrap
0157             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0158             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0159             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0160             ui = handle_expand( agrids{dims(4)}, newgrids{dims(4)} );
0161         else
0162             xi = newgrids{dims(1)};
0163             yi = newgrids{dims(2)};
0164             zi = newgrids{dims(3)};
0165             ui = newgrids{dims(4)};
0166         end
0167         %                                                                     %&%
0168         if isempty(xi)||isempty(yi)||isempty(zi)||isempty(ui)                    %&%
0169             error(errid,'Empty grid found in *newgrid* for non-singleton dimension.');%&%
0170         end                                                                   %&%
0171         %
0172         [xi,yi,zi,ui] = ndgrid( xi, yi, zi, ui );
0173         %
0174         B = interpn( agrids{dims(1)}, agrids{dims(2)}, ...
0175             agrids{dims(3)}, agrids{dims(4)}, A, xi, yi, zi, ui, iopt );
0176         
0177     elseif dim == 5
0178         %
0179         if extrap
0180             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0181             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0182             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0183             ui = handle_expand( agrids{dims(4)}, newgrids{dims(4)} );
0184             vi = handle_expand( agrids{dims(5)}, newgrids{dims(5)} );
0185         else
0186             xi = newgrids{dims(1)};
0187             yi = newgrids{dims(2)};
0188             zi = newgrids{dims(3)};
0189             ui = newgrids{dims(4)};
0190             vi = newgrids{dims(5)};
0191         end
0192         %                                                                     %&%
0193         if isempty(xi)||isempty(yi)||isempty(zi)||isempty(ui)||isempty(vi)        %&%
0194             error(errid,'Empty grid found in *newgrid* for non-singleton dimension.');%&%
0195         end                                                                   %&%
0196         %
0197         [xi,yi,zi,ui,vi] = ndgrid( xi, yi, zi, ui, vi );
0198         %
0199         B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0200             agrids{dims(4)}, agrids{dims(5)}, A, xi, yi, zi, ui, vi, iopt );
0201         
0202     else                                                             %&%
0203         error( errid,'Interpolation above 5D is not handled.' );             %&%
0204     end
0205     
0206     
0207     %- Post-processing for *extrap*
0208     %
0209     if extrap  &&  dim < outdim
0210         % Make extrapolation for A's singleton dimensions
0211         map = ones( 1, max([2 outdim]) ); % Repmat demands at least two values
0212         lg  = ones( 1, max([2 outdim]) );
0213         for d = 1 : outdim
0214             lg(d) = max([ 1 length( newgrids{d} ) ] );  % The 1 to handle empty grids
0215             if ~any( dims == d )  &&  lg(d) > 1          % for singleton dimensions
0216                 map(d) = length( newgrids{d} );
0217             end
0218         end
0219         if any( map > 1 )
0220             B = repmat( B, map );
0221         end
0222         
0223         % Full dimension is not obtained above for some cases of singleton
0224         % dimensions. Fixed by a reshape
0225         bsize = size(B);
0226         if length(bsize) < length(lg)  ||  any( bsize ~= lg )
0227             B = reshape( B, lg );
0228         end
0229     end
0230 else
0231 
0232     % ------------
0233     % The else contains the same as in the if environment except but without
0234     % assertions.
0235     
0236     outdim = length(newgrids);
0237     %- Check a or determine a's effective dimensions
0238     %
0239     if extrap    % Ignore singleton dimensions
0240         dims = find( size( A ) > 1 );
0241         A    = getdims( A, dims );
0242     else
0243         dims = 1:outdim;
0244     end
0245     
0246     dim = length( dims );
0247     
0248     if dim == 0
0249         B = A;
0250     elseif dim == 1
0251         if extrap
0252             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0253         else
0254             xi = newgrids{dims(1)};
0255         end
0256         B = interp1( agrids{dims(1)}, A, xi, iopt );
0257         
0258     elseif dim == 2
0259         if extrap
0260             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0261             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0262         else
0263             xi = newgrids{dims(1)};
0264             yi = newgrids{dims(2)};
0265         end
0266         [xi,yi] = ndgrid( xi, yi );
0267         B = interpn( agrids{dims(1)}, agrids{dims(2)}, A, xi, yi, iopt );
0268         
0269     elseif dim == 3
0270         
0271         if extrap
0272             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0273             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0274             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0275         else
0276             xi = newgrids{dims(1)};
0277             yi = newgrids{dims(2)};
0278             zi = newgrids{dims(3)};
0279         end
0280         [xi,yi,zi] = ndgrid( xi, yi, zi );
0281         B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0282             A, xi, yi, zi, iopt );
0283         
0284     elseif dim == 4
0285         if extrap
0286             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0287             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0288             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0289             ui = handle_expand( agrids{dims(4)}, newgrids{dims(4)} );
0290         else
0291             xi = newgrids{dims(1)};
0292             yi = newgrids{dims(2)};
0293             zi = newgrids{dims(3)};
0294             ui = newgrids{dims(4)};
0295         end
0296         [xi,yi,zi,ui] = ndgrid( xi, yi, zi, ui );
0297         B = interpn( agrids{dims(1)}, agrids{dims(2)}, ...
0298             agrids{dims(3)}, agrids{dims(4)}, A, xi, yi, zi, ui, iopt );
0299         
0300     elseif dim == 5
0301         if extrap
0302             xi = handle_expand( agrids{dims(1)}, newgrids{dims(1)} );
0303             yi = handle_expand( agrids{dims(2)}, newgrids{dims(2)} );
0304             zi = handle_expand( agrids{dims(3)}, newgrids{dims(3)} );
0305             ui = handle_expand( agrids{dims(4)}, newgrids{dims(4)} );
0306             vi = handle_expand( agrids{dims(5)}, newgrids{dims(5)} );
0307         else
0308             xi = newgrids{dims(1)};
0309             yi = newgrids{dims(2)};
0310             zi = newgrids{dims(3)};
0311             ui = newgrids{dims(4)};
0312             vi = newgrids{dims(5)};
0313         end
0314         [xi,yi,zi,ui,vi] = ndgrid( xi, yi, zi, ui, vi );
0315         B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0316             agrids{dims(4)}, agrids{dims(5)}, A, xi, yi, zi, ui, vi, iopt );
0317         
0318         %- Post-processing for *extrap*
0319         %
0320         if extrap  &&  dim < outdim
0321             % Make extrapolation for A's singleton dimensions
0322             map = ones( 1, max([2 outdim]) ); % Repmat demands at least two values
0323             lg  = ones( 1, max([2 outdim]) );
0324             for d = 1 : outdim
0325                 lg(d) = max([ 1 length( newgrids{d} ) ] );  % The 1 to handle empty grids
0326                 if ~any( dims == d )  &&  lg(d) > 1          % for singleton dimensions
0327                     map(d) = length( newgrids{d} );
0328                 end
0329             end
0330             if any( map > 1 )
0331                 B = repmat( B, map );
0332             end
0333             
0334             % Full dimension is not obtained above for some cases of singleton
0335             % dimensions. Fixed by a reshape
0336             bsize = size(B);
0337             if length(bsize) < length(lg)  ||  any( bsize ~= lg )
0338                 B = reshape( B, lg );
0339             end
0340         end
0341         
0342     end
0343     
0344 end
0345 
0346 
0347 %---
0348 function xi = handle_expand(x,xi)
0349 %
0350 
0351 % use find instead of indexing because the code relies on i1 sometimes being = []
0352 
0353 v1 = min( x );
0354 i1 =  find( xi<v1 ) ;
0355 %
0356 v2 = max( x );
0357 i2 = find( xi>v2 );
0358 %
0359 if ~isempty(v1)
0360     xi(i1) = v1;
0361 end
0362 %
0363 if ~isempty(v2)
0364     xi(i2) = v2;
0365 end
0366 %
0367 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005