/[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.3 - (hide annotations) (download)
Thu Apr 14 18:50:14 2011 UTC (14 years, 3 months ago) by roquet
Branch: MAIN
Changes since 1.2: +19 -23 lines
speed optimizations

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

  ViewVC Help
Powered by ViewVC 1.1.22