/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_IO_external/profiles_read_odv.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_IO_external/profiles_read_odv.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Sat May 7 20:48:47 2011 UTC (14 years, 2 months ago) by roquet
Branch: MAIN
Changes since 1.3: +31 -19 lines
simplify profiles_IO_external interface:
init loading: dataset=profiles_read_xxx(dataset,nf,0);
then, to load m-th profile: profile=profiles_read_xxx(dataset,nf,m);

1 roquet 1.2 function [varargout]=profile_read_odv(dataset,nf,m,varargin)
2 roquet 1.1 % read seal data in the odv spreadsheet format
3     %
4     % if m=0 :
5     % return the information on the nf-th file referenced in
6     % dataset.fileInList.
7 roquet 1.4 % dataset=profile_read_odv(dataset,nf,0);
8     % dataset.nprofiles: number of profiles in the nf-th file
9     % dataset.Iprofile: index referencing the profile for each dataline
10 roquet 1.1 %
11     % if m~=0 :
12     % return the m-th profile from the nf-th file referenced in
13     % dataset.fileInList.
14 roquet 1.4 % profileCur=profile_read_odv(dataset,nf,m);
15 roquet 1.1 %
16     % profileCur = (m-th profile of nf-th file)
17     % pnum_txt: '5900841'
18     % ymd: 20060101
19     % hms: 3207
20     % lat: -46.709
21     % lon: 137.501
22     % direc: 1
23     % t: [1x115 single]
24     % s: [1x115 single]
25     % p: [1x115 single]
26     % t_ERR: [1x115 single]
27     % s_ERR: [1x115 single]
28     % PorZisBAD: 0
29     %
30 roquet 1.4 % IMPORTANT: Use three global variables: pnum Date and Data.
31 roquet 1.1
32 roquet 1.4 global pnum Date Data
33 roquet 1.1
34     fileIn=[dataset.dirIn dataset.fileInList(nf).name];
35    
36     if m==0;
37    
38     % return the number of profiles in the file
39     % also build an index of first-line prosition of profiles: index_prof
40     fid_cur=fopen(fileIn,'rt');
41     nprofiles=0;index_prof=[];
42    
43     % jump comment line at the beginning
44     line_cur = fgetl(fid_cur);
45     while ~feof(fid_cur) & strcmp(line_cur(1:2),'//')
46     line_cur = fgetl(fid_cur);
47     continue
48     end
49    
50     % analyse data description line
51     I=[0 find(double(line_cur)==9) length(line_cur)+1];
52     format=[];
53     % column index for (in this order) :
54     % Cruise, Station, Date, Lon, Lat, Dep, Dep_qv, Temp,
55     % Temp_qv, Sal, Sal_qv
56     for ii=1:length(I)-1,
57     field=line_cur(I(ii)+1:I(ii+1)-1);
58     field=field(1:min([6 length(field)]));
59     switch field
60     case 'Cruise'
61     format(1)=ii;
62     case 'Statio'
63     format(2)=ii;
64     case 'yyyy-m'
65     format(3)=ii;
66     case {'Longit'}
67     format(4)=ii;
68     case {'Latitu'}
69     format(5)=ii;
70     case {'Depth '}
71     format([6 7])=[ii ii+1];
72     case {'Temper'}
73     format([8 9])=[ii ii+1];
74     case {'Salini'}
75     format([10 11])=[ii ii+1];
76     end
77     end
78    
79     % read end of the textfile
80 roquet 1.2 nlines=0;nbuffer=100000;niter=0;
81 roquet 1.1 while ~feof(fid_cur);
82 roquet 1.2 niter=niter+1;
83     data_odv=cell(length(format)-1,nbuffer);
84     for ii=1:nbuffer,
85     tline=fgets(fid_cur);
86 roquet 1.1 nlines=nlines+1;
87     I_sep=[0 find(double(tline)==9)];
88 roquet 1.2 data_odv{1,ii}=[ tline(I_sep(format(1))+1:I_sep(format(1)+1)-1) '//' ...
89     tline(I_sep(format(2))+1:I_sep(format(2)+1)-1) ];
90     data_odv{2,ii}=tline(I_sep(format(3))+1:I_sep(format(3)+1)-1);
91     for kk=4:length(format),
92     I=I_sep(format(kk))+1:I_sep(format(kk)+1)-1;
93     if length(I)>1,
94     data_odv{kk-1,ii}=tline(I);
95     elseif isempty(I),
96     data_odv{kk-1,ii}='99999';
97     else
98     val=tline(I);
99     if val==0,
100     data_odv{kk-1,ii}='99999';
101     else
102     data_odv{kk-1,ii}=val;
103     end
104     end
105 roquet 1.1 end
106     if feof(fid_cur), break, end
107     end
108 roquet 1.2 data_odv=data_odv(:,1:ii);
109     data=data_odv(3:length(format)-1,:);
110     datanum=sscanf(sprintf('%s,',data{:}), '%g,');
111     datanum(datanum==99999)=NaN;
112 roquet 1.4
113     %no salinity case
114     if dataset.inclS,
115     datanum=reshape(datanum(:),8,ii);
116     else
117     datanum=reshape(datanum(:),6,ii);
118     end
119 roquet 1.2
120     fprintf('%d lines\n',nlines);
121     if niter==1,
122 roquet 1.3 pnum=data_odv(1,:);
123 roquet 1.2 Date=data_odv(2,:);
124     Data=datanum;
125 roquet 1.1 else
126 roquet 1.3 pnum=[pnum data_odv(1,:)];
127 roquet 1.2 Date=[Date data_odv(2,:)];
128     Data=[Data datanum];
129 roquet 1.1 end
130    
131     end
132 roquet 1.2 fclose(fid_cur);
133    
134 roquet 1.1 % get rid of empty indicator string
135 roquet 1.3 [pnum_txt_list,Iprof]=unique(pnum);
136     for kk=1:length(pnum),
137 roquet 1.2 if strcmp(pnum_txt_list{kk},'//'),
138     pnum_txt_list=pnum_txt_list(setdiff(1:length(pnum_txt_list),kk));
139 roquet 1.1 Iprof(kk)=[];
140     break
141     end
142     end
143    
144     % sort index of profile
145     [Iprof_sort,J]=sort(Iprof);
146 roquet 1.2 pnum_txt_sort=pnum_txt_list(J);
147 roquet 1.1
148     % write index of profile
149     Iprof_sort=[Iprof_sort nlines+1];
150 roquet 1.3 Iprofile=zeros(nlines,1);
151 roquet 1.1 for kk=1:length(Iprof_sort)-1,
152     for ii=Iprof_sort(kk):Iprof_sort(kk+1)-1,
153 roquet 1.3 Iprofile(ii)=kk;
154 roquet 1.1 end
155     end
156 roquet 1.4 dataset.Iprofile=Iprofile;
157 roquet 1.2
158 roquet 1.1 % number of profiles
159 roquet 1.4 dataset.nprofiles=length(Iprof);
160 roquet 1.1
161 roquet 1.4 varargout(1) = {dataset};
162 roquet 1.1
163    
164     else;%if m==0;
165    
166     % profile coordinates
167 roquet 1.4 Iprof=find(dataset.Iprofile==m);
168 roquet 1.3 pnum_txt=pnum{Iprof(1)};
169     lon=Data(1,Iprof(1));
170     lat=Data(2,Iprof(1));
171     strdate=Date{Iprof(1)};
172 roquet 1.1 if length(strdate)<10, varargout = {[]}; return; end
173     ymd=str2num(strdate([1:4 6:7 9:10]));
174     hms=0;
175     if length(strdate)==19, hms=str2num(strdate([12:13 15:16 18:19])); end
176     if length(strdate)==16, hms=str2num(strdate([12:13 15:16]))*100; end
177     if length(strdate)==13, hms=str2num(strdate([12:13]))*10000; end
178    
179     %case when necessary information missing: dont retrieve profile
180     if isempty(lon)|isempty(lat)|isempty(ymd), varargout = {[]}; return; end
181     if lon < 0; lon=lon+360; end;
182    
183     % get T/S data
184 roquet 1.3 z=Data(3,Iprof);
185     z_qv=Data(4,Iprof);z_qv(isempty(z_qv))=1;
186     t=Data(5,Iprof);
187     t_qv=Data(6,Iprof);t_qv(isempty(t_qv))=1;
188 roquet 1.4 if dataset.inclS,
189     s=Data(7,Iprof);
190     s_qv=Data(8,Iprof);s_qv(isempty(s_qv))=1;
191     else
192     s=t*NaN; s_qv=t*NaN; s_ERR=t*NaN;
193     end
194    
195 roquet 1.1 I=find(isnan(z)|z_qv~=0);
196 roquet 1.4 z(I)=[]; t(I)=[]; t_qv(I)=[];
197     t(t_qv~=0)=NaN;
198     t_ERR=t*0;
199     if dataset.inclS,
200     s_qv(I)=[]; s(I)=[];
201     s(s_qv~=0)=NaN; s_ERR=s*0;
202     end
203 roquet 1.1
204     direc=2;
205     PorZisBAD=0;
206    
207     profileCur.pnum_txt=pnum_txt;
208     profileCur.ymd=ymd; profileCur.hms=hms;
209     profileCur.lat=lat; profileCur.lon=lon;
210     profileCur.direc=direc;
211     profileCur.PorZisBAD=PorZisBAD;
212     profileCur.t=t;
213     profileCur.s=s;
214     profileCur.z=z;
215     profileCur.t_ERR=t_ERR;
216     profileCur.s_ERR=s_ERR;
217    
218     varargout = {profileCur};
219    
220     end;
221    
222    
223    

  ViewVC Help
Powered by ViewVC 1.1.22