read_saphir_l1

 function read_saphir_l1: reads HDF5 files with SAPHIR leve1 1 data and
 provides a structure containing the useful L1 data and attributes

 call: dataL1 = read_saphir_l1(filename,[selected_entries, mode])

 filename is the name of the hdf5 file or a gzipped version thereof
 desired_entries is optional if you do not want to extract all variables
 but only a selection. See variable 'strmember' for entries.

 relate variable names structure vs hdf5

 Initial version mostly by Mathias Milz, adapted by Gerrit Holl

 'mode' can be 'l1a1' or 'l1a2'.

0001 function saphirL1=read_saphir_l1(filename,selected_entries, mode)
0002 % function read_saphir_l1: reads HDF5 files with SAPHIR leve1 1 data and
0003 % provides a structure containing the useful L1 data and attributes
0004 %
0005 % call: dataL1 = read_saphir_l1(filename,[selected_entries, mode])
0006 %
0007 % filename is the name of the hdf5 file or a gzipped version thereof
0008 % desired_entries is optional if you do not want to extract all variables
0009 % but only a selection. See variable 'strmember' for entries.
0010 %
0011 % relate variable names structure vs hdf5
0012 %
0013 % Initial version mostly by Mathias Milz, adapted by Gerrit Holl
0014 %
0015 % 'mode' can be 'l1a1' or 'l1a2'.
0016 
0017 % $Id: read_saphir_l1.m 8924 2014-07-30 07:50:32Z mmilz $
0018 
0019 errorid = 'Function read_saphir_l1:';
0020 
0021 % structure variable names for output
0022 strmember= {'IncidencAngle',...
0023     'LatNadir',...
0024     'LatScan',...
0025     'LonNadir',...
0026     'LonScan',...
0027     'QFCh1',...
0028     'QFCh2',...
0029     'QFCh3',...
0030     'QFCh4',...
0031     'QFCh5',...
0032     'QFCh6',...
0033     'QFScan',...
0034     'ScanTimestart',...
0035     'ScanGain',...
0036     'ScanHLTemp',...
0037     'ScanNumber',...
0038     'ScanOffset',...
0039     'TbCh1',...
0040     'TbCh2',...
0041     'TbCh3',...
0042     'TbCh4',...
0043     'TbCh5',...
0044     'TbCh6'};
0045 
0046 % variable names for h5 internal structure
0047 h5member= {'/ScienceData/IncidenceAngle_Samples',...
0048     '/ScienceData/Latitude_Nadir',...
0049     '/ScienceData/Latitude_Samples',...
0050     '/ScienceData/Longitude_Nadir',...
0051     '/ScienceData/Longitude_Samples',...
0052     '/ScienceData/QF_Samples_S1',...
0053     '/ScienceData/QF_Samples_S2',...
0054     '/ScienceData/QF_Samples_S3',...
0055     '/ScienceData/QF_Samples_S4',...
0056     '/ScienceData/QF_Samples_S5',...
0057     '/ScienceData/QF_Samples_S6',...
0058     '/ScienceData/SAPHIR_QF_scan',...
0059     '/ScienceData/Scan_FirstSampleAcqTime',...
0060     '/ScienceData/Scan_Gain',...
0061     '/ScienceData/Scan_HotLoadTemperature',...
0062     '/ScienceData/Scan_Number',...
0063     '/ScienceData/Scan_Offset',...
0064     '/ScienceData/TB_Samples_S1',...
0065     '/ScienceData/TB_Samples_S2',...
0066     '/ScienceData/TB_Samples_S3',...
0067     '/ScienceData/TB_Samples_S4',...
0068     '/ScienceData/TB_Samples_S5',...
0069     '/ScienceData/TB_Samples_S6'};
0070 
0071 switch nargin
0072     case 0
0073         error(errorid,'You need at least Filename as input');
0074     case 1
0075         % read all variables
0076         useentries = 1:length(strmember);
0077     case {2, 3}
0078         %read only selected variables given in 'selected_entries'
0079         [~,useentries] = ismember(selected_entries,strmember);
0080 end
0081 
0082 if nargin < 3
0083     [~,name,~] = fileparts(filename);
0084     switch name(8:11)
0085         case 'L1A_'
0086             mode = 'l1a1';
0087         case 'L1A2'
0088             mode = 'l1a2';
0089         otherwise
0090             error(['atmlab:' mfilename ':modeunknown'], ...
0091                 'Can''determine mode. Set mode manually to ''l1a1'' or ''l1a2''');
0092     end
0093 end
0094 
0095 switch mode
0096     case 'l1a1'
0097         % above is good
0098     case 'l1a2'
0099         h5member = strrep(h5member, 'IncidenceAngle_Samples', 'Incidence_Angle');
0100         h5member = strrep(h5member, 'TB_Samples', 'TB_Pixels');
0101         h5member = strrep(h5member, 'QF_Samples', 'QF_Pixels');
0102         h5member = strrep(h5member, 'Samples', 'pixels');
0103         h5member = strrep(h5member, 'Sample', 'Pixel');
0104     otherwise
0105         error(['atmlab:' mfilename ':invalidmode'], ...
0106             'Expected mode ''l1a1'' or ''l1a2'', got %s', mode);
0107 end
0108 
0109 % uncompress file if necessary
0110 try
0111     info=hdf5info(filename);
0112 catch ME
0113     switch ME.identifier
0114         case 'MATLAB:imagesci:deprecatedHDF5:libraryError' % probably compressed
0115             filename = uncompress(filename, atmlab('WORK_AREA'), ...
0116                 struct('unidentified', 'error'));
0117             % if we reach here, it means filename uncompressed is not the
0118             % original filename!
0119             cleanupObj = onCleanup(@()delete(filename));
0120             info = h5info(filename);
0121         otherwise
0122             ME.rethrow();
0123     end
0124 end
0125 %Begin MMILZ
0126 % read version from info and adjust h5member names (Latitude_pixels -> Latitude_Pixels
0127 % (1.05->1.06))
0128 info=h5info(filename); %Call necessary to enable the following tastks
0129 natt=length(info.Groups.Attributes);
0130 for iat=1:natt
0131     if strcmp(info.Groups.Attributes(iat).Name,'Product_Identification')
0132         fname=info.Groups.Attributes(iat).Value;
0133         vers=fname(13:16);
0134         if strcmp(vers,'1.06')
0135             %fprintf(1,'Ping\n');
0136             h5member = strrep(h5member, 'Latitude_pixels', 'Latitude_Pixels');
0137             h5member = strrep(h5member, 'Longitude_pixels', 'Longitude_Pixels');
0138         end
0139     end
0140 end
0141 clear natt;
0142 clear iat;
0143 clear fname;
0144 clear vers;
0145 %END MMILZ
0146 
0147 % Begin MMILZ 20140730
0148 % Catch the change in Variable naming Icidence_angle -->
0149 % IncidenceAngle_Pixels and adapt the fieldnames in
0150 ndatasets=length(info.Groups.Datasets);
0151 for idatas=1:ndatasets
0152     if strcmp(info.Groups.Datasets(idatas).Name,'IncidenceAngle_Pixels')
0153         h5member = strrep(h5member, 'Incidence_Angle', 'IncidenceAngle_Pixels');
0154     end
0155 end
0156 clear ndatasets;
0157 clear idatas;
0158 %%MM %END MMILZ 20140730
0159 
0160 
0161 for ie=1:length(useentries)
0162     % read data
0163     % use h5read, not hdf5read, as to not get annoying hdf5.h5string object
0164     %    saphirL1.data.(strmember{useentries(ie)}) = hdf5read(filename,h5member{useentries(ie)});
0165     saphirL1.data.(strmember{useentries(ie)}) = h5read(filename,h5member{useentries(ie)});
0166     %fprintf(1,'Read Variable: %s\n',strmember{useentries(ie)});
0167     % read data attributes
0168     %fname={info.GroupHierarchy.Groups.Datasets(useentries(ie)).Attributes.Shortname};
0169     fname = {info.Groups.Datasets(useentries(ie)).Attributes.Name};
0170     for i =1:length(fname)
0171         %fprintf(1,'i: %g -- %s\n',i,fname{i});
0172         try
0173             %val = info.GroupHierarchy.Groups.Datasets(useentries(ie)).Attributes(i).Value.Data;
0174             val = info.Groups.Datasets(useentries(ie)).Attributes(i).Value;
0175             
0176             saphirL1.attr.(strmember{useentries(ie)}).(strrep(fname{i},'_',''))= val;
0177         catch ME
0178             logtext(atmlab('ERR'), ...
0179                 'Error in reading data attribute %s for field %s: %s\n', ...
0180                 fname{i}, strmember{useentries(ie)}, ME.message);
0181         end
0182     end
0183     % Get absolute values by applying offset and scale factor
0184     if isfield(saphirL1.attr.(strmember{useentries(ie)}),'scalefactor')
0185         tmpscalefactor=str2num(saphirL1.attr.(strmember{useentries(ie)}).scalefactor);
0186         saphirL1.data.(strmember{useentries(ie)}) = double(saphirL1.data.(strmember{useentries(ie)})) * ...
0187             double(tmpscalefactor);
0188         % Do nothing
0189         % else
0190         %    Do something
0191     end
0192     
0193     if isfield(saphirL1.attr.(strmember{useentries(ie)}),'addoffset')
0194         tmpaddoffset=str2num(saphirL1.attr.(strmember{useentries(ie)}).addoffset);
0195         saphirL1.data.(strmember{useentries(ie)}) = double(saphirL1.data.(strmember{useentries(ie)})) + ...
0196             tmpaddoffset;
0197         % Do nothing
0198         % else
0199         %    Do something
0200     end
0201 end
0202 % Include quality flags.
0203 % If all data are read, use QFlags from structure
0204 % If QFlags are not read, read QFs explicitly and apply
0205 % Channel 1
0206 for ie=6:11
0207     if isfield(saphirL1.data,strmember{ie+12})
0208         if ~isfield(saphirL1.data,strmember{ie})
0209             %            saphirL1.data.(strmember{ie})=  hdf5read(filename,h5member{ie});
0210             saphirL1.data.(strmember{ie})=  h5read(filename,h5member{ie});
0211         end
0212         [QFCh1validtmp,lsChtmp,lsmixChtmp]=getvalidpixels(saphirL1.data.(strmember{ie}));
0213         % set invalid entries to NaN
0214         saphirL1.data.(strmember{ie+12})(~QFCh1validtmp) = NaN;
0215     end
0216 end
0217 %Land Sea? Quality Flag Scan? --> if this is required, the according HDF5
0218 %Entrie can be read and used strmember{12}
0219 end
0220 
0221 function [QFvalid,landsea,landseamix]=getvalidpixels(arrayQFCh)
0222 % getvalidpixels scans all uint16 bits of the qualityflags to reurn only
0223 % valid pixel
0224 %
0225 % getvalidpixels extracts the individual bits from the uint16 Quality Flag
0226 % array of each channel and determines the validity of the pixels.
0227 %
0228 % input: arrayQFCh: array with uint16 containing the bitwise quality information
0229 % output: QFvalid: array with true for valid bixels and false for invalid pixels, affected by any undesired error source.
0230 %         landsea: mask defining land or sea surface
0231 %         landseamix: mask defining mixes surfaces
0232 %
0233 % By Mathias Milz
0234 
0235 %
0236 % Access individual bits using 'bitget(A,BIT)'
0237 %
0238 % Bits 1, 2: Ice all values accepted
0239 flag1 = bitget(arrayQFCh,1) | ~bitget(arrayQFCh,1) | bitget(arrayQFCh,2)...
0240     | ~bitget(arrayQFCh,2);
0241 % Bits 3 = 0 blank
0242 flag2 = ~bitget(arrayQFCh,3);
0243 % BITS 4: Interpolation quality: should be 0
0244 flag3 = ~bitget(arrayQFCh,4);
0245 % BITS 5: ColdSky count error: should be 0
0246 flag4 = ~bitget(arrayQFCh,5);
0247 % BITS 6: Hoot Count error: should be 0
0248 flag5 = ~bitget(arrayQFCh,6);
0249 % BITS 7 and 8 : Calibration Flag: should be 0
0250 flag6 = ~bitget(arrayQFCh,7) & ~bitget(arrayQFCh,8);
0251 % BITS 9: Geo-location stimaton should be 0
0252 flag7 = ~bitget(arrayQFCh,9);
0253 % BITS 10: Level-0 Count poor value: should be 0
0254 flag8 = ~bitget(arrayQFCh,10);
0255 % BITS 11: Level-0 Count saturated: should be 0
0256 flag9 = ~bitget(arrayQFCh,11);
0257 % BITS 12: ON OFF Channel flag: should be 0
0258 flag10 = ~bitget(arrayQFCh,12);
0259 % BITS 13: SurfaceType: 0 Sea, 1 Land, both accepted
0260 flag11 = ~bitget(arrayQFCh,13) | bitget(arrayQFCh,13);
0261 landsea = bitget(arrayQFCh,13);
0262 % BITS 14: LandSea contamination: 0 no contamination, 1 contamination, both accepted
0263 flag12 = ~bitget(arrayQFCh,14) | bitget(arrayQFCh,14);
0264 landseamix = bitget(arrayQFCh,14);
0265 % BITS 15: Sun Glint: both accepted
0266 flag13 = ~bitget(arrayQFCh,15) | bitget(arrayQFCh,15);
0267 % BITS 16: TB validy: should be 0
0268 flag14 = ~bitget(arrayQFCh,16);
0269 QFvalid = flag1 & flag2 & flag3 & flag4 & flag5 & flag6 & flag7 & flag8 & ...
0270     flag9 & flag10 & flag11 & flag12 & flag13 & flag14;
0271 end