colloc_process_data_cpr_mhs Gather collocation data to store for CPR+MHS
After collocating, and gathering information about the collocations, we
want to collect some actual data: in this case, brightness temperatures.
Because this m-file has the same format as other colloc_process_data_* files
(so that the caller doesn't need to know what is being collocated), some of
the arguments are ignored.
FORMAT
M = process_data_cpr_mhs(collocations, M_c, sat1, date1, data1, ...
sat2, date2, data2)
IN
collocations Nx4 matrix (ignored) as returned by collocate
M_c Nxp matrix As returned by process_cpr_mhs
sat1 string name of 1st satellite
date1 1x5 vector (ignored) starting time of 1st satellite
data1 structure As returned by read_cpr
sat2 string name of 2nd satellite
date2 1x5 vector As date1, but for secondary granule.
data2 structure As returned by read_mhs
OUT
M Nxp matrix Matrix containing p fields of information for
all N collocations.
$Id$0001 function M = colloc_process_data_cpr_mhs(~, M_c, ~, ~, data1, sat2, date2, data2) 0002 0003 % colloc_process_data_cpr_mhs Gather collocation data to store for CPR+MHS 0004 % 0005 % After collocating, and gathering information about the collocations, we 0006 % want to collect some actual data: in this case, brightness temperatures. 0007 % 0008 % Because this m-file has the same format as other colloc_process_data_* files 0009 % (so that the caller doesn't need to know what is being collocated), some of 0010 % the arguments are ignored. 0011 % 0012 % FORMAT 0013 % 0014 % M = process_data_cpr_mhs(collocations, M_c, sat1, date1, data1, ... 0015 % sat2, date2, data2) 0016 % 0017 % IN 0018 % 0019 % collocations Nx4 matrix (ignored) as returned by collocate 0020 % M_c Nxp matrix As returned by process_cpr_mhs 0021 % sat1 string name of 1st satellite 0022 % date1 1x5 vector (ignored) starting time of 1st satellite 0023 % data1 structure As returned by read_cpr 0024 % sat2 string name of 2nd satellite 0025 % date2 1x5 vector As date1, but for secondary granule. 0026 % data2 structure As returned by read_mhs 0027 % 0028 % OUT 0029 % 0030 % M Nxp matrix Matrix containing p fields of information for 0031 % all N collocations. 0032 % 0033 % $Id$ 0034 0035 % prepare 0036 c = colloc_constants('cols_cpr_mhs'); 0037 n = size(M_c, 1); 0038 M = nan*zeros(n, c.data.NCOLS); 0039 0040 % row and column numbers 0041 r1 = M_c(:, c.overlap.C_I); 0042 r2 = M_c(:, c.overlap.B_I); 0043 c2 = M_c(:, c.overlap.B_C); 0044 0045 %% cloudsat data 0046 M(:, c.data.ROIWP) = data1.RO_ice_water_path(r1); 0047 M(:, c.data.dROIWP) = data1.RO_ice_water_path_uncertainty(r1); 0048 M(:, c.data.IOROIWP) = data1.IO_RO_ice_water_path(r1); 0049 M(:, c.data.dIOROIWP) = data1.IO_RO_ice_water_path_uncertainty(r1); 0050 0051 %% AMSUB/MHS data 0052 0053 % index for direct addressing 0054 i2 = sub2ind(size(data2.lat), r2, c2); 0055 0056 % reshape so that I can use direct addressing for brightness temperatures 0057 tb = reshape(data2.tb, [numel(data2.lat) 5]); 0058 0059 M(:, c.data.MHS) = tb(i2, :); 0060 0061 %% AMSUA, HIRS 0062 % those require reading additional files 0063 0064 % TODO/FIXME: optimise this, I already read for the earlier processing 0065 amsua_data = read_granule(sat2, 'amsua', date2, false, true); 0066 hirs_data = read_granule(sat2, 'hirs', date2, false, true); 0067 0068 if ~isempty(amsua_data) 0069 % AMSU_A 0070 ai = sub2ind(size(amsua_data.lat), M_c(:, c.overlap.A_I), M_c(:, c.overlap.A_C)); 0071 tb = reshape(amsua_data.tb, [numel(amsua_data.lat) 15]); 0072 M(:, c.data.AMSU_A) = tb(ai, :); 0073 end 0074 0075 if ~isempty(hirs_data) 0076 % HIRS 0077 hi = sub2ind(size(hirs_data.lat), M_c(:, c.overlap.H_I), M_c(:, c.overlap.H_C)); 0078 tb = reshape(hirs_data.tb, [numel(hirs_data.lat) 20]); 0079 M(:, c.data.HIRS) = tb(hi, :); 0080 end 0081 0082 %% MSPPS 0083 0084 % Apply scalefactor. 0085 % From 0086 % http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c9/sec96-3.htm 0087 % divide by the scalefactor to get kg/m^2 0088 % we want g/m^2 so we multiply by 1000/scalefactor 0089 % TODO: convert to specialised reading function 0090 try 0091 mspps_file = find_datafile_by_date(date2, sat2, 'mspps', 120); % 120 seconds tolerance 0092 logtext(atmlab('OUT'), ... 0093 'Reading mspps\n'); 0094 scalefactor = 1000/double(cell2mat(hdfread(mspps_file, 'IWP_SCAL'))); 0095 MSPPS_IWP_raw = hdfread(mspps_file, 'IWP'); 0096 MSPPS_IWP = ((MSPPS_IWP_raw<0) + (MSPPS_IWP_raw>0)*scalefactor) .* double(MSPPS_IWP_raw); 0097 assert(size(MSPPS_IWP, 1)==size(data2.tb, 1), ... 0098 'atmlab:colloc_process_data_cpr_mhs', ... 0099 'MSPPS has not the same number of rows as AMSU. Skipping'); 0100 M(:, c.data.MSPPS_IWP) = MSPPS_IWP(i2); 0101 catch ME 0102 switch (ME.identifier) 0103 case {'MATLAB:HDF:invalidFile', 'atmlab:find_datafile_by_date', ... 0104 'atmlab:colloc_process_data_cpr_mhs'} 0105 logtext(atmlab('ERR'), ... 0106 'WARNING: Unable to process mspps: %s\n', ME.message); 0107 otherwise 0108 ME.rethrow(); 0109 end 0110 end 0111 0112