function [out,weights] = resample_geodata(in,field,newgrid,weights)
% 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], {[newlon],[newlat]} (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 [lonstep, latstep]
%                 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'},...
%                        {[ newlons ],[ newlats ]} );
%
%                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 $

errID = ['atmlab:' mfilename ':badInput'];
assert(nargin>2,errID,'Incorrect number of input arguments')
if ~iscell(field), field={field}; end

% rename the structure so nothing is lost after this function
out = in;

% CHECK input
sz                = check_input(in,field,newgrid,errID);

% ASSEMBLE_newgrid
[newlon,newlat]   = assemble_newgrid(in.lon,in.lat,newgrid,errID);

% STANDARDIZE "new grid" internally
gridoptions=struct('lat_descend',false,'lon_descend',false,'lon360',false,'silent',true);
[rg_flags,tmp]    = standardize_geodata(struct('lon',newlon,'lat',newlat),gridoptions);
newlon=tmp.lon;
newlat=tmp.lat;
clear tmp

% LEAVE function if nothing needs to be done
if isequal(in.lat(:),newlat(:)) && isequal(in.lon(:),newlon(:))
    warning(errID,['input grid is identical to requested new grid.\n'...
        'Leaving data untouched'])
    out = in; return
end

% CHECK WHICH test
out.method        = which_test(in.lat,in.lon,newlat,newlon);

% PRINT to screen what is about to happen
if strcmp(out.method,'shift grid')
    logtext(1,'shifting the grid, keeping the same resolution\n')
else
    logtext(1,'Changing grid from %.3fx%.3f to %.3fx%.3f deg by method: %s\n',...
        mean(abs(diff(in.lat))),mean(abs(diff(in.lon))),...
        mean(diff(newlat)),mean(diff(newlon)),....
        out.method)
end

if ~strcmp(out.method,'interpolate') && isfield(in,'interpolation')
    warning(errID,'Will not interpolate, but instead use method: %s',out.method)
elseif strcmp(out.method,'interpolate')
    if isfield(in,'interpolation')
        intm = in.interpolation;
    else
        intm = 'linear';
    end
end

if strcmp(out.method,'coarse grid with area weighting') && ...
    (~exist('weights','var') || isempty(weights))
    % do this only once. This saves time if there are more than one field
    [~,tmp]   = standardize_geodata(struct('lat',in.lat,'lon',in.lon),gridoptions);
    logtext(1,'getting weights\n')
    weights = get_weights_for_coarser_grid_withAreaweighting(struct('lat',tmp.lat,'lon',tmp.lon),newlat,newlon);
end
%% LOOP over FIELDS to regrid

for F = field
    logtext(1,'Processing field: "%s"\n',F{1})
    fsz = size(in.(F{1}));

    % COLLAPSE data onto the 3rd dimesion
    data = in.(F{1})(:,:,:);
    
    % STANDARDIZE "input" data internally
    gridoptions=struct('lat_descend',false,'lon_descend',false,'lon360',false,'dimorder','latlon','silent',true);
    [flags,tmp]   = standardize_geodata(struct('lat',in.lat,'lon',in.lon,'data',data),gridoptions);
    lat  = tmp.lat;
    lon  = tmp.lon;
    data = tmp.data;
    clear tmp
    
    % CHECK for every field
    assert(isequal(fsz(1:2),sz(1:2)),...
        errID,'All input fields for processing don''t have matching lat lon dims')
    
    % USE appropriate METHOD
    switch out.method
        case 'interpolate'
            % successfully tested for 4D data and lon=0:360, lat = end:-1:1
            data = interpolate_grid(struct('lat',lat,'lon',lon,'data',data,'method',intm),newlat,newlon);
        case 'coarse grid'
            % successfully tested for 4D data and lon=0:360, lat = end:-1:1
            data = coarser_grid(struct('lat',lat,'lon',lon,'data',data),newlat,newlon);
        case 'coarse grid with area weighting'
            % successfully tested for 4D data and lon=0:360, lat = end:-1:1
            data = coarser_grid_withAreaweighting(struct('lat',lat,'lon',lon,'data',data),weights,newlat,newlon);
        case 'shift grid'
            % successfully: regular, 4D, 
            data = interpolate_grid(struct('lat',lat,'lon',lon,'data',data,'method',intm),newlat,newlon);
        case 'interp and averaging'
            data = fine_and_coarse(struct('lat',lat,'lon',lon,'data',data),newlat,newlon,true);
        case 'averaging and iterp'
            data = fine_and_coarse(struct('lat',lat,'lon',lon,'data',data),newlat,newlon,false);
        otherwise
            error('atmlab:resample_geodata:Bug','Shouldn''t be here')
    end
    
    % REORDER data to original sorting order
    if flags.lat_descend
        data = data(end:-1:1,:,:);
    end
    if flags.lon_descend
        data = data(:,end:-1:1,:);
    end
    if rg_flags.lon360 % the picked newgrid decides...
        [~,lnindex] = sort(newlon+(newlon<0)*360);
        data = data(:,lnindex,:);
    end
    if strcmp(flags.dimorder,'lonlat')
        data = permute(data,[2,1,3]);
    end

    % RESHAPE data back to original format
    out.(F{1}) = reshape(data,[size(data,1),size(data,2),fsz(3:end)]);
