-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_zamg.m
216 lines (202 loc) · 8.59 KB
/
read_zamg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
function [stationNo, dateAll, dataAll, varargout] = read_zamg(dirName, completeData)
% This function reads observational or simulated data in the format used by
% ZAMG.
% The simulated data is complete data (completeData==1), the observational
% data is not complete. To speed up reading, the function creates a matlab
% datafile with the data (matlabDataFileDirName) at the end and if this file is present at the
% beginning, it will read this matlab file. To override this, delete the
% matlab datafile or set readAnew = 1.
% Note that DecimalDate is given in years and is not exact, should only be used for plotting.
if ( ( nargin < 1 ) || isempty(dirName) )
number = [1, 2, 3, 4, 5, 7];
number = number(6);
dirName = ['/data/zamg_humidity/Netzwerke/Netzwerk', num2str(number), '_voll_hom_vv/'];
end
if isnumeric(dirName)
number = [1, 2, 3, 4, 5, 7];
number = number(dirName);
dirName = ['/data/zamg_humidity/Netzwerke/Netzwerk', num2str(number), '_voll_hom_vv/'];
end
if ( ( nargin < 2 ) || isempty(completeData) )
completeData = 0;
end
if ( completeData )
[stationNo, dateAll, dataAll] = read_zamg_complete(dirName);
else
matlabDataFileDirName = fullfile(dirName, 'data.mat');
readAnew = 0;
if ( (exist(matlabDataFileDirName, 'file') > 0 ) && (readAnew == 0) )
load(matlabDataFileDirName)
else
files = dir(fullfile(dirName, 'Complete_Rel_Feuchte_*.txt'));
noFiles = numel(files);
fileContainsMetaData = 0;
% First make vectors containing dates running for the full period of the network
% Read files to determine begin and end year (only full years are used)
beginDate = +inf;
endDate = -inf;
for iFile = 1:noFiles
fileName = fullfile(dirName, files(iFile).name);
fid = fopen(fileName);
tline = fgetl(fid); %#ok<NASGU>
tline = fgetl(fid);
delimiters = [9:13 32]; % White space characters
string_vector=split_string(tline, delimiters);
frewind(fid);
if ( numel(string_vector) > 3 )
if ( iFile == 1 )
fileContainsMetaData =1;
else
if ( fileContainsMetaData == 0 )
error('Not all files seems to have meta data, reading routine not written for that.\n')
end
end
end
if ( fileContainsMetaData )
data = textscan(fid, '%d %d %f %d', 'HeaderLines', 1);
else
data = textscan(fid, '%d %d %f', 'HeaderLines', 1);
end
fclose(fid);
datetmp=data{2};
if ( datetmp(1) < beginDate )
beginDate = datetmp(1);
end
if ( datetmp(end) > endDate )
endDate = datetmp(end);
end
end
% a=0;
% Compute vector with the years, taking leap years into account
beginYearAll = floor(beginDate/1e4);
endYearAll = floor(endDate/1e4);
noValuesYear = [365 366];
noValuesAll = 0;
for iYear = beginYearAll:endYearAll
noValuesAll = noValuesAll + noValuesYear(is_leap(iYear)+1);
% fprintf(1, '%d, %d\n', iYear, noValuesAll)
end
% Compute the dates within the year
[yearAll, julianAll] = generateYearDay(noValuesAll, beginYearAll);
monthAll = zeros(noValuesAll, 1);
dayAll = zeros(noValuesAll, 1);
decYearAll = zeros(noValuesAll, 1);
for iDay = 1:noValuesAll
[day, month] = julian2date(julianAll(iDay), yearAll(iDay));
monthAll(iDay) = month;
dayAll(iDay) = day;
decYearAll(iDay) = calcDecYear(yearAll(iDay), julianAll(iDay));
end
dateAll.day = dayAll;
dateAll.month = monthAll;
dateAll.year = yearAll;
dateAll.julianDay = julianAll;
dateAll.decYear = decYearAll;
% Now the vectors with the years should be ready
% Read data again and insert the data at the right position of the full vector
dataAll = NaN*zeros(noValuesAll, noFiles);
if ( fileContainsMetaData )
metaAll = NaN*zeros(noValuesAll, noFiles);
end
halfADay = 0.5/365.24;
stationNo = zeros(1,noFiles);
for iFile = 1:noFiles
fileName = fullfile(dirName, files(iFile).name);
fid = fopen(fileName);
if ( fileContainsMetaData )
data = textscan(fid, '%d %s %f %d', 'HeaderLines', 1);
else
data = textscan(fid, '%d %s %f', 'HeaderLines', 1);
end
fclose(fid);
stationNoFile = data{1}(1);
stationNo(iFile) = stationNoFile;
datetmp=data{2};
dataFile=data{3};
if ( fileContainsMetaData )
metaFile=data{4};
end
% In network4 there are some values with NA at the end, which are not
% read by Matlab, which is okay as they are the last values. However
% because of this the stationno and date columns are one element longer
if (numel(datetmp) == numel(dataFile)+1 )
datetmp = datetmp(1:numel(dataFile));
end
if (sum(isfinite(dataFile)) <= 365+365+365+366 )
warning('Station no. %d contains less than 4 years of data.\n', iFile);
end
% Convert date to usefull 3 column values
year = zeros(size(datetmp));
month = zeros(size(datetmp));
day = zeros(size(datetmp));
julianDay = zeros(size(datetmp));
decYearFile = zeros(size(datetmp));
noValues = numel(year);
for iDate = 1:noValues
temp = datetmp{iDate};
tempstr = temp(1:4);
year(iDate) = str2num(tempstr); %#ok<*ST2NM>
tempstr = temp(5:6);
month(iDate) = str2num(tempstr); %#ok<*ST2NM>
tempstr = temp(7:8);
day(iDate) = str2num(tempstr); %#ok<*ST2NM>
julianDay(iDate) = date2julian(day(iDate), month(iDate), year(iDate));
decYearFile(iDate) = calcDecYear(year(iDate), julianDay(iDate));
end % for iDate
% Cut away the last part of the data to make the dataset contain only full years of data.
dayCounter = noValues;
lastDay = julianDay(1)-1;
if ( lastDay == 0 )
lastDay = noValuesYear(is_leap(iYear)+1);
end
while ( julianDay(dayCounter) ~= lastDay )
dayCounter = dayCounter - 1;
end
dataFile = dataFile(1:dayCounter);
if ( fileContainsMetaData )
metaFile = metaFile(1:dayCounter);
end
noValues = dayCounter;
% Check if the dates match
index = find(decYearAll-decYearFile(1)+halfADay>0);
index = index(1);
% noValues = numel(datetmp);
for iDate = 1:noValues
if ( (decYearAll(index+iDate-1)-decYearFile(iDate)) ~= 0 )
fprintf(1, '%d, %d\n', decYearAll(index+iDate-1), decYearFile(iDate) )
fprintf(1, '%d\n', iDate)
end
end
dataAll(index:index+noValues-1,iFile) = dataFile;
if ( fileContainsMetaData )
metaAll(index:index+noValues-1,iFile) = metaFile;
end
if 0
figure(1)
imagesc(1:noFiles, decYearAll, dataAll),
colorbar
end
a=0; %#ok<NASGU>
end % for iFile
sumData = sum(dataAll, 2);
if (sum(isfinite(sumData)) <= 365 )
error('There is no period of longer than one year where all stations have data. There may not be overlapping data for every pair. Implement check for that.\n')
end
if ( exist(matlabDataFileDirName, 'file') == 0 )
if ( exist('metaAll', 'var') > 0 )
save(matlabDataFileDirName, 'stationNo', 'dateAll', 'dataAll', 'metaAll');
else
save(matlabDataFileDirName, 'stationNo', 'dateAll', 'dataAll');
end
end
end % if exists matlabDataFileDirName
end % if complete data
a=0; %#ok<NASGU>
dataAll(dataAll<-999) = NaN;
if ( nargout > 3 )
if ( exist('metaAll', 'var') > 0 )
varargout{1} = metaAll;
else
varargout{1} = NaN;
end
end