Merge pull request #57 from heikekonow/issue-56-time_correction_for_r…

…adiometer

Issue 56 time correction for radiometer

findRadiometerErrorindex.m

@@ -0,0 +1,29 @@
+function [idx] = findRadiometerErrorindex(time, frequency, varargin)
+
+% Convert time format
+if isdatetime(time)
+    time = datenum(time);
+elseif nargin>2
+    time = datenum(time, varargin{1});
+end
+
+% Check time format
+if ~(all(time>700000) && all(time<800000))
+    error('Time must either be given as serial date number, datetime, or string with string format as third input.')
+end
+
+% Get string for frequency
+freqString = getHAMPfrequencyString(frequency);
+
+% Get date as string
+day = datestr(time(1), 'yyyymmdd');
+
+% Get path to measurement file for current date
+filepath = listFiles([getPathPrefix getCampaignFolder(day) 'radiometer/' freqString '/*' day(3:end) '*BRT*NC'], 'full', 'mat');
+
+% Read time and convert to serial date number
+t = time2001_2sdn(ncread(filepath, 'time'));
+
+% Find indices of given time
+idx(1) = find(t>=time(1), 1, 'first');
+idx(2) = find(t<=time(2), 1, 'last');

functions/cellflat.m

@@ -0,0 +1,59 @@
+function out = cellflat(celllist)
+% CELLFLAT is a helper function to flatten nested cell arrays. 
+% 
+% CELLFLAT(celllist) searches every cell element in cellist and put them on
+% the top most level. Therefore, CELLFLAT linearizes a cell array tree
+% structure. 
+%
+% Example: cellflat({[1 2 3], [4 5 6],{[7 8 9 10],[11 12 13 14 15]},{'abc',{'defg','hijk'},'lmnop'}}) 
+% 
+% Output: 
+%Columns 1 through 7
+%     [1x3 double]    [1x3 double]    [1x4 double]    [1x5 double]    'abc'    'defg'    'hijk'
+%   Column 8 
+%     'lmnop'
+%
+% cellflat(({{1 {2 3}} 'z' {'y' 'x' 'w'} {4 @iscell 5} 6}) )
+% Output: 
+% [1]    [2]    [3]    'z'    'y'    'x'    'w'    [4]    @iscell    [5]    [6]
+%
+% Version: 1.0
+% Author: Yung-Yeh Chang, Ph.D. (yungyeh@hotmail.com)
+% Date: 12/31/2014
+% Copyright 2015, Yung-Yeh Chang, Ph.D.
+% See Also: cell
+%
+% Copyright (c) 2017, Yung-Yeh
+% All rights reserved.
+% 
+% Redistribution and use in source and binary forms, with or without
+% modification, are permitted provided that the following conditions are met:
+% 
+% * Redistributions of source code must retain the above copyright notice, this
+%   list of conditions and the following disclaimer.
+% 
+% * Redistributions in binary form must reproduce the above copyright notice,
+%   this list of conditions and the following disclaimer in the documentation
+%   and/or other materials provided with the distribution
+% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+% FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+if ~iscell(celllist)
+    error('CELLFLAT:ImproperInputAugument','Input argument must be a cell array');
+end
+out = {};
+for idx_c = 1:numel(celllist)
+    if iscell(celllist{idx_c})
+        out = [out cellflat(celllist{idx_c})];
+    else
+        out = [out celllist(idx_c)];
+    end
+end

functions/convertRadarErrorTimes.m

@@ -70,6 +70,11 @@ function convertRadarErrorTimes(campaign)
 
     % Copy errors to variables
     errorsDay = errors{indDay,2};
+
+    % If variable is not cell, convert to cell
+    if ~iscell(errorsDay)
+        errorsDay = {errorsDay};
+    end
 
     % Read time from unified data
     timeUni = ncread(filepathnc,'time');
