source: issm/trunk/src/m/model/plot/processdata.m@ 7098

function [data datatype]=processdata(md,data,options);
%PROCESSDATA - process data to be plotted
%
%   datatype = 1 -> elements
%   datatype = 2 -> nodes
%   datatype = 3 -> node quivers
%   datatype = 4 -> patch
%
%   Usage:
%      [data datatype]=processdata(md,data,options);
%
%   See also: PLOTMODEL, PROCESSMESH

%check format
if (iscell(data) | isempty(data) | length(data)==0 | (length(data)==1 & ~isstruct(data) & isnan(data))),
        error('plotmodel error message: data provided is empty');
end

%Process Patch
if isstruct(data) 
        if (isfield(data,'index') & isfield(data,'value')),
                if data.interpolation(1)==P1Enum(),
                        data=data.value;
                        data=data';
                        data=data(:);
                        datatype=4;
                elseif data.interpolation(1)==P0Enum(),
                        data=data.value;
                        datatype=5;
                else
                        error(['interpolation ' data.interpolation(1) ' not supported yet']);
                end
        else
                error('structure other than Patch not supported yet');
        end
else
        %initialize datatype
        datatype=0;
end

%get datatype
datasize=size(data);

%non patch processing
if datatype~=4 & datatype~=5,

        %transpose data if necessary
        if (size(data,2) > size(data,1)),
                data=data';
        end
        datasize=size(data);

        %convert to double if necessary
        if ~isnumeric(data);
                disp('processdata info message: data is not numeric (logical?). Converted to double');
                data=double(data);
        end

        %check length
        if datasize(1)~=md.numberofgrids & datasize(1)~=md.numberofelements & datasize(1)~=md.numberofgrids*6 & (md.dim==3 & ~(datasize(1)==md.numberofelements2d | datasize(1)==md.numberofgrids2d))
                error('plotmodel error message: data not supported yet');
        end

        %quiver?
        if datasize(2)>1,
                datatype=3;

                %check number of columns, add zeros if necessary,
                if (md.dim==3)
                        if datasize(2)==2,
                                data=[data, zeros(datasize(1),1)];
                        elseif datasize(2)~=3,
                                error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot');
                        end
                %elseif ((md.dim==2) & datasize(2)~=2),
                %       error('plotmodel error message: data provided should have 2 columns for quiver plot, and 1 for regular plot');
                end
        end

        %treat the case datasize(1)=6*grids
        if datasize(1)==6*md.numberofgrids
                %keep the only norm of data
                data1=data(1:6:md.numberofgrids*6,:);
                data2=data(2:6:md.numberofgrids*6,:);
                data=sqrt(data1.^2+data2.^2);
                datasize(1)=md.numberofgrids;
                %---> go to grid data
        end

        %treat the case datasize(1)=grids2d
        if (md.dim==3 & datasize(1)==md.numberofgrids2d),
                data=project3d(md,data,'node');
                datasize(1)=md.numberofgrids;
                %---> go to grid data
        end

        %treat the case datasize(1)=grids2d
        if (md.dim==3 & datasize(1)==md.numberofelements2d),
                data=project3d(md,data,'element');
                datasize(1)=md.numberofelements;
                %---> go to grid data
        end

        %smoothing?
        if exist(options,'smooth')
                data=averaging(md,data,getfieldvalue(options,'smooth'));
                datasize(1)=md.numberofgrids;
                %---> go to grid data
        end
end

%element data
if (datasize(1)==md.numberofelements & datasize(2)==1),

        %Initialize datatype if non patch
        if datatype~=4 & datatype~=5,
                datatype=1;
        end

        %ice sheet only?
        if getfieldvalue(options,'iceshelf',1)==0,
                data(find(md.elementoniceshelf),:)=NaN;
        end
        %ice shelf only?
        if getfieldvalue(options,'icesheet',1)==0,
                data(find(~md.elementoniceshelf),:)=NaN;
        end
        %no water?
        if getfieldvalue(options,'water',1)==0,
                data(find(md.elementonwater),:)=NaN;
        end
        %log?
        if exist(options,'log'),
                bounds=getfieldvalue(options,'caxis',[min(data(:)) max(data(:))]);
                data(find(data<bounds(1)))=bounds(1);
                if any(data<=0),
                        error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
                end
                pos=find(~isnan(data));
                data(pos)=log(data(pos))/log(getfieldvalue(options,'log'));
        end
end

%grid data
if (datasize(1)==md.numberofgrids & datasize(2)==1),
        datatype=2;
        %ice sheet only?
        if getfieldvalue(options,'iceshelf',1)==0,
                pos=find(md.gridoniceshelf);
                data(pos,:)=NaN;
        end
        %ice shelf only?
        if getfieldvalue(options,'icesheet',1)==0,
                pos=find(md.gridonicesheet);
                data(pos,:)=NaN;
        end
        %no water?
        if getfieldvalue(options,'water',1)==0,
                pos=find(md.gridonwater);
                data(pos,:)=NaN;
        end
        %log?
        if exist(options,'log'),
                %if any(data<=0),
                %       error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
                %end
                data=log(data)/log(getfieldvalue(options,'log'));
        end

end

%layer projection? 
if getfieldvalue(options,'layer',0)>=1,
        data=project2d(md,data,getfieldvalue(options,'layer')); %project onto 2d mesh
end

%control arrow density if quiverplot
if datatype==3 & exist(options,'density')
        databak=data;
        data=NaN*ones(datasize);
        density=getfieldvalue(options,'density');
        data(1:density:end,:)=databak(1:density:end,:);
        clear databak
end