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
if strcmp(COORD,'depth');
  %KK is ecco v4 level index; kk is corresponding MITprof level indices
  RC=squeeze(rdmds([myenv.gcmfaces_dir '/sample_input/GRIDv4/RC']));
  RF=squeeze(rdmds([myenv.gcmfaces_dir '/sample_input/GRIDv4/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;
end;
%
date0=datenum(1950,1,1); date1=datenum(2049,12,31);%choice of time range

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