Home > atmlab > deprecated > collocate_granule.m

collocate_granule

PURPOSE ^

collocate_granule Collocate granules from sat/sensor pairs

SYNOPSIS ^

function [overlap, data, meandata] =collocate_granule(sat1, sensor1, sat2, sensor2, date1, oneday)

DESCRIPTION ^

 collocate_granule Collocate granules from sat/sensor pairs

 Given the names of the satellite and the sensors, as well as the starting
 date/time for the primary satellite/sensor-pair, this will find all
 collocations between the two. It finds the granule for the primary
 sat/sensor-pair, all the overlapping granules for the secondary
 sat/sensor-pair, and collocates with each of those granules. It then
 collects additional information and, if applicable, calculates mean
 values for each of the collocations (for example, for CloudSat/MHS, it
 return a matrix with one row for each CloudSat pixel, as well as a matrix
 with one row for each MHS pixel).

 A special value for 'sat2' can be 'poes', which is short for 'every POES
 satellite'.

 The argument 'oneday' specifies whether data should be limited to the
 starting date for sat1 (oneday=1), the day after this (oneday=2), or no
 limitation (oneday=0, default).

 FORMAT

   [overlap, data, meandata] = ...
       collocate_granule(sat1, sensor1, sat2, sensor2, date1, oneday)

 IN

   sat1        string      Primary satellite (example: 'CloudSat')
   sensor1     string      Primary sensor (example: 'cpr')
   sat2        string      Secondary satellite (example: 'noaa18')
   sensor2     string      Secondary sensor (example: 'mhs')
   date1       vector      Date vector for starting time sat1/sensor1
   oneday      number      Limit data to one day.
                           1=start date for sat1, 2=day after

 OUT

   overlap     structure   With fields for each satellite (normally just 
                           one, but for 'poes' multiple), collocation
                           information along with times, distances, etc.
   data        structure   With fields for each satellite, data
                           corresponding to the collocation information,
                           IWP, brightness temperatures, etc.
   meandata    structure   With fields for each satellite, for instruments
                           with very different footprint sizes, statistics
                           for the smaller within the larger: number of
                           points, mean, standard deviation, etc.
                           Otherwise empty.

 $Id$

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

collocate_granule.m

SOURCE CODE ^