end

out.lat = newlat;
out.lon = newlon;

% REORDER lat lons if need be
if rg_flags.lat_descend
    out.lat = out.lat(end:-1:1);
end
if rg_flags.lon_descend
    out.lon = out.lon(end:-1:1);
end
if rg_flags.lon360
    out.lon = sort(out.lon+(out.lon<0)*360);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUNFUNCTIONS bellow here
%
function sz = check_input(in,field,newgrid,errID)
%% CHECK_input assorted list of input checks

assert(isstruct(in) && ischar(field{1}) && ...
    (isnumeric(newgrid) || iscell(newgrid)),...
    errID,'Incorrect input argument types')

if iscell(newgrid)
    assert(~any([isempty(newgrid{1}),isempty(newgrid{2})]),...
        errID,'Lat or Lon vector is empty')
    
end

assert(~any(~isfield(in,{'lat','lon'})),errID,'"lat,lon" are required fields')
assert(all(ismember(field,fieldnames(in))),errID,'One or more of the specified data fields is missing');
assert(~isequal([size(in.(field{1}),1),size(in.(field{1}),2)],size(in.lat),size(in.lon)),...
    errID,'function only works on regular grids, i.e. data(lat,lon,...) or data(lon,lat,...)')
% assert(max(diff(in.lat))-min(diff(in.lat))<0.1&&max(diff(in.lon))-min(diff(in.lon))<0.1,...
%     errID,'The input data must be equidegree (constant spacing)')
assert( ( isequal(in.lat,sort(in.lat)) || isequal(in.lat,sort(in.lat,'descend')) ) && ... % lats
    ( isequal(in.lon,sort(in.lon)) || isequal(in.lon,sort(in.lon,'descend')) ),...        % lons
    errID,'both lat & lon must be ascending or descending')
sz = size(in.(field{1}));
assert(isequal([length(in.lat),length(in.lon)],[sz(1),sz(2)]) || ...
    isequal([length(in.lat),length(in.lon)],[sz(2),sz(1)]),...
    errID, 'data must have the dimensions data(lat,lon,.....), or data(lon,lat,.....)');

if ~iscell(newgrid)
    assert(min(in.lat)>-90 & max(in.lat)<90 ,errID,'Latitudes are not centers')
    if min(in.lon)<0
        assert(~(all(ismember([-180,180],in.lon))),errID,'Longitudes are probably not centers')
    else
        assert(~(all(ismember([0,360],in.lon))),errID,'Longitudes are probably not centers')
    end
else
   assert(length(newgrid)==2,...
        errID,'if newgrid is a cell, it must contatin 2 vectors {[LATs] , [LONs]}')
    
end


function [newlon,newlat] = assemble_newgrid(LON,LAT,newgrid,errID)
%% ASSEMBLE_newgrid
% in:    LAT             orig latvector
%        LON             orig lonvector
%        newgrid         numeric or cell
%
% out:  newlat,newlon vectors to map to, and logical descn flagging is the
%       are in ascending order or not