@@ -80,7 +85,7 @@ function convertRadarErrorTimes(campaign)
     timeErrorFlag = zeros(size(timeUni));
 
     % If errors array is not empty
-    if sum(cellfun(@isempty,errorsDay))==0
+    if sum(cellfun(@isempty,errorsDay))==0 && ~isempty(timeRaw)
 
         % Loop all errors for current date
         for k=1:length(errorsDay)

functions/convertRadiometerErrorTimesSingleChannel.m

@@ -141,4 +141,5 @@ function convertRadiometerErrorTimesSingleChannel(campaign)
 checkandcreate(basefolder, 'aux')
 
 % Save flags to file
-save([basefolder 'aux/errorFlagRadiometer.mat'],'errorFlagSingleChannel','frequencySingleChannel','dateSingleChannel','instrSingleChannel','-append')
+% save([basefolder 'aux/errorFlagRadiometer.mat'],'errorFlagSingleChannel','frequencySingleChannel','dateSingleChannel','instrSingleChannel','-append')
+save([basefolder 'aux/errorFlagRadiometer.mat'],'errorFlagSingleChannel','frequencySingleChannel','dateSingleChannel','-append')

functions/getHAMPfrequencyString.m

@@ -0,0 +1,22 @@
+%% freqString = getHAMPfrequencyString(frequency)
+%
+% Use this function to get the module name string for a given frequency
+%
+% Explanation for different radiometer modules
+% 1: 183,   f>180           (G band)
+% 2: 11990, f>=90 & f<180   (WF band)
+% 3: KV,    f<90            (KV band)
+
+function freqString = getHAMPfrequencyString(frequency)
+
+
+% Translate frequency to string from table above
+if frequency >= 180
+    freqString = '183';
+elseif frequency>=90 && frequency<180
+    freqString = '11990';
+elseif frequency < 90
+    freqString = 'KV';
+else
+    error('Frequency not found')
+end

functions/listFiles.m

@@ -111,7 +111,8 @@
     end
 
     if any(strcmp(varargin, 'mat'))
-        Files = [Files{:}];
+%         Files = [Files{:}];
+        Files = cell2mat(Files);
     end
 end
 

functions/prepareMat2NetCDF.m

@@ -45,6 +45,10 @@
 sizes(strcmp(varnames,'extra_info')) = [];
 % Delete variable name of variable 'extra_info'
 varnames(strcmp(varnames,'extra_info')) = [];
+% Delete size information of variable 'corrCommentString'
+sizes(strcmp(varnames,'corrCommentString')) = [];
+% Delete variable name of variable 'corrCommentString'
+varnames(strcmp(varnames,'corrCommentString')) = [];
 
 % Concatenate dimension sizes and variable size 
 % Here, just the order is changed so that the dimensions will be written

functions/radiometerErrorsSingleChannelLookup.m

@@ -2,6 +2,9 @@
 % Numbers represent indices in raw measurement data and are not yet
 % converted to times
 %
+% Indices can be found by using findRadiometerErrorindex for given time
+% intervals.
+%
 % Note individual indices or intervals in the form [ii jj]
 
 function errors = radiometerErrorsSingleChannelLookup
@@ -15,7 +18,7 @@
           '20200122',    23.84,     {[115703 115760], [119740 119780]};...
           '20200207',    90,     {[94442 94449]};...
           '20200130',    190.81,     {[44438 63007]};...
-          '',    190.81,     {[]};...
+          '20200211',    23.84,     {[22207 22240]};...
           '',    190.81,     {[]};...
           '',    190.81,     {[]};...
     };

functions/radiometerTimeOffset.m

