Home > atmlab > geographical > resample_geodata.m

resample_geodata

PURPOSE ^

RESAMPLE_GEODATA regrid the data to a new gridsize or fitting to new lat lons

SYNOPSIS ^

function out = resample_geodata(in,field,newgrid)

DESCRIPTION ^

 RESAMPLE_GEODATA regrid the data to a new gridsize or fitting to new lat lons

 PURPOSE:  To resample GRIDDED data according a given gridsize or lat/lon
           vectors. The data is regridded using matlabs INTERPOLATION
           functions if the grid resolution should be FINER and the data
           is RESAMPLED using area wieghting averages using a local plane
           assumption if the grid resolution should be COARSER.



 IN:

   in          struct      contains: lat,lon,field1,etc (maybe more fields)
                           opt: interpolation = 'str'; default='linear'
                                for interpolation method in interp2.

   field       %s          {'field1',etc}; The field/s  you want to regrid

   newgrid     %f, or [%f %f], {[newlat],[newlon]} (e.g. newlat = -90:5:90)

              a) Either make a sqaure grid using a scalar %f. e.g 1deg ->
                 2deg (Assumes centered data!)

              b) or a 2 element vector [%f %f] means [latstep, lonstep]
                 e.g. 1*1 -> 1.25*2.75 (Assumes centered data!)

              c) or {[newlat],[newlon]} if you want to match to a predefined
              lat lon grid. (Use this if your data is NOT centered)

 OUT:         struct       Original structure 'in' with lat, lon, and specified fields
                           replaced with their regridded versions


 NOTE:     - Interpolation is linear, unless specified by in.interpolation.
           - With a combination of conditions
                 1) regridding to a very coarse resolution
                 2) the finer grid resolution doesn't cleanly fit in the larger grid.
             the local plane assumption will cause some errors, the remedy is
             to include area weighting depending on latitude. This is 
             yet to be implemented. 
           - The input data and the newgrid (if full vectors) must be
           equidegree (homogeneous spacing)
           - If you use new grid = [latstep, lonstep], both have to be higher or lower
           than the resolution of the input data, otherwise the
           interpolation function is applied also to the coarser grid, but
           this will throw a warning.

         Important:  (but does not apply is you provide new lat lons in regrid{} )
                   - It is assumed that the grid values are centered!!
                   - For the interpolation/resampling it is always assumed
                   that the whole world is in the data.


 USAGE:         out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},3);
         e.g.   out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},[3,5]);
                out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},...
                        {[ newlats ],[ newlons ]} );

                1st example: Someone wants to resample 7 data fields
                in one go to a 3x3 deg grid (the common lat lons are in the structure)
                2nd example: Same as above but to a 3x5 deg grid
                3rd example: Same as above but match the data to certain
                             lat lons

 IMPORTANT: data fields must have the dimensions data(lat,lon,.....), or
                                                 data(lon,lat,.....)

 Created by Salomon Eliasson
 $Id: regrid_dataset.m 6417 2011-04-21 06:44:37Z seliasson $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

resample_geodata.m

SOURCE CODE ^

