Home > atmlab > deprecated > colloc_process_cpr_mhs.m

colloc_process_cpr_mhs

PURPOSE ^

colloc_process_cpr_mhs Gather collocation info to store for CloudSat+MHS

SYNOPSIS ^

function M = colloc_process_cpr_mhs(collocations,~, date_cpr, data_cpr, name_mhs, date_mhs, data_mhs)

DESCRIPTION ^

 colloc_process_cpr_mhs Gather collocation info to store for CloudSat+MHS

 After the collocating process is done by 'collocate', one wants to
 collect some additional information: longitudes, latitudes, times,
 viewing angles, distances, intervals, etc. This function takes care of
 that in the case of CPR/MHS.

 All colloc_process_* functions are automatically called by
 collocate_granule and thus have the same interface. Not all arguments are
 used by all.
 
 FORMAT

   M = colloc_process_cpr_mhs(collocations, ...
       name_cpr, date_cpr, data_cpr, ...
       name_mhs, date_mhs, data_mhs)

 IN

   collocations    Nx4 matrix  Collocations as returned by collocate
   name_cpr        string      Name of satellite carrying CPR ('cloudsat')
   date_cpr        vector      Date vector for CPR starting time
   data_cpr        structure   Data for CPR, as returned by read_cpr
   name_mhs        string      Name of MHS-carrying satellite
   date_mhs        vector      Date vector for MHS starting time
   data_mhs        structure   Data for MHS, as returned by read_mhs

 OUT

   M               matrix      Matrix with one row for each collocation
                               and columns as in
                               colloc_constants('cols_cpr_mhs').overlap

 $Id$

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

colloc_process_cpr_mhs.m

SOURCE CODE ^