@@ -0,0 +1,118 @@
+function [time, timeOffsetRange, timeOffsetValue, commentstr] = ...
+                            radiometerTimeOffset(flightdate, frequency, time)
+
+timeOffsetsAll = {...
+    '20200119',    'HAMP-WF',    ':',          -141;
+    '20200122',    'none',       ':',             0;
+    '20200124',    'HAMP-KV',    ':',            -1;
+    '20200124',    'HAMP-WF',    '1:20501',      -2;
+    '20200124',    'HAMP-G',     '1:21001',      -7;
+    '20200124',    'HAMP-G',     '21001:29401',  -3;
+    '20200126',    'HAMP-KV',    '1:23351',      -2;
+    '20200126',    'HAMP-WF',    '1:23351',       1;
+    '20200126',    'HAMP-G',     '23701:',       -3;
+    '20200128',    'none',       ':',             0;
+    '20200130',    'none',       ':',             0;
+    '20200131',    'HAMP-WF',    ':',            -1;
+    '20200131',    'HAMP-G',     '21201:',       -3;
+    '20200202',    'HAMP-KV',    '4001:6211',    -2;
+    '20200202',    'HAMP-WF',    '1:2001',       -3;
+    '20200202',    'HAMP-WF',    '2001:4001',    -1;
+    '20200202',    'HAMP-WF',    '4001:6255',    -3;
+    '20200202',    'HAMP-WF',    '6255:6721',    -2;
+    '20200202',    'HAMP-WF',    '6721:',        -1;
+    '20200202',    'HAMP-G',     ':',            -2;
+    '20200205',    'HAMP-KV',    ':',            -3;
+    '20200205',    'HAMP-WF',    ':',            -4;
+    '20200207',    'HAMP-KV',    '21301:23261',  -2;
+    '20200207',    'HAMP-WF',    '1:22401',       2;
+    '20200209',    'HAMP-G',     ':',            -2;
+    '20200211',    'HAMP-KV',    '1:22901',      -2;
+    '20200211',    'HAMP-G',     '21001:23501',  -5;
+    '20200213',    'none',       ':',             0;
+    '20200215',    'HAMP-KV',    '20484:',       -1;
+    '20200215',    'HAMP-WF',    '2034:28757',    2;
+    '20200218',    'HAMP-WF',    ':',            -1;
+    '20200218',    'HAMP-G',     ':',            -4;
+    };
+
+% Copy data to variables
+offsetDates = timeOffsetsAll(:,1);
+offsetModules = timeOffsetsAll(:,2);
+
+% Find date index
+indexDate = strcmp(flightdate, offsetDates);
+
+% Explanation for different radiometer modules
+% 1: 183,   f>180           (G band)
+% 2: 11990, f>=90 & f<180   (WF band)
+% 3: KV,    f<90            (KV band)
+
+% Translate frequency to string from table above
+if frequency >= 180
+    freqString = 'HAMP-G';
+elseif frequency>=90 && frequency<180
+    freqString = 'HAMP-WF';
+elseif frequency < 90
+    freqString = 'HAMP-KV';
+else
+    error('Frequency not found')
+end
+
+% Find frequency index
+indexFrequency = strcmp(freqString, offsetModules);
+
+% Find row where frequency and data match given date and frequency
+indexEntry = indexFrequency & indexDate;
+
+if any(indexEntry)
+
+    % List Modules in variable
+    modules = offsetModules(indexEntry);
+
+    % Extract entries from table
+    timeOffsetRange = timeOffsetsAll(indexEntry, 3);
+    timeOffsetValue = cell2mat(timeOffsetsAll(indexEntry, 4));
+
+    %%%%%%%%%%%
+    % Apply time offset values
+
+    commentstr = cell(1, length(timeOffsetRange));
+    for i=1:length(timeOffsetRange)
+
+        % Get colon position from string
+        colPos = regexp(timeOffsetRange{i}, ':');
+
+        % Analyse time offset index
+        if strcmp(timeOffsetRange{i}, ':')          % ':'
+            ind(1) = 1;
+            ind(2) = length(time);
+
+        elseif strncmp(timeOffsetRange{i}, ':', 1)  % ':yyyy'
+            ind(1) = 1;
+            a = timeOffsetRange{i}(2:end);
+            ind(2) = str2double(a);
+
+        elseif colPos==length(timeOffsetRange{i})
+            a = timeOffsetRange{i}(1:colPos-1);     % 'xxxx:'
+
+            ind(1) = str2double(a);
+            ind(2) = length(time);
+
+        else                                        % 'xxxx:yyyy'
+            ind{1} = timeOffsetRange{i}(1:colPos-1);
+            ind{2} = timeOffsetRange{i}(colPos+1:end);
+            ind = cellfun(@str2double, ind);
+        end
+
+        % Apply offset to time array
+        time(ind(1):ind(2)) = time(ind(1):ind(2)) + timeOffsetValue(i) ./24./60./60;
+        clear ind
+
+        commentstr{i} = [modules{i} ': (' timeOffsetRange{i} ') ' num2str(timeOffsetValue(i)) ' s'];
+    end
+else
+    commentstr = {};
+    timeOffsetRange = {};
+    timeOffsetValue = {};
+end