if iscell(newgrid)
    % E.g. for mapping one dataset to anothers lat/lons
    newlon = newgrid{1};
    newlat = newgrid{2};
    assert(max(diff(newlat))-min(diff(newlat))<0.1&&max(diff(newlon))-min(diff(newlon))<0.1,...
        errID,'The input newgrid must be equidegree (constant spacing)')
    
else
    % make newlat newlon based on range given by in.lat,in.lon (centered)
    g=newgrid; if isscalar(g), g = [g,g]; end
    dn = mean(abs(diff(LON))); dt = mean(abs(diff(LAT)));
    
    % LATITUDES
    if min(LAT) == LAT(1)
        newlat = LAT(1)-dt/2+g(1)/2 : g(1) : LAT(end)+dt/2-g(1)/2;
        flt = newlat(end)~=LAT(end)+dt/2-g(1)/2;
    else
        newlat = LAT(1)+dt/2-g(1)/2 : -g(1) : LAT(end)-dt/2+g(1)/2;
        flt = newlat(end)~=LAT(end)-dt/2+g(1)/2;
    end
    
    % LONGITUDES
    if min(LON) == LON(1)
        newlon = LON(1)-dn/2+g(2)/2 : g(2) : LON(end)+dn/2-g(2)/2;
        fln = newlon(end)~=LON(end)+dn/2-g(2)/2;
    else
        newlon = LON(1)+dn/2-g(2)/2 : -g(2) : LON(end)-dn/2+g(2)/2;
        fln = newlon(end)~=LON(end)-dn/2+g(2)/2;
    end
    
    % do they FIT ?
    if flt
        warning(errID,['Requested latitude range does not fit into original latitude range\n',...
            'Suggest to give the latitude vector direcly in newgrid = {[lats][lons]}'])
    end
    if fln
        warning(errID,['Requested longitude range does not fit into original longitude range\n',...
            'Suggest to give the longitude vector direcly in newgrid = {[lats][lons]}'])
    end
end

function method = which_test(LAT,LON,newlat,newlon)
%% CHECK_WHICH_test_to_use
% IN:  lat lons of both grids to decide which method to use
% OUT: string name of the method to use

% first test is to see if the edges of the input grid match the edges of
% the new grid

eps = 1e-5; %this corresponds to 1 meter at the equator

% diff averages
dnlt = mean(abs(diff(newlat)));
dnln = mean(abs(diff(newlon)));
dlt  = mean(abs(diff(LAT)));
dln  = mean(abs(diff(LON)));

% temporay vectors
nltb   = newlat - dnlt/2;
nltt   = newlat + dnlt/2;
nlnb   = newlon - dnln/2;
nlnt   = newlon + dnln/2;
ltb    = LAT    - dlt/2;
ltt    = LAT    + dlt/2;
lnb    = LAT    - dln/2;
lnt    = LON    + dln/2;

% LOGICALS

%to a finer grid
cond1   = dlt > dnlt; %length(LAT) < length(newlat);
cond2   = dln > dnln; %length(LON) < length(newlon);

%to a coarser grid
cond3   = dlt < dnlt; %length(LAT) > length(newlat);
cond4   = dln < dnln; %length(LON) > length(newlon);

% same resolution
cond5   = abs(dlt-dnlt)<eps; %length(LAT) == length(newlat);
cond6   = abs(dln-dnln)<eps; %length(LON) == length(newlon);

% exact fit
onLatEdges = ismember(nltb(1),ltb) && ismember(nltt(end),ltt);
onLonEdges = ismember(nlnb(1),lnb) && ismember(nlnt(end),lnt);
fitlt      = dnlt/dlt == floor(dnlt/dlt);         %dnlt = n*dlt; n interger
fitln      = dnln/dln == floor(dnln/dln);         %dnln = n*dln; n interger

if (cond1 && ~cond4) || (cond2 && ~cond3)
    % to a FINER grid
    
    method = 'interpolate';
elseif (cond3 && ~cond2) || (cond4 && ~cond1)
    % to a COARSER grid

    if fitlt && fitln && onLatEdges && onLonEdges
    
        % fit together nicely, this is much faster
        method = 'coarse grid';
    else
        % overlapping gridboxes
        method = 'coarse grid with area weighting';
    end
