Add batch processing

atovs_read_find.m

@@ -7,7 +7,18 @@
 % IN    file_name   Name of an MHS data file.
 %
 
-function [hdrinfo, data, err] = atovs_read_headerMHS19_all2(filename)
+function [ret, hdrinfo, data, err] = atovs_read_headerMHS19_all2(filename)
+ret.mhsfile = filename;
+ret.x = -1;
+ret.y = -1;
+ret.maximum = -1;
+ret.lon=-1;
+ret.lat=-1;
+ret.height=-1;
+ret.year = -1;
+ret.day = -1;
+ret.hour = -1;
+ret.minute = -1;
 mean_scannumber = 1000;
 ang_diam = 0.5;
 % we want to read this many bytes (MHS Data Set Header Record size)
@@ -124,6 +135,10 @@ end
 % Das gesuchte Maximum muss mindestens 120 Scans von Anfang und Ende des Erdumlaufs entfernt sein.
 intr_start = 121;
 intr_end = nlines_read-131;
+if intr_end < intr_start
+    disp('Not enough data in file')
+    return
+end
 y1 = double(squeeze(data.dsvcounts11));
 y2 = double(squeeze(data.dsvcounts21));
 y3 = double(squeeze(data.dsvcounts31));
@@ -138,19 +153,29 @@ M(1) = max(A(1,:));
 M(2) = max(A(2,:));
 M(3) = max(A(3,:));
 M(4) = max(A(4,:));
-maximum = max(max(M));
-[x,y]=find(A==maximum);
-longitude=(double(data.lonfov45(y))+double(data.lonfov46(y)))/2E4; % longitude of nadir at maximum of light curve
-latitude=(double(data.latfov45(y))+double(data.latfov46(y)))/2E4; % latitude of nadir at maximum of light curve
-height=double(data.height(y))/10; % altitude of satellite at maximum of light curve
-Z0 = ['Map Ref ',num2str(round(longitude,1)),'deg ',num2str(round(latitude,1)),'deg'];
-disp(Z0)
-Z1 = ['Altitude of Satellite = ',num2str(height),' km'];
-disp(Z1)
-Z2 = ['Maximum counts found in DSV no. ',num2str(x),' for scan ID ',num2str(y)];
-disp(Z2)
-Z3 = ['Date (yyyy-DDD   HR:MN) ',num2str(data.year(y)),'-',num2str(data.day(y)),'   ',num2str(fix(data.time(y)/3600)),':',num2str(fix(data.time(y)/60-60*fix(data.time(y)/3600)))];
-disp(Z3)
+maximum = max(max(M))';
+[ret.x,ret.y]=find(A==maximum);
+ret.maximum = repmat(maximum, size(ret.x));
+% ret.x = ret.x(1);
+% ret.y = ret.y(1);
+ret.lon=((double(data.lonfov45(ret.y))+double(data.lonfov46(ret.y)))/2E4)'; % longitude of nadir at maximum of light curve
+ret.lat=((double(data.latfov45(ret.y))+double(data.latfov46(ret.y)))/2E4)'; % latitude of nadir at maximum of light curve
+ret.height=(double(data.height(ret.y))/10)'; % altitude of satellite at maximum of light curve
+ret.year = data.year(ret.y)';
+ret.day = data.day(ret.y)';
+ret.hour = fix(data.time(ret.y)/3600)';
+ret.minute = fix(data.time(ret.y)/60-60*fix(data.time(ret.y)/3600))';
+ret.mhsfile = repmat(ret.mhsfile, size(ret.x));
+% for i = 1:size(ret.x, 1)
+% Z0 = ['Map Ref ',num2str(round(ret.lon(i),1)),'deg ',num2str(round(ret.lat(i),1)),'deg'];
+% disp(Z0)
+% Z1 = ['Altitude of Satellite = ',num2str(ret.height(i)),' km'];
+% disp(Z1)
+% Z2 = ['Maximum counts found in DSV no. ',num2str(ret.x(i)),' for scan ID ',num2str(ret.y(i))];
+% disp(Z2)
+% Z3 = ['Date (yyyy-DDD   HR:MN) ',num2str(ret.year(1)),'-',num2str(ret.day(1)),'   ',num2str(ret.hour(i)),':',num2str(ret.minute(i))];
+% disp(Z3)
+% end
 end
 %big endian
  function value = extract_uint16(bytearr, startbyte, endbyte);

process_mhs.m

+
+function result = process_mhs(folder, threshold)
+folder = fullfile(folder);
+mhsglob = struct2table(dir([folder '/**/*NSS*.gz']));
+
+% We have two directories per month (name and number) which is why
+% `dir` returns two paths for each file. We need to get rid of the
+% duplicates
+[~, unique_ind] = unique(mhsglob.name);
+mhsfiles = mhsglob.name(unique_ind);
+mhsfolders = mhsglob.folder(unique_ind);
+
+disp(['Found ' num2str(size(mhsfiles, 1)) ' MHS files'])
+
+result = table();
+for i = 1:size(mhsfiles, 1)
+    disp(['Processing ' mhsfiles{i}])
+    mhsfile = [mhsfolders{i} '/' mhsfiles{i}];
+    try
+        r = struct2table(atovs_read_find(mhsfile));
+        result = vertcat(result, r);
+    catch ME
+        warning(['Error processing file ' mhsfile])
+        disp(ME)
+    end
+end
+
+if size(result, 1) == 0
+    warning('Nothing to write')
+    return
+end
+
+parts = strsplit(folder, '/');
+outname = ['output/' strjoin(parts(find(parts=="mhs")+1:end), '-') '-threshold' num2str(threshold)];
+[~, ~] = mkdir('output');
+idx_match = result.maximum >= threshold;
+
+matches = result(idx_match, :);
+misses = result(~idx_match, :);
+if size(matches, 1)
+    filename = [outname '-match.csv'];
+    writetable(matches, filename);
+    disp(['Wrote ' num2str(size(matches, 1)) ' entries to ' filename])
+end
+if size(misses)
+    filename = [outname '-miss.csv'];
+    writetable(misses, filename);
+    disp(['Wrote ' num2str(size(misses, 1)) ' entries to ' filename])
+end