0001 function out = resample_geodata(in,field,newgrid)
0002 % RESAMPLE_GEODATA regrid the data to a new gridsize or fitting to new lat lons
0003 %
0004 % PURPOSE:  To resample GRIDDED data according a given gridsize or lat/lon
0005 %           vectors. The data is regridded using matlabs INTERPOLATION
0006 %           functions if the grid resolution should be FINER and the data
0007 %           is RESAMPLED using area wieghting averages using a local plane
0008 %           assumption if the grid resolution should be COARSER.
0009 %
0010 %
0011 %
0012 % IN:
0013 %
0014 %   in          struct      contains: lat,lon,field1,etc (maybe more fields)
0015 %                           opt: interpolation = 'str'; default='linear'
0016 %                                for interpolation method in interp2.
0017 %
0018 %   field       %s          {'field1',etc}; The field/s  you want to regrid
0019 %
0020 %   newgrid     %f, or [%f %f], {[newlat],[newlon]} (e.g. newlat = -90:5:90)
0021 %
0022 %              a) Either make a sqaure grid using a scalar %f. e.g 1deg ->
0023 %                 2deg (Assumes centered data!)
0024 %
0025 %              b) or a 2 element vector [%f %f] means [latstep, lonstep]
0026 %                 e.g. 1*1 -> 1.25*2.75 (Assumes centered data!)
0027 %
0028 %              c) or {[newlat],[newlon]} if you want to match to a predefined
0029 %              lat lon grid. (Use this if your data is NOT centered)
0030 %
0031 % OUT:         struct       Original structure 'in' with lat, lon, and specified fields
0032 %                           replaced with their regridded versions
0033 %
0034 %
0035 % NOTE:     - Interpolation is linear, unless specified by in.interpolation.
0036 %           - With a combination of conditions
0037 %                 1) regridding to a very coarse resolution
0038 %                 2) the finer grid resolution doesn't cleanly fit in the larger grid.
0039 %             the local plane assumption will cause some errors, the remedy is
0040 %             to include area weighting depending on latitude. This is
0041 %             yet to be implemented.
0042 %           - The input data and the newgrid (if full vectors) must be
0043 %           equidegree (homogeneous spacing)
0044 %           - If you use new grid = [latstep, lonstep], both have to be higher or lower
0045 %           than the resolution of the input data, otherwise the
0046 %           interpolation function is applied also to the coarser grid, but
0047 %           this will throw a warning.
0048 %
0049 %         Important:  (but does not apply is you provide new lat lons in regrid{} )
0050 %                   - It is assumed that the grid values are centered!!
0051 %                   - For the interpolation/resampling it is always assumed
0052 %                   that the whole world is in the data.
0053 %
0054 %
0055 % USAGE:         out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},3);
0056 %         e.g.   out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},[3,5]);
0057 %                out = resample_geodata(in,{'CC','HCC','MCC','LCC','P','T','Q'},...
0058 %                        {[ newlats ],[ newlons ]} );
0059 %
0060 %                1st example: Someone wants to resample 7 data fields
0061 %                in one go to a 3x3 deg grid (the common lat lons are in the structure)
0062 %                2nd example: Same as above but to a 3x5 deg grid
0063 %                3rd example: Same as above but match the data to certain
0064 %                             lat lons
0065 %
0066 % IMPORTANT: data fields must have the dimensions data(lat,lon,.....), or
0067 %                                                 data(lon,lat,.....)
0068 %
0069 % Created by Salomon Eliasson
0070 % $Id: regrid_dataset.m 6417 2011-04-21 06:44:37Z seliasson $
0071 
0072 errID = ['atmlab:' mfilename ':badInput'];
0073 assert(nargin==3,errID,'Incorrect number of input arguments')
0074 if ~iscell(field), field={field}; end
0075 
0076 % rename the structure so nothing is lost after this function
0077 out = in;
0078 
0079 % CHECK input
0080 sz                                  = check_input(in,field,newgrid,errID);
0081 
0082 % ASSEMBLE_newgrid
0083 [newlat,newlon]                     = assemble_newgrid(in.lat,in.lon,newgrid,errID);
0084 
0085 % STANDARDIZE "new grid" internally
0086 [rg_flags,newlat,newlon]            = standardize_geodata(newlat,newlon);
0087 
0088 % LEAVE function if nothing needs to be done
0089 if isequal(in.lat(:),newlat(:)) && isequal(in.lon(:),newlon(:))
0090     warning(errID,['input grid is identical to requested new grid.\n'...
0091         'Leaving data untouched'])
0092     out = in; return
0093 end
0094 
0095 % CHECK WHICH test
0096 out.method                          = which_test(in.lat,in.lon,newlat,newlon);
0097 
0098 % PRINT to screen what is about to happen
0099 if strcmp(out.method,'shift grid')
0100     logtext(1,'shifting the grid, keeping the same resolution\n')
0101 else
0102     logtext(1,'Changing grid from %.3fx%.3f to %.3fx%.3f deg by method: %s\n',...
0103         mean(abs(diff(in.lat))),mean(abs(diff(in.lon))),...
0104         mean(diff(newlat)),mean(diff(newlon)),....
0105         out.method)
0106 end
0107 
0108 if ~strcmp(out.method,'interpolate') && isfield(in,'interpolation')
0109     warning(errID,'Will not interpolate, but instead use method: %s',out.method)
0110 elseif strcmp(out.method,'interpolate')
0111     if isfield(in,'interpolation')
0112         intm = in.interpolation;
0113     else
0114         intm = 'linear';
0115     end
0116 end
0117 
0118 %% LOOP over FIELDS to regrid
0119 for F = field
0120     fsz = size(in.(F{1}));
0121 
0122     % COLLAPSE data onto the 3rd dimesion
0123     data = in.(F{1})(:,:,:);
0124     
0125     % STANDARDIZE "input" data internally
0126     [flags,lat,lon,data]  = standardize_geodata(in.lat,in.lon,data);
0127 
0128     % CHECK for every field
0129     assert(isequal(fsz(1:2),sz(1:2)),...
0130         errID,'All input fields for processing don''t have matching lat lon dims')
0131    
0132     
0133     % USE appropriate METHOD
0134     switch out.method
0135         case 'interpolate'
0136             % successfully tested for 4D data and lon=0:360, lat = end:-1:1
0137             data = interpolate_grid(struct('lat',lat,'lon',lon,'data',data,'method',intm),newlat,newlon);
0138         case 'coarse grid'
0139             % successfully tested for 4D data and lon=0:360, lat = end:-1:1
0140             data = coarser_grid(struct('lat',lat,'lon',lon,'data',data),newlat,newlon);
0141         case 'coarse grid with area weighting'
0142             % successfully tested for 4D data and lon=0:360, lat = end:-1:1
0143             data = coarser_grid_withAreaweighting(struct('lat',lat,'lon',lon,'data',data),newlat,newlon);
0144         case 'shift grid'
0145             % successfully: regular, 4D,
0146             data = interpolate_grid(struct('lat',lat,'lon',lon,'data',data,'method',intm),newlat,newlon);
0147         case 'interp and averaging'
0148             data = fine_and_coarse(struct('lat',lat,'lon',lon,'data',data),newlat,newlon,true);
0149         case 'averaging and iterp'
0150             data = fine_and_coarse(struct('lat',lat,'lon',lon,'data',data),newlat,newlon,false);
0151         otherwise
0152             error('atmlab:resample_geodata:Bug','Shouldn''t be here')
0153     end
0154     
0155     % REORDER data to original sorting order
0156     if flags.lat_descend
0157         data = data(end:-1:1,:,:);
0158     end
0159     if flags.lon_descend
0160         data = data(:,end:-1:1,:);
0161     end
0162     if rg_flags.lon360 % the picked newgrid decides...
0163         [~,lnindex] = sort(newlon+(newlon<0)*360);
0164         data = data(:,lnindex,:);
0165     end
0166     if flags.permute
0167         data = permute(data,[2,1,3]);
0168     end
0169 
0170     % RESHAPE data back to original format
0171     out.(F{1}) = reshape(data,[size(data,1),size(data,2),fsz(3:end)]);
0172 end
0173 
0174 out.lat = newlat;
0175 out.lon = newlon;
0176 
0177 % REORDER lat lons if need be
0178 if rg_flags.lat_descend
0179     out.lat = out.lat(end:-1:1);
0180 end
0181 if rg_flags.lon_descend
0182     out.lon = out.lon(end:-1:1);
0183 end
0184 if rg_flags.lon360
0185     out.lon = sort(out.lon+(out.lon<0)*360);
0186 end
0187 
0188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0189 % SUNFUNCTIONS bellow here
0190 %
0191 function sz = check_input(in,field,newgrid,errID)
0192 %% CHECK_input assorted list of input checks
0193 
0194 assert(isstruct(in) && ischar(field{1}) && ...
0195     (isnumeric(newgrid) || iscell(newgrid)),...
0196     errID,'Incorrect input argument types')
0197 
0198 if iscell(newgrid)
0199     assert(~any([isempty(newgrid{1}),isempty(newgrid{2})]),...
0200         errID,'Lat or Lon vector is empty')
0201     
0202 end
0203 
0204 assert(~any(~isfield(in,{'lat','lon'})),errID,'"lat,lon" are required fields')
0205 assert(all(ismember(field,fieldnames(in))),errID,'One or more of the specified data fields is missing');
0206 assert(~isequal([size(in.(field{1}),1),size(in.(field{1}),2)],size(in.lat),size(in.lon)),...
0207     errID,'function only works on regular grids, i.e. data(lat,lon,...) or data(lon,lat,...)')
0208 assert(max(diff(in.lat))-min(diff(in.lat))<0.1&&max(diff(in.lon))-min(diff(in.lon))<0.1,...
0209     errID,'The input data must be equidegree (constant spacing)')
0210 assert( (isequal(in.lat,sort(in.lat)) || isequal(in.lat,sort(in.lat,'descend')) ) && ...
0211     (isequal(in.lon,sort(in.lon)) || isequal(in.lon,sort(in.lon,'descend')) ),...
0212     errID,'both lat & lon must be ascending or descending')
0213 sz = size(in.(field{1}));
0214 assert(isequal([length(in.lat),length(in.lon)],[sz(1),sz(2)]) || ...
0215     isequal([length(in.lat),length(in.lon)],[sz(2),sz(1)]),...
0216     errID, 'data must have the dimensions data(lat,lon,.....), or data(lon,lat,.....)');
0217 
0218 if ~iscell(newgrid)
0219     assert(min(in.lat)>-90 & max(in.lat)<90 ,errID,'Latitudes are not centers')
0220     if min(in.lon)<0
0221         assert(~(all(ismember([-180,180],in.lon))),errID,'Longitudes are probably not centers')
0222     else
0223         assert(~(all(ismember([0,360],in.lon))),errID,'Longitudes are probably not centers')
0224     end
0225 else
0226    assert(length(newgrid)==2,...
0227         errID,'if newgrid is a cell, it must contatin 2 vectors {[LATs] , [LONs]}')
0228     
0229 end
0230 
0231 
0232 function [newlat,newlon] = assemble_newgrid(LAT,LON,newgrid,errID)
0233 %% ASSEMBLE_newgrid
0234 % in:    LAT             orig latvector
0235 %        LON             orig lonvector
0236 %        newgrid         numeric or cell
0237 %
0238 % out:  newlat,newlon vectors to map to, and logical descn flagging is the
0239 %       are in ascending order or not
0240 
0241 if iscell(newgrid)
0242     % E.g. for mapping one dataset to anothers lat/lons
0243     newlat = newgrid{1};
0244     newlon = newgrid{2};
0245     assert(max(diff(newlat))-min(diff(newlat))<0.1&&max(diff(newlon))-min(diff(newlon))<0.1,...
0246         errID,'The input newgrid must be equidegree (constant spacing)')
0247     
0248 else
0249     % make newlat newlon based on range given by in.lat,in.lon (centered)
0250     g=newgrid; if isscalar(g), g = [g,g]; end
0251     dn = mean(abs(diff(LON))); dt = mean(abs(diff(LAT)));
0252     
0253     % LATITUDES
0254     if min(LAT) == LAT(1)
0255         newlat = LAT(1)-dt/2+g(1)/2 : g(1) : LAT(end)+dt/2-g(1)/2;
0256         flt = newlat(end)~=LAT(end)+dt/2-g(1)/2;
0257     else
0258         newlat = LAT(1)+dt/2-g(1)/2 : -g(1) : LAT(end)-dt/2+g(1)/2;
0259         flt = newlat(end)~=LAT(end)-dt/2+g(1)/2;
0260     end
0261     
0262     % LONGITUDES
0263     if min(LON) == LON(1)
0264         newlon = LON(1)-dn/2+g(2)/2 : g(2) : LON(end)+dn/2-g(2)/2;
0265         fln = newlon(end)~=LON(end)+dn/2-g(2)/2;
0266     else
0267         newlon = LON(1)+dn/2-g(2)/2 : -g(2) : LON(end)-dn/2+g(2)/2;
0268         fln = newlon(end)~=LON(end)-dn/2+g(2)/2;
0269     end
0270     
0271     % do they FIT ?
0272     if flt
0273         warning(errID,['Requested latitude range does not fit into original latitude range\n',...
0274             'Suggest to give the latitude vector direcly in newgrid = {[lats][lons]}'])
0275     end
0276     if fln
0277         warning(errID,['Requested longitude range does not fit into original longitude range\n',...
0278             'Suggest to give the longitude vector direcly in newgrid = {[lats][lons]}'])
0279     end
0280 end
0281 
0282 function method = which_test(LAT,LON,newlat,newlon)
0283 %% CHECK_WHICH_test_to_use
0284 % IN:  lat lons of both grids to decide which method to use
0285 % OUT: string name of the method to use
0286 
0287 % first test is to see if the edges of the input grid match the edges of
0288 % the new grid
0289 
0290 eps = 1e-5; %this corresponds to 1 meter at the equator
0291 
0292 % diff averages
0293 dnlt = mean(diff(newlat));
0294 dnln = mean(diff(newlon));
0295 dlt  = mean(diff(LAT));
0296 dln  = mean(diff(LON));
0297 
0298 % temporay vectors
0299 nltb   = newlat - dnlt/2;
0300 nltt   = newlat + dnlt/2;
0301 nlnb   = newlon - dnln/2;
0302 nlnt   = newlon + dnln/2;
0303 ltb    = LAT    - dlt/2;
0304 ltt    = LAT    + dlt/2;
0305 lnb    = LAT    - dln/2;
0306 lnt    = LON    + dln/2;
0307 
0308 % LOGICALS
0309 
0310 %to a finer grid
0311 cond1   = dlt > dnlt; %length(LAT) < length(newlat);
0312 cond2   = dln > dnln; %length(LON) < length(newlon);
0313 
0314 %to a coarser grid
0315 cond3   = dlt < dnlt; %length(LAT) > length(newlat);
0316 cond4   = dln < dnln; %length(LON) > length(newlon);
0317 
0318 % same resolution
0319 cond5   = abs(dlt-dnlt)<eps; %length(LAT) == length(newlat);
0320 cond6   = abs(dln-dnln)<eps; %length(LON) == length(newlon);
0321 
0322 % exact fit
0323 onLatEdges = ismember(nltb(1),ltb) && ismember(nltt(end),ltt);
0324 onLonEdges = ismember(nlnb(1),lnb) && ismember(nlnt(end),lnt);
0325 fitlt      = dnlt/dlt == floor(dnlt/dlt);         %dnlt = n*dlt; n interger
0326 fitln      = dnln/dln == floor(dnln/dln);         %dnln = n*dln; n interger
0327 
0328 if (cond1 && ~cond4) || (cond2 && ~cond3)
0329     % to a FINER grid
0330     
0331     method = 'interpolate';
0332 elseif (cond3 && ~cond2) || (cond4 && ~cond1)
0333     % to a COARSER grid
0334 
0335     if fitlt && fitln && onLatEdges && onLonEdges
0336     
0337         % fit together nicely, this is much faster
0338         method = 'coarse grid';
0339     else
0340         % overlapping gridboxes
0341         method = 'coarse grid with area weighting';
0342     end
0343 elseif cond5 && cond6 
0344     % SAME resolution but DIFFERENT grid
0345     method = 'shift grid';
0346 elseif cond1 && cond4
0347     % to a COMBINATION of both methods (INTERP & AVERAGING)
0348     method = 'interp and averaging';
0349 elseif cond2 && cond3
0350     % to a COMBINATION of both methods (AVERAGING & INTERP)
0351     method = 'averaging and iterp';
0352 else
0353     error('atmlab:resample_geodata:Bug','Shouldn''t be here')
0354 end
0355 
0356 function DATA = interpolate_grid(in,newlat,newlon)
0357 %% INTERPOLATE_GRID interpolate data to a smaller grid.
0358 
0359 % EXPAND the new lat lons so they encompass orig latlon
0360 data = in.data;
0361 
0362 LON=in.lon(:)';
0363 LAT=in.lat(:)';
0364 dn = mean(diff(LON));
0365 dt = mean(diff(LAT));
0366 
0367 while LON(1)>=newlon(1)
0368     LON = [LON(1)-dn,LON,LON(end)+dn];
0369     data = [data(:,end,:),data,data(:,1,:)];
0370 end
0371 while LAT(1)>=newlat(1)
0372     LAT = [LAT(1)-dt,LAT,LAT(end)+dt];
0373     data = [data(1,:,:);data;data(end,:,:)];
0374 end
0375 
0376 % do the INTERPOLATION
0377 [X,Y] = meshgrid(LON,LAT);
0378 [XI,YI] = meshgrid(newlon,newlat);
0379 if islogical(data)
0380     DATA = false([length(newlat),length(newlon),size(data,3)]);
0381 else
0382     DATA = zeros([length(newlat),length(newlon),size(data,3)],class(data));
0383 end
0384 
0385 fprintf('Using interpolation method: %s\n',in.method)
0386 for i = 1:size(data,3)
0387     DATA(:,:,i) = interp2(X,Y,data(:,:,i),XI,YI,in.method);
0388 end
0389 
0390 function data = coarser_grid(in,lat,lon)
0391 %% coarser_grid by averaging.
0392 %
0393 % Purpose: To make gridded data more sparse if the grids fit nicely into each other.
0394 % NOTE: Biases may be introduced because I am using nanmean. The solution
0395 %       is to also output the number of values used in the averaging for
0396 %       each grid (not doing that right now. Add feature if needed)
0397 
0398 % PUT lat/lons LAST and make lat lons both size(lat,lon)
0399 tmpdata = permute(in.data,[3,1,2]); 
0400 tmplon  = repmat(in.lon,length(in.lat),1);
0401 tmplat  = repmat(in.lat,1,length(in.lon));
0402 
0403 % COLLAPSE data (transposed), lat, lon to simulate UNGRIDDED data
0404 tmpdata = binning_fast(struct('data',tmpdata(:,:)','lon',tmplon(:),'lat',tmplat(:),...
0405     'newlon',lon,'newlat',lat));
0406 
0407 % TAKE the grid MEAN for each level (e.g. {[4x288]} -> {[1x288]}).
0408 if any(isnan(tmpdata(:))), warning(['atmlab:' mfilename ':calculation'],...
0409         'There are NaNs in the data, and these may lead to biases on the coarser grid'); end
0410 tmpdata = binned_statistics(tmpdata,{@nanmean});
0411 
0412 if ndims(in.data)==2
0413     data = cell2mat(tmpdata.mean);
0414 elseif ndims(in.data)==3
0415     % make into MATRIX(lat,lat,...)
0416     % this loop is because cell2mat doesn't work if cells contain vectors. It
0417     % doesn't work in conjuction with repmat, it just scrambled
0418     data = zeros(length(lat),length(lon),size(in.data,3));
0419     for i = 1:length(lat)
0420         for j = 1:length(lon)
0421             data(i,j,:) = tmpdata.mean{i,j};
0422         end
0423     end
0424 end
0425 
0426 function data = coarser_grid_withAreaweighting(in,lat,lon)
0427 %% COARSER_GRID make a coarse grid from a finer grid
0428 % method: Uses area weighting for grids that don't fit nicely into each
0429 % other and assumes local plane assumption, which is good as long as the
0430 % grids are not too large (see help resample_data).
0431 % NOTE: Biases may be introduced because I am using nanmean. The solution
0432 %       is to also output the number of values used in the averaging for
0433 %       each grid (not doing that right now. Add feature if needed)
0434 
0435 
0436 % preallocate
0437 if ~islogical(in.data)
0438     data = zeros([length(lat),length(lon),size(in.data,3)],class(in.data));
0439 else data = false([length(lat),length(lon),size(in.data,3)]);
0440 end
0441 
0442 polyWeightSum = @(d,w,i)(nansum(d(i).*w(i))./nansum(w(i)));
0443 
0444 % DETERMINE new grid BOUNDARIES from grid spacing
0445 gn = mean(diff(lon)); gt = mean(diff(lat));
0446 dn = mean(diff(in.lon)); dt = mean(diff(in.lat));
0447 if isnan(gt), gt = (in.lat(end)-in.lat(1) + dt);end %incase there's only one latitude
0448 if isnan(gn), gn = (in.lon(end)-in.lon(1) + dn);end %incase there's only one longitude
0449 
0450 lattop = lat+gt/2+dt;
0451 latbot = lat-gt/2;
0452 lontop = lon+gn/2+dn;
0453 lonbot = lon-gn/2;
0454 
0455 %LOOP over larger grid
0456 for ii = 1:length(lat)
0457     for jj=1:length(lon)
0458         
0459         % GET lats lons in or tangent to grid(i,j)
0460         ilon = in.lon >= lonbot(jj) & in.lon < lontop(jj);
0461         ilat = in.lat >= latbot(ii) & in.lat < lattop(ii);
0462         
0463         % find WEIGHTS
0464         % A = [x,y,width,height] x,y = bottom corner
0465         A = [lonbot(jj),latbot(ii),gn,gt];
0466         [x,y] = meshgrid(in.lon(ilon)-dn/2,in.lat(ilat)-dt/2);
0467         B = [ x(:), y(:), repmat(dn,length(y(:)),1), repmat(dt,length(x(:)),1)];
0468         weights = rectint(A,B);
0469 
0470         % get data SUBSET
0471         d = in.data(ilat,ilon,:);
0472         
0473         % COLLAPSE data on the lat/lon dimensions to match output from rectint
0474         d = permute(d,[3,1,2]);
0475         
0476         for kk = 1:size(d,1)
0477             data(ii,jj,kk)    = polyWeightSum(d(kk,:),weights,~isnan(d(kk,:)));
0478         end
0479     end
0480 end
0481 
0482 function data = fine_and_coarse(in,newlat,newlon,method)
0483 
0484 data = interpolate_grid(in,newlat,newlon);
0485 warning('atmlab:resample_geodata:incomplete',...
0486     ['Using interpolation (interp2) for grids that are finer in one dimension '...
0487     'and coarser in the other. This is NO GOOD for the dimension that is coarser, '...
0488     'because of linear interpolation.'])
0489 return
0490 
0491 % Code to fix this problem incomplete below:
0492 
0493 % %'interp and averaging' = true
0494 % sz = size(in.data);
0495 % if method
0496 %     lp = [sz(2),sz(1)];
0497 % else lp = [sz(1),sz(2)];
0498 % end
0499 %
0500 % % first INTERPOLATE
0501 % if method
0502 %     data = zeros(length(newlat),sz(2),sz(3),class(in.data));
0503 % else data = zeros(sz(1),length(newlon),sz(3),class(in.data));
0504 % end
0505 %
0506 % for i = 1:lp(1)
0507 %     if method
0508 %         data(:,i,:) = interp1(in.lat,squeeze(in.data(:,i,:)),newlat,'linear','extrap');
0509 %     else
0510 %         data(i,:,:) = interp1(in.lon,squeeze(in.data(i,:,:)),newlon,'linear','extrap');
0511 %     end
0512 % end
0513 %
0514 % % then COARSEGRID
0515 % for i = 1:lp(2)
0516 %     if method
0517 %         data(i,:,:) = ;
0518 %     else
0519 %         %
0520 %     end
0521 % end

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