profilesMatlabProcessing/profiles_devel/profiles_subgrid_stats.m

function [myStat]=profiles_subgrid_stats(KK,VV,TYPE,SUB,COORD);
%KK level choice
%VV variable choice
%TYPE 'obs' 'estim' or 'anom'
%SUB is the subsampling rate (box width=SUB)
%COORD is the coordinate type

%global variables:
gcmfaces_global;
global lon lat obs point;
global COORDOld; if isempty(COORDOld); COORDOld='x'; end;
global VVOld; if isempty(VVOld); VVOld='x'; end;
global TYPEOld; if isempty(TYPEOld); TYPEOld='x'; end;
global KKOld; if isempty(KKOld); KKOld=0; end;

%choice of time/depth ranges
%RC=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RC'));
%RF=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RF'));
%depth0=-RF(1:end-1); depth0(2:end-1)=depth0(1:end-2);
%depth1=-RF(2:end); depth1(2:end-1)=depth1(3:end);
%if KK>=1;
%    depth0=depth0(KK); depth1=depth1(KK);%choice of depth range
%else;
%    error('not implemented');
%end;
%
date0=datenum(1950,1,1); date1=datenum(2049,12,31);%choice of time range

%get the grid:
dirGrid='/net/nares/raid10/gforget/2012julyIters/GRID/';
grid_load(dirGrid,5,'compact');
gcmfaces_bindata;

%get the data:
dirData='./';
listData=dir([dirData 'argo_2may13_set*.nc']);
listData={listData(:).name};
suffOut='argo';

%
test1=~strcmp(VV,VVOld)|~strcmp(TYPE,TYPEOld)|KK~=KKOld|~strcmp(COORD,COORDOld);
if test1;
    if strcmp(TYPE,'anom');
    [MITprof]=MITprof_stats_load(dirData,listData,VV,1);
    elseif strcmp(TYPE,'estim');
    [MITprof]=MITprof_stats_load(dirData,listData,VV,['prof_' VV 'estim']);
    else;
    [MITprof]=MITprof_stats_load(dirData,listData,VV,['prof_' VV]);
    end;
    %in DRHODR case, switch to log10 :
    if strcmp(VV,'DRHODR'); MITprof.prof=log10(MITprof.prof); end;
    %mask out values that are not in year range:
    ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
    MITprof.prof(ii,:)=NaN;
    %restrict to depth range of interest:
    %kk=find(MITprof.prof_depth>=depth0&MITprof.prof_depth<=depth1);
    kk=KK;
    lon=MITprof.prof_lon; lat=MITprof.prof_lat; obs=MITprof.prof(:,kk);
end;
%
KKOld=KK;
VVOld=VV;
TYPEOld=TYPE;
COORDOld=COORD;

%get indices in full grid: (not necessarily the ecco_v4 one)
point=gcmfaces_bindata(lon,lat);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%check nearest neighbor mapping:
% [tmp1,tmp2]=gcmfaces_bindata(lon,lat,lon);
% figureL; qwckplot(log10(tmp2));
% figureL; qwckplot(tmp1./tmp2);

%map of full grid indices: (consistent with prof_point2)
indGrid=convert2array(mygrid.XC);
indGrid(:)=[1:length(indGrid(:))];
indGrid=convert2array(indGrid);

%reduce grid and map reduced grid indices:
indBox=mygrid.XC;
if SUB==0;
    %global computation
    indBox(:)=1;