elseif cond5 && cond6 
    % SAME resolution but DIFFERENT grid
    method = 'shift grid';
elseif cond1 && cond4
    % to a COMBINATION of both methods (INTERP & AVERAGING)
    method = 'interp and averaging';
elseif cond2 && cond3
    % to a COMBINATION of both methods (AVERAGING & INTERP)
    method = 'averaging and iterp';
else
    error('atmlab:resample_geodata:Bug','Shouldn''t be here')
end

function DATA = interpolate_grid(in,newlat,newlon)
%% INTERPOLATE_GRID interpolate data to a smaller grid.

% EXPAND the new lat lons so they encompass orig latlon
data = in.data;

LON=in.lon(:)';
LAT=in.lat(:)';
dn = mean(diff(LON));
dt = mean(diff(LAT));

% make sure they have the same data class
newlon=cast(newlon,class(in.lon));
newlat=cast(newlat,class(in.lat));

while LON(1)>=newlon(1)
    LON = [LON(1)-dn,LON,LON(end)+dn]; %#ok<*AGROW>
    data = [data(:,end,:),data,data(:,1,:)];
end
while LAT(1)>=newlat(1)
    LAT = [LAT(1)-dt,LAT,LAT(end)+dt];
    data = [data(1,:,:);data;data(end,:,:)];
end

% do the INTERPOLATION
[X,Y] = meshgrid(LON,LAT);
[XI,YI] = meshgrid(newlon,newlat);
if islogical(data)
    DATA = false([length(newlat),length(newlon),size(data,3)]);
else
    DATA = zeros([length(newlat),length(newlon),size(data,3)],class(data));
end

logtext(1,'Using interpolation method: %s\n',in.method)
for i = 1:size(data,3)
    DATA(:,:,i) = interp2(X,Y,data(:,:,i),XI,YI,in.method);
end

function data = coarser_grid(in,lat,lon)
%% coarser_grid by averaging.
%
% Purpose: To make gridded data more sparse if the grids fit nicely into each other.
% NOTE: Biases may be introduced because I am using nanmean. The solution
%       is to also output the number of values used in the averaging for
%       each grid (not doing that right now. Add feature if needed)

% PUT lat/lons LAST and make lat lons both size(lat,lon)
tmpdata = permute(in.data,[3,1,2]); 
tmplon  = repmat(in.lon,length(in.lat),1);
tmplat  = repmat(in.lat,1,length(in.lon));

% COLLAPSE data (transposed), lat, lon to simulate UNGRIDDED data
tmpdata = binning_fast(struct('data',tmpdata(:,:)','lon',tmplon(:),'lat',tmplat(:),...
    'newlon',lon,'newlat',lat));

% TAKE the grid MEAN for each level (e.g. {[4x288]} -> {[1x288]}).
if any(isnan(tmpdata(:))), warning(['atmlab:' mfilename ':calculation'],...
        'There are NaNs in the data, and these may lead to biases on the coarser grid'); end
tmpdata = binned_statistics(tmpdata,{@nanmean});

if ismatrix(in.data)
    data = cell2mat(tmpdata.mean);
elseif ndims(in.data)==3
    % make into MATRIX(lat,lat,...)
    % this loop is because cell2mat doesn't work if cells contain vectors. It
    % doesn't work in conjuction with repmat, it just scrambled
    data = zeros(length(lat),length(lon),size(in.data,3));
    for i = 1:length(lat)
        for j = 1:length(lon)
            data(i,j,:) = tmpdata.mean{i,j};
        end
    end
end

function weights = get_weights_for_coarser_grid_withAreaweighting(in,newlat,newlon)

% DETERMINE new grid BOUNDARIES from grid spacing
gn = mean(diff(newlon));    gt = mean(diff(newlat));
dn = mean(diff(in.lon)); dt = mean(diff(in.lat));
if isnan(gt), gt = (in.lat(end)-in.lat(1) + dt);end %incase there's only one latitude
if isnan(gn), gn = (in.lon(end)-in.lon(1) + dn);end %incase there's only one longitude

lattop = newlat+gt/2+dt;
latbot = newlat-gt/2;
lontop = newlon+gn/2+dn;
lonbot = newlon-gn/2;