0001 function [overlap, data, meandata] = ...
0002     collocate_granule(sat1, sensor1, sat2, sensor2, date1, oneday)
0003 
0004 % collocate_granule Collocate granules from sat/sensor pairs
0005 %
0006 % Given the names of the satellite and the sensors, as well as the starting
0007 % date/time for the primary satellite/sensor-pair, this will find all
0008 % collocations between the two. It finds the granule for the primary
0009 % sat/sensor-pair, all the overlapping granules for the secondary
0010 % sat/sensor-pair, and collocates with each of those granules. It then
0011 % collects additional information and, if applicable, calculates mean
0012 % values for each of the collocations (for example, for CloudSat/MHS, it
0013 % return a matrix with one row for each CloudSat pixel, as well as a matrix
0014 % with one row for each MHS pixel).
0015 %
0016 % A special value for 'sat2' can be 'poes', which is short for 'every POES
0017 % satellite'.
0018 %
0019 % The argument 'oneday' specifies whether data should be limited to the
0020 % starting date for sat1 (oneday=1), the day after this (oneday=2), or no
0021 % limitation (oneday=0, default).
0022 %
0023 % FORMAT
0024 %
0025 %   [overlap, data, meandata] = ...
0026 %       collocate_granule(sat1, sensor1, sat2, sensor2, date1, oneday)
0027 %
0028 % IN
0029 %
0030 %   sat1        string      Primary satellite (example: 'CloudSat')
0031 %   sensor1     string      Primary sensor (example: 'cpr')
0032 %   sat2        string      Secondary satellite (example: 'noaa18')
0033 %   sensor2     string      Secondary sensor (example: 'mhs')
0034 %   date1       vector      Date vector for starting time sat1/sensor1
0035 %   oneday      number      Limit data to one day.
0036 %                           1=start date for sat1, 2=day after
0037 %
0038 % OUT
0039 %
0040 %   overlap     structure   With fields for each satellite (normally just
0041 %                           one, but for 'poes' multiple), collocation
0042 %                           information along with times, distances, etc.
0043 %   data        structure   With fields for each satellite, data
0044 %                           corresponding to the collocation information,
0045 %                           IWP, brightness temperatures, etc.
0046 %   meandata    structure   With fields for each satellite, for instruments
0047 %                           with very different footprint sizes, statistics
0048 %                           for the smaller within the larger: number of
0049 %                           points, mean, standard deviation, etc.
0050 %                           Otherwise empty.
0051 %
0052 % $Id$
0053 
0054 % FIXME: update to new-style
0055 warning(['atmlab:' mfilename], 'old style function, being phased out, use OO way');
0056 
0057 if ~exist('oneday', 'var')
0058     oneday = 0;
0059 end
0060 
0061 %% output
0062 
0063 fid = atmlab('OUT');
0064 eid = atmlab('ERR');
0065 
0066 
0067 %% get handles for reading, processing
0068 
0069 dummy = @(varargin) (-1);
0070 fhp = colloc_constants(['process_' sensor1 '_' sensor2]);
0071 
0072 try
0073     fhpd = colloc_constants(['process_data_', sensor1 '_' sensor2]);
0074     fhpm = colloc_constants(['process_meandata_', sensor1 '_' sensor2]);
0075 catch ME
0076     switch (ME.identifier)
0077         case 'MATLAB:nonExistentField'
0078             if any(strfind(ME.message, '_meandata_'))
0079                 fhpm = dummy;
0080             else % failed already at _data_
0081                 fhpd = dummy;
0082                 fhpm = dummy;
0083             end
0084         otherwise
0085             ME.rethrow();
0086     end
0087 end
0088 
0089 
0090 %% set 'date2_all', 'overlap', and 'sat2_all_names'
0091 % find date/times for which to search for secondary granules
0092 % Also pass on the primary satellite, because for poes/poes, we don't want
0093 % overlaps where sat1==sat2
0094 
0095 % overlap check 1: primary filename, secondary filename
0096 [date2_all, sat2_all_names] = overlap_granule(sat1, sensor1, ...
0097     date1, sat2, sensor2);
0098 if size(date2_all, 1)==0
0099     error('atmlab:collocate_granule:noother', ...
0100         'No overlapping granules found for [%s] %s/%s with %s/%s', ...
0101         num2str(date1), sat1, sensor1, sat2, sensor2);
0102 end
0103 
0104 % pre-allocate overlap-structure
0105 % TODO: rather than [], have good width already
0106 uniq_names = unique(sat2_all_names);
0107 overlap = cell2struct(repmat({[]}, [length(uniq_names) 1]), uniq_names);
0108 
0109 data = overlap; % same pre-allocation
0110 meandata = overlap; % this one as well
0111 M_c = overlap; M_d = overlap; M_m = overlap; % and those
0112 
0113 %% read data for sat1/sensor1
0114 % Here, I decide to keep the doubles for now. Simple is better than
0115 % complex. My collocation output are rows and columns, and if I already
0116 % remove doubles /before/ collocation, the row-numbers will be incorrect or
0117 % at least need to be corrected for. This causes more problems than it
0118 % solves, so I will remove the doubles in the postprocessing instead.
0119 try
0120     data1 = read_granule(sat1, sensor1, date1, false); % keep doubles for now
0121 catch ME
0122     switch ME.identifier
0123         case {'MATLAB:load:couldNotReadFile', 'atmlab:invalid_data', 'MATLAB:imagesci:hdfinfo:fileOpen'}
0124             logtext(eid, 'Collocating failed: %s\n', ME.message);
0125             return
0126         otherwise
0127             ME.rethrow();
0128     end
0129 end
0130 
0131 overlap.version{1} = data1.version;
0132 overlap.version{2} = '?';
0133 
0134 if isempty(data1.time)
0135     logtext(fid, 'Primary is empty, nothing to collocate\n');
0136     overlap.version{2} = 'N/A';
0137     return
0138 end
0139 
0140 % keep effective range, check time overlap before actually reading the
0141 % secondary data, because the primary might be sparse. E.g. when judging
0142 % from the filenames the granules contain overlapping periods, but if the
0143 % primary in reality only contains data for a certain part of this period
0144 % (e.g. in the case of collocations), should check the primary data with
0145 % the secondary filename in order not to read files one knows will not
0146 % contain overlap anyway.
0147 data1.time_orig = data1.time; % store because of subsequent unify_time_axis can lead to errors
0148 data1.eff_range = double(data1.epoch) + [min(data1.time)-colloc_config('interval'), max(data1.time)+colloc_config('interval')];
0149 
0150 logtext(fid, 'Found %d other granules to collocate with\n', size(date2_all, 1));
0151 
0152 
0153 %% loop through all the granules to collocate with
0154 
0155 for i = 1:size(date2_all, 1)
0156     datecell = num2cell(date2_all(i, 1:5));
0157     satname = sat2_all_names{i};
0158     
0159     %% TODO/FIXME: GET RID OF UGLY SPECIAL CASE HACK!
0160     if strcmp(sat2, 'poes') && any(strcmp(sensor2, {'amsub', 'mhs'}))
0161         sens = amsub_or_mhs(satname);
0162     else
0163         sens = sensor2;
0164     end
0165     %% END OF UGLY SPECIAL CASE HACK!
0166     %% Continue sensibly.
0167     
0168     logtext(atmlab('OUT'), 'Collocating with %s %s %04d-%02d-%02d %02d:%02d\n', ...
0169         satname, sens, datecell{1:5});
0170     
0171     % Overlap check 2; primary data with secondary filename
0172     data2_start = date2unixsecs(datecell{1:5});
0173     data2_end = data2_start + datasets_constants(['granule_duration_' sensor2]);
0174     
0175     if (data2_start > data1.eff_range(2)) || (data2_end < data1.eff_range(1))
0176         % no possible time overlap; data1 range already compensated for
0177         % colloc_config('interval'); no use even reading granule2
0178         logtext(atmlab('OUT'), 'Not reading, no overlap with primary\n');
0179         continue
0180     end
0181 
0182     try
0183         data2 = read_granule(satname, sens, date2_all(i, 1:5), false); % keep doubles for now
0184     catch ME
0185         switch ME.identifier
0186             case {'atmlab:find_datafile_by_date', 'atmlab:invalid_data'}
0187                 logtext(eid, 'Error in reading datafile %4d-%02d-%02d %02d:%02d: %s. SKIPPING\n', ...
0188                     datecell{:}, ME.message);
0189                 continue
0190             otherwise
0191                 ME.rethrow();
0192         end
0193     end
0194     
0195     overlap.version{2} = data2.version;
0196     data2.time_orig = data2.time; % store, see note at data1.time_orig
0197     if isempty(data2.time)
0198         logtext(fid, 'Secondary is empty, skipping\n');
0199         continue
0200     end
0201     
0202     switch sign(data1.epoch - data2.epoch)
0203         case -1 % data1 epoch is earliest
0204             data2.time = data2.time + double(data2.epoch - data1.epoch);
0205             data2.epoch = data1.epoch;
0206         case 0 % do nothing
0207         case 1 % data2 epoch is earliest
0208             data1.time = data1.time + double(data1.epoch - data2.epoch);
0209             data1.epoch = data2.epoch;
0210         otherwise
0211             error('atmlab:collocate_granule', 'Reached impossible place. Bug.');
0212     end
0213 
0214     % overlap check 3: primary data, secondary data
0215     [iv1, iv2] = find_common_time(data1.time, data2.time, colloc_config('interval'));
0216     % limit to one day (if needed)
0217     switch oneday
0218         case 1
0219             % only consider data on the same date as the start for data1
0220             if data1.time(1) < 86400
0221                 sameday = data1.time < 86400;
0222             else
0223                 sameday = data1.time < 2*86400;
0224             end
0225             if ~all(sameday)
0226                 logtext(fid, 'taking only %d/%d that are on the same day\n', ...
0227                     sum(sameday), length(data1.time));
0228             end
0229             iv1 = iv1 & sameday;
0230         case 2
0231             % only consider data on date after the start date for data1
0232             sameday = data1.time > 86400;
0233             iv1 = iv1 & sameday;            
0234             if ~all(sameday)
0235                 logtext(fid, 'taking only %d/%d that are on the next day\n', ...
0236                     sum(sameday), length(data1.time));
0237             end
0238     end
0239     if ~(any(iv1) && any(iv2)) % no overlap
0240         logtext(fid, 'no time overlap, so nothing to collocate\n');
0241         continue
0242     end
0243     
0244     % perform collocations
0245     collocations = collocate(data1.time(iv1), data1.lat(iv1, :), data1.lon(iv1, :), ...
0246         data2.time(iv2), data2.lat(iv2, :), data2.lon(iv2, :), ...
0247         colloc_config('distance'), colloc_config('interval'));
0248     if any(collocations)
0249         % compensate for the fact that we passed only a subset to collocate
0250         collocations(:, 1) = collocations(:, 1) + find(iv1, 1, 'first') - 1;
0251         collocations(:, 3) = collocations(:, 3) + find(iv2, 1, 'first') - 1;
0252     else
0253         logtext(fid, 'No collocations\n');
0254         continue
0255     end
0256     % process data
0257     % should return a matrix with info
0258     logtext(fid, 'Collecting info for %d collocations\n', size(collocations, 1));
0259     M_c.(satname) = fhp(collocations, ...
0260         sat1, date1, data1, ...
0261         satname, date2_all(i, :), data2);
0262 
0263     if nargout > 1 && ~isequal(fhpd, dummy) % also data
0264         M_d.(satname) = fhpd(collocations, M_c.(satname), ...
0265             sat1, date1, data1, ...
0266             satname, date2_all(i, :), data2);
0267     end
0268     
0269     if nargout > 2 && ~isequal(fhpm, dummy) % also meandata
0270         M_m.(satname) = fhpm(collocations, M_c.(satname), M_d.(satname), date1, data1, date2_all(i, :), data2);
0271     end
0272     
0273     [overlap, data, meandata] = colloc_concatenate_colloc_data_mean(...
0274         overlap, M_c, data, M_d, meandata, M_m);
0275     
0276     logtext(fid, 'Info collected\n');
0277 end
0278 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005