else;
    %regional computation
    boxMax=0;
    for iFace=1:mygrid.nFaces;
        tmp1=indBox{iFace};
        tmp3=ceil([1:size(tmp1,1)]'/SUB);
        tmp4=(ceil([1:size(tmp1,2)]/SUB)-1)*max(tmp3);
        tmp5=tmp3*ones(1,size(tmp1,2))+ones(size(tmp1,1),1)*tmp4;
        indBox{iFace}=tmp5+boxMax;
        boxMax=boxMax+max(tmp5(:));
    end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%main computational loop: stats for each region in indBox
%
indBox=convert2array(indBox);
box=repmat(indBox(point),[1 size(obs,2)]);
box(isnan(obs))=NaN;
boxList=unique(indBox(:));
boxList=boxList(find(~isnan(boxList)));
%
tmp1=convert2array(0*mygrid.XC);
myStat.mea=tmp1; myStat.prc90=tmp1; myStat.med=tmp1; myStat.prc10=tmp1; myStat.nb=tmp1;
myStat.std=tmp1; myStat.iqr=tmp1; myStat.mad=tmp1;
%
for ii=boxList';
    if mod(ii,1000)==0;
        [KK ii length(boxList)]
    end;
    jj=find(box==ii);
    if length(jj)>=10;
        tmpStat=myStats(obs(jj));%need reduced params
        jj=find(indBox==ii);
        tmpList=fieldnames(tmpStat);
        for pp=1:length(tmpList);
          eval(['myStat.' tmpList{pp} '(jj)=tmpStat.' tmpList{pp} ';']);
        end;
    end;
end;
%
tmpList=fieldnames(myStat);
for pp=1:length(tmpList); 
  eval(['myStat.' tmpList{pp} '=convert2array(myStat.' tmpList{pp} ');']); 
end;

eval(['save ' dirData 'stats/' VV '_k' num2str(KK) '_' num2str(SUB) '.mat myStat;']);

function [myStat]=myStats(obs);

myStat.nb=sum(~isnan( obs ));

myStat.mea=mean(obs);%sample mean
myStat.prc10=prctile(obs,10);
myStat.med=median(obs);
myStat.prc90=prctile(obs,90);

myStat.std=std(obs);%sample standard deviation
myStat.iqr=0.7413*iqr(obs);%intequartile range estimate of std
myStat.mad=1.4826*mad(obs,1);%median absolute difference estimate of std

1	gforget	1.3	function [myStat]=profiles_subgrid_stats(KK,VV,TYPE,SUB,COORD);
2			%KK level choice
3			%VV variable choice
4			%TYPE 'obs' 'estim' or 'anom'
5			%SUB is the subsampling rate (box width=SUB)
6			%COORD is the coordinate type
7	gforget	1.1
8			%global variables:
9			gcmfaces_global;
10			global lon lat obs point;
11	gforget	1.3	global COORDOld; if isempty(COORDOld); COORDOld='x'; end;
12			global VVOld; if isempty(VVOld); VVOld='x'; end;
13			global TYPEOld; if isempty(TYPEOld); TYPEOld='x'; end;
14			global KKOld; if isempty(KKOld); KKOld=0; end;
15	gforget	1.1
16			%choice of time/depth ranges
17	gforget	1.3	%RC=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RC'));
18			%RF=squeeze(rdmds('/net/weddell/raid3/gforget/ecco_v4/GRID/RF'));
19			%depth0=-RF(1:end-1); depth0(2:end-1)=depth0(1:end-2);
20			%depth1=-RF(2:end); depth1(2:end-1)=depth1(3:end);
21			%if KK>=1;
22			% depth0=depth0(KK); depth1=depth1(KK);%choice of depth range
23			%else;
24			% error('not implemented');
25			%end;
26	gforget	1.1	%
27			date0=datenum(1950,1,1); date1=datenum(2049,12,31);%choice of time range
28
29			%get the grid:
30	gforget	1.3	dirGrid='/net/nares/raid10/gforget/2012julyIters/GRID/';
31			grid_load(dirGrid,5,'compact');
32			gcmfaces_bindata;
33	gforget	1.1
34			%get the data:
35	gforget	1.3	dirData='./';
36			listData=dir([dirData 'argo_2may13_set*.nc']);
37			listData={listData(:).name};
38			suffOut='argo';
39
40	gforget	1.1	%
41	gforget	1.3	test1=~strcmp(VV,VVOld)\|~strcmp(TYPE,TYPEOld)\|KK~=KKOld\|~strcmp(COORD,COORDOld);
42	gforget	1.1	if test1;
43	gforget	1.3	if strcmp(TYPE,'anom');
44			[MITprof]=MITprof_stats_load(dirData,listData,VV,1);
45			elseif strcmp(TYPE,'estim');
46			[MITprof]=MITprof_stats_load(dirData,listData,VV,['prof_' VV 'estim']);
47	gforget	1.2	else;
48	gforget	1.3	[MITprof]=MITprof_stats_load(dirData,listData,VV,['prof_' VV]);
49	gforget	1.2	end;
50	gforget	1.1	%in DRHODR case, switch to log10 :
51	gforget	1.3	if strcmp(VV,'DRHODR'); MITprof.prof=log10(MITprof.prof); end;
52	gforget	1.1	%mask out values that are not in year range:
53			ii=find(MITprof.prof_date<date0\|MITprof.prof_date>date1);
54			MITprof.prof(ii,:)=NaN;
55			%restrict to depth range of interest:
56	gforget	1.3	%kk=find(MITprof.prof_depth>=depth0&MITprof.prof_depth<=depth1);
57			kk=KK;
58	gforget	1.1	lon=MITprof.prof_lon; lat=MITprof.prof_lat; obs=MITprof.prof(:,kk);
59			end;
60			%
61	gforget	1.3	KKOld=KK;
62			VVOld=VV;
63			TYPEOld=TYPE;
64			COORDOld=COORD;
65	gforget	1.1
66			%get indices in full grid: (not necessarily the ecco_v4 one)
67			point=gcmfaces_bindata(lon,lat);
68
69			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70
71			%check nearest neighbor mapping:
72			% [tmp1,tmp2]=gcmfaces_bindata(lon,lat,lon);
73			% figureL; qwckplot(log10(tmp2));
74			% figureL; qwckplot(tmp1./tmp2);
75
76			%map of full grid indices: (consistent with prof_point2)
77			indGrid=convert2array(mygrid.XC);
78			indGrid(:)=[1:length(indGrid(:))];
79			indGrid=convert2array(indGrid);
80
81			%reduce grid and map reduced grid indices:
82			indBox=mygrid.XC;
83	gforget	1.3	if SUB==0;
84	gforget	1.1	%global computation
85			indBox(:)=1;
86			else;
87			%regional computation
88			boxMax=0;
89			for iFace=1:mygrid.nFaces;
90			tmp1=indBox{iFace};
91	gforget	1.3	tmp3=ceil([1:size(tmp1,1)]'/SUB);
92			tmp4=(ceil([1:size(tmp1,2)]/SUB)-1)*max(tmp3);
93	gforget	1.1	tmp5=tmp3ones(1,size(tmp1,2))+ones(size(tmp1,1),1)tmp4;
94			indBox{iFace}=tmp5+boxMax;
95			boxMax=boxMax+max(tmp5(:));
96			end;
97			end;
98
99			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100
101			%main computational loop: stats for each region in indBox
102			%
103			indBox=convert2array(indBox);
104			box=repmat(indBox(point),[1 size(obs,2)]);
105			box(isnan(obs))=NaN;
106			boxList=unique(indBox(:));
107			boxList=boxList(find(~isnan(boxList)));
108			%
109			tmp1=convert2array(0*mygrid.XC);
110			myStat.mea=tmp1; myStat.prc90=tmp1; myStat.med=tmp1; myStat.prc10=tmp1; myStat.nb=tmp1;
111	gforget	1.2	myStat.std=tmp1; myStat.iqr=tmp1; myStat.mad=tmp1;
112	gforget	1.1	%
113			for ii=boxList';
114			if mod(ii,1000)==0;
115	gforget	1.3	[KK ii length(boxList)]
116	gforget	1.1	end;
117			jj=find(box==ii);
118			if length(jj)>=10;
119	gforget	1.2	tmpStat=myStats(obs(jj));%need reduced params
120	gforget	1.1	jj=find(indBox==ii);
121	gforget	1.2	tmpList=fieldnames(tmpStat);
122			for pp=1:length(tmpList);
123			eval(['myStat.' tmpList{pp} '(jj)=tmpStat.' tmpList{pp} ';']);
124			end;
125	gforget	1.1	end;
126			end;
127			%
128	gforget	1.2	tmpList=fieldnames(myStat);
129			for pp=1:length(tmpList);
130			eval(['myStat.' tmpList{pp} '=convert2array(myStat.' tmpList{pp} ');']);
131			end;
132
133	gforget	1.3	eval(['save ' dirData 'stats/' VV '_k' num2str(KK) '_' num2str(SUB) '.mat myStat;']);
134
135	gforget	1.2	function [myStat]=myStats(obs);
136
137			myStat.nb=sum(~isnan( obs ));
138
139			myStat.mea=mean(obs);%sample mean
140			myStat.prc10=prctile(obs,10);
141			myStat.med=median(obs);
142			myStat.prc90=prctile(obs,90);
143	gforget	1.1
144	gforget	1.2	myStat.std=std(obs);%sample standard deviation
145			myStat.iqr=0.7413*iqr(obs);%intequartile range estimate of std
146			myStat.mad=1.4826*mad(obs,1);%median absolute difference estimate of std
147	gforget	1.1