weights = cell(length(newlat),length(newlon));

%LOOP over larger grid
for ii = 1:length(newlat)
    for jj=1:length(newlon)
        
        % GET lats lons in or tangent to grid(i,j)
        ilon = in.lon >= lonbot(jj) & in.lon < lontop(jj);
        ilat = in.lat >= latbot(ii) & in.lat < lattop(ii);
        
        % find WEIGHTS
        % A = [x,y,width,height] x,y = bottom corner
        A = [lonbot(jj),latbot(ii),gn,gt];
        [x,y] = meshgrid(in.lon(ilon)-dn/2,in.lat(ilat)-dt/2);
        B = [ x(:), y(:), repmat(dn,length(y(:)),1), repmat(dt,length(x(:)),1)];
        weights{ii,jj} = rectint(A,B);

    end
end

function data = coarser_grid_withAreaweighting(in,weights,lat,lon)
%% COARSER_GRID make a coarse grid from a finer grid
% method: Uses area weighting for grids that don't fit nicely into each
% other and assumes local plane assumption, which is good as long as the
% grids are not too large (see help resample_data).
% NOTE: Biases may be introduced because I am using nanmean. The solution
%       is to also output the number of values used in the averaging for
%       each grid (not doing that right now. Add feature if needed)


% preallocate
if ~islogical(in.data)
    data = zeros([length(lat),length(lon),size(in.data,3)],class(in.data));
else
    data = false([length(lat),length(lon),size(in.data,3)]);
end

polyWeightSum = @(d,w,i)(nansum(single(d(i)).*w(i))./nansum(w(i)));

% DETERMINE new grid BOUNDARIES from grid spacing
gn = mean(diff(lon)); gt = mean(diff(lat));
dn = mean(diff(in.lon)); dt = mean(diff(in.lat));
if isnan(gt), gt = (in.lat(end)-in.lat(1) + dt);end %incase there's only one latitude
if isnan(gn), gn = (in.lon(end)-in.lon(1) + dn);end %incase there's only one longitude

lattop = lat+gt/2+dt;
latbot = lat-gt/2;
lontop = lon+gn/2+dn;
lonbot = lon-gn/2;

%LOOP over larger grid
for ii = 1:length(lat)
    for jj=1:length(lon)
        
        % GET lats lons in or tangent to grid(i,j)
        ilon = in.lon >= lonbot(jj) & in.lon < lontop(jj);
        ilat = in.lat >= latbot(ii) & in.lat < lattop(ii);

        % get data SUBSET
        d = in.data(ilat,ilon,:);
        
        % COLLAPSE data on the lat/lon dimensions to match output from rectint
        d = permute(d,[3,1,2]);
        
        for kk = 1:size(d,1)
            data(ii,jj,kk)    = polyWeightSum(d(kk,:),weights{ii,jj},~isnan(d(kk,:)));
        end
    end
end

% function data = fine_and_coarse(in,newlat,newlon,method)
% 
% data = interpolate_grid(in,newlat,newlon);
% warning('atmlab:resample_geodata:incomplete',...
%     ['Using interpolation (interp2) for grids that are finer in one dimension '...
%     'and coarser in the other. This is NO GOOD for the dimension that is coarser, '...
%     'because of linear interpolation.'])
% 
% return

% Code to fix this problem incomplete below:

% %'interp and averaging' = true
% sz = size(in.data);
% if method
%     lp = [sz(2),sz(1)];
% else lp = [sz(1),sz(2)];
% end
% 
% % first INTERPOLATE
% if method
%     data = zeros(length(newlat),sz(2),sz(3),class(in.data));
% else data = zeros(sz(1),length(newlon),sz(3),class(in.data));
% end
% 
% for i = 1:lp(1)
%     if method
%         data(:,i,:) = interp1(in.lat,squeeze(in.data(:,i,:)),newlat,'linear','extrap');
%     else
%         data(i,:,:) = interp1(in.lon,squeeze(in.data(i,:,:)),newlon,'linear','extrap');
%     end
% end
% 
% % then COARSEGRID
% for i = 1:lp(2)
%     if method
%         data(i,:,:) = ;
%     else
%         %
%     end
% end