0001 function M = colloc_process_cpr_mhs(collocations, ...
0002     ~, date_cpr, data_cpr, name_mhs, date_mhs, data_mhs)
0003 
0004 % colloc_process_cpr_mhs Gather collocation info to store for CloudSat+MHS
0005 %
0006 % After the collocating process is done by 'collocate', one wants to
0007 % collect some additional information: longitudes, latitudes, times,
0008 % viewing angles, distances, intervals, etc. This function takes care of
0009 % that in the case of CPR/MHS.
0010 %
0011 % All colloc_process_* functions are automatically called by
0012 % collocate_granule and thus have the same interface. Not all arguments are
0013 % used by all.
0014 %
0015 % FORMAT
0016 %
0017 %   M = colloc_process_cpr_mhs(collocations, ...
0018 %       name_cpr, date_cpr, data_cpr, ...
0019 %       name_mhs, date_mhs, data_mhs)
0020 %
0021 % IN
0022 %
0023 %   collocations    Nx4 matrix  Collocations as returned by collocate
0024 %   name_cpr        string      Name of satellite carrying CPR ('cloudsat')
0025 %   date_cpr        vector      Date vector for CPR starting time
0026 %   data_cpr        structure   Data for CPR, as returned by read_cpr
0027 %   name_mhs        string      Name of MHS-carrying satellite
0028 %   date_mhs        vector      Date vector for MHS starting time
0029 %   data_mhs        structure   Data for MHS, as returned by read_mhs
0030 %
0031 % OUT
0032 %
0033 %   M               matrix      Matrix with one row for each collocation
0034 %                               and columns as in
0035 %                               colloc_constants('cols_cpr_mhs').overlap
0036 %
0037 % $Id$
0038 
0039 % prepare
0040 c = colloc_constants('cols_cpr_mhs');
0041 c = c.overlap;
0042 n = size(collocations, 1);
0043 M = nan*zeros(n, c.NCOLS);
0044 if n==0
0045     return % don't bother
0046 end
0047 
0048 % convert to cell array for easy passing into date2unixsecs
0049 date_cpr_cell = num2cell(date_cpr);
0050 date_mhs_cell = num2cell(date_mhs);
0051 
0052 % row and column numbers
0053 cr = collocations(:, 1);
0054 ar = collocations(:, 3);
0055 ac = collocations(:, 4);
0056 
0057 % index for direct addressing
0058 ai = sub2ind(size(data_mhs.lat), ar, ac);
0059 
0060 %% cloudsat info
0061 
0062 % lat/long
0063 M(:, c.C_LONG) = data_cpr.lon(cr);
0064 M(:, c.C_LAT) = data_cpr.lat(cr);
0065 
0066 % time
0067         
0068 M(:, c.C_START, :) = round(date2unixsecs(date_cpr_cell{1:5}));
0069 M(:, c.C_TIME, :) = data_cpr.epoch + data_cpr.time(cr);
0070 M(:, c.C_I, :) = cr;
0071 
0072 %% amsub/mhs
0073 
0074 M(:, c.B_LONG) = data_mhs.lon(ai);
0075 M(:, c.B_LAT) = data_mhs.lat(ai);
0076 M(:, c.B_START) = round(date2unixsecs(date_mhs_cell{1:5}));
0077 M(:, c.B_TIME) = data_mhs.epoch + data_mhs.time(ar);
0078 M(:, c.B_I) = ar;
0079 M(:, c.B_C) = ac;
0080 M(:, c.B_DIST) = sphdist(data_cpr.lat(cr), data_cpr.lon(cr), ...
0081     data_mhs.lat(ai), data_mhs.lon(ai), constants('EARTH_RADIUS')/1e3);
0082 M(:, c.B_INT) = M(:, c.B_TIME) - M(:, c.C_TIME);
0083 M(:, c.B_LZA) = data_mhs.lza(ai);
0084 M(:, c.B_LAA) = data_mhs.laa(ai);
0085 
0086 %% amsua
0087 
0088 amsua = read_granule(name_mhs, 'amsua', date_mhs, false);
0089 
0090 [aah aaw] = size(amsua.lon); % amsu-a width, amsu-a height
0091 M(:, c.A_I) = min(max(round(ar/3), 1), aah); % amsu-a row
0092 M(:, c.A_C) = min(max(round(ac/3), 1), aaw); % amsu-a column
0093 aai = sub2ind(size(amsua.lat), M(:, c.A_I), M(:, c.A_C));
0094 M(:, c.A_LONG) = amsua.lon(aai);
0095 M(:, c.A_LAT) = amsua.lat(aai);
0096 M(:, c.A_START) = M(:, c.B_START);
0097 M(:, c.A_TIME) = data_mhs.epoch + ...
0098     amsua.time(M(:, c.A_I));
0099 M(:, c.A_DIST) = sphdist(data_cpr.lat(cr), data_cpr.lon(cr), ...
0100     amsua.lat(aai), amsua.lon(aai), constants('EARTH_RADIUS')/1e3);
0101 M(:, c.A_INT) = M(:, c.A_TIME) - M(:, c.C_TIME);
0102 
0103 %% hirs
0104 
0105 hirs = read_granule(name_mhs, 'hirs', date_mhs, false, true);
0106 
0107 
0108 if ~isempty(hirs)
0109     % TODO: optimise this bit of code
0110     % find hirs scanline (row) closest in time for each mhs scanline
0111     hr = arrayfun(...
0112         @(v) mini(abs(hirs.time-v)), ...
0113         data_mhs.time(M(:, c.B_I)));
0114     % find hirs pixel (col) in scanline closest in distance to CPR-point
0115     hc = arrayfun(...
0116         @(i) mini(...
0117         (hirs.lat(hr(i), :)-M(i, c.C_LAT)).^2 + ...
0118         (hirs.lon(hr(i), :)-M(i, c.C_LONG)).^2 ...
0119         ), 1:n);
0120     % find index from (row, column)
0121     hi = sub2ind(size(hirs.lat), hr, hc');
0122     % ready to add information to matrix
0123     M(:, c.H_LONG) = hirs.lon(hi);
0124     M(:, c.H_LAT) = hirs.lat(hi);
0125     M(:, c.H_START) = M(:, c.B_START);
0126     M(:, c.H_TIME) = data_mhs.epoch + hirs.time(hr);
0127     M(:, c.H_I) = hr;
0128     M(:, c.H_C) = hc;
0129     M(:, c.H_DIST) = sphdist(M(:, c.C_LAT), M(:, c.C_LONG), ...
0130         M(:, c.H_LAT), M(:, c.H_LONG), constants('EARTH_RADIUS')/1e3);
0131     M(:, c.H_INT) = M(:, c.H_TIME) - M(:, c.C_TIME);
0132 end
0133 
0134 %% others
0135 
0136 %M(:, c.SZA) = circ_rad2ang(ignoreNaN([data_mhs.sza(ai) amsua.sza(aai) hirs.sza(hi)], ...
0137 %    @circ_meand, 2));
0138 %M(:, c.SAA) = circ_rad2ang(ignoreNaN([data_mhs.saa(ai) amsua.saa(aai) hirs.saa(hi)], ...
0139 %    @circ_meand, 2));
0140 % assert(all(circ_distd(data_mhs.sza(ai), amsua.sza(aai))<0.5), ...
0141 %     'atmlab:colloc_process_cpr_mhs', ...
0142 %     'solar zenith angles MHS/AMSUA do not match. Bug?');
0143 
0144 %% sort by: MHS row, then MHS column, then CloudSat pixel
0145 
0146 M = sortrows(M, [c.B_I c.B_C c.C_I]);
0147 
0148 %% remove doubles
0149 
0150 try
0151     wrongrows = M(:, c.B_I)< granule_first_line(name_mhs, amsub_or_mhs(name_mhs), date_mhs);
0152     M(wrongrows, :) = [];
0153 catch ME
0154     switch ME.identifier
0155         case 'atmlab:granule_first_line'
0156             logtext(atmlab('ERR'), ...
0157                 'Warning: Unable to remove doubles: %s\n', ...
0158                 ME.message);
0159         otherwise
0160             ME.rethrow();
0161     end
0162 end
0163 
0164 end
0165 
0166 function i = mini(v)
0167 % mini return index for minimum of v (arrayfun wants 1 output)
0168 [~, i] = min(v);
0169 end

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