functions/removeTBErrors.m

@@ -10,7 +10,7 @@
 
 if any(strcmp(extractfield(vars, 'name'), 'errorFlagSingleChannel'))
 
-    load(datafile, 'dateSingleChannel', 'errorFlagSingleChannel', 'frequencySingleChannel', 'instrSingleChannel')
+    load(datafile, 'dateSingleChannel', 'errorFlagSingleChannel', 'frequencySingleChannel')
 
     singleChannelError = true;
 end

functions/run_unifyGrid.m

@@ -1,5 +1,6 @@
 function run_unifyGrid(version, subversion, flightdates_use, comment, contact, altitudeThreshold, ...
-                        rollThreshold, radarmask,  radarClutter, missingvalue, fillvalue, filenameprefix)
+                        rollThreshold, radarmask,  radarClutter, correctRadiometerTime, ...
+                        missingvalue, fillvalue, filenameprefix)
 
 tic 
 %% Switches 
@@ -71,7 +72,7 @@ function run_unifyGrid(version, subversion, flightdates_use, comment, contact, a
 
         % Radiometer
         unifyGrid_radiometer(pathtofolder,flightdates_use{i},uniTime,radiometerVars,...
-            altitudeThreshold, rollThreshold, missingvalue, fillvalue)
+            altitudeThreshold, rollThreshold, missingvalue, fillvalue, correctRadiometerTime)
 
         % Radar
         unifyGrid_radar(pathtofolder,flightdates_use{i},uniHeight,uniTime,radarVars)
@@ -94,7 +95,7 @@ function run_unifyGrid(version, subversion, flightdates_use, comment, contact, a
 
 %% Export to netcdf
 
-instr = {'radar','bahamas','radiometer','dropsondes'};
+instr = {'radiometer','radar','bahamas','dropsondes'};
 % instr = {'bahamas','radar','radiometer'};
 % instr = {'radar'};
 % instr = {'bahamas'};
@@ -219,6 +220,11 @@ function run_unifyGrid(version, subversion, flightdates_use, comment, contact, a
                     removeClutter(outfile, missingvalue, fillvalue)
                 end 
 
+                %% Add comment about radiometer time correction
+                if correctRadiometerTime && strcmp(instr{j}, 'radiometer')
+                    addRadiometerTimeComment(outfile, infile)
+                end
+
             else
                 disp(['No ' instr{j} ' data found'])
             end
@@ -429,3 +435,9 @@ function removeSideLobes(outfile, rollThreshold, fillvalue, radarmask)
         ncwriteatt(outfile, 'data_flag', 'long_name', longNameAtt)
     end
 end
+
+function addRadiometerTimeComment(outfile, infile)
+    string = load(infile, 'corrCommentString');
+
+    ncwriteatt(outfile, 'tb', 'comment', string.corrCommentString)
+end