Home > atmlab > datasets > read_isccp.m

read_isccp

PURPOSE ^

% isccpreadnative

SYNOPSIS ^

function S = read_isccp(file, opt)

DESCRIPTION ^

% isccpreadnative
 Purpose: read ISSCP binary data (based on the fortran reading rountines)

 IN:    file   %s           ISCCP file (fullpath)
        opt    struct       contains     e.g.
                               'dataset'    'dx'
                               'cols'       [32,54,...]      (D2 dataset only)

 OUT:   S      structure with all the data from the file

 USAGE: S = read_isccp (file, opt)

 Salomon Eliasson
 $Id: read_isccp.m 8372 2013-04-24 06:23:43Z seliasson $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

read_isccp.m

SOURCE CODE ^

0001 function S = read_isccp(file, opt)
0002 %% isccpreadnative
0003 % Purpose: read ISSCP binary data (based on the fortran reading rountines)
0004 %
0005 % IN:    file   %s           ISCCP file (fullpath)
0006 %        opt    struct       contains     e.g.
0007 %                               'dataset'    'dx'
0008 %                               'cols'       [32,54,...]      (D2 dataset only)
0009 %
0010 % OUT:   S      structure with all the data from the file
0011 %
0012 % USAGE: S = read_isccp (file, opt)
0013 %
0014 % Salomon Eliasson
0015 % $Id: read_isccp.m 8372 2013-04-24 06:23:43Z seliasson $
0016 
0017 file = uncompress(file);
0018 
0019 fid = fopen (file,'rb','b');
0020 cleanObj = onCleanup(@() fclose(fid));
0021 logtext(atmlab('OUT'),'Processing file: %s\n',file)
0022 switch opt.dataset
0023     case 'dx'
0024         
0025         %Each data record consists of a
0026         %prefix area containing location information followed by a data area containing the packed pixel data,
0027         %followed by padding to the end of the record (pad value = 255). A data record contains a variable
0028         %number of whole pixels. All pixel data are reported as single byte values from 0 to 255, with 255
0029         %reserved to represent missing data. Contents of the header record are given in Table 2.3.3 and of the
0030         %data record in Table 2.3.4.
0031         
0032         rlength = 30720; % record length [bytes]
0033         % ------------------------
0034         % GET HEADER
0035         % ------------------------
0036         %Year, Month, Day, UTC (1-8), Satellite ID code, Satellite
0037         %Type Code (SATTYP), Number of channels present (NCHANS), Night Image flag
0038         %logtext(atmlab('OUT'),'Processing record #1: HEADER\n')
0039         info = str2num(fread (fid, 80, 'uint8=>char')'); %#ok<ST2NM>
0040         
0041         S.global_attributes = struct('year',info(1),'month',info(2),'day',info(3),'utc',info(4),...
0042             'satID',info(5),'SATTYP',info(6),'NCHANS',info(7),'nightFlag',info(8),...
0043             'production_time',strtrim(fread (fid, 80, 'uint8=>char')'),...
0044             'production_date',strtrim(fread (fid, 80, 'uint8=>char')'));
0045         clear info
0046         % move to the first byte of line 101 (skipping blanks)
0047         [~]=fseek(fid,100*80,-1);
0048         
0049         % Title, units, format for selected DX variables
0050         S.global_attributes.titleETC = fread (fid, 80*(1+(135-101)),'uint8=>char')';
0051         S.attributes = DX_attributes();
0052         
0053         % move to the first record of the data (skipping blanks)
0054         [~]=fseek(fid,rlength,-1); % i.e., end of first record
0055         
0056         % ----------------
0057         % READ DATA
0058         % ----------------
0059         
0060         tic
0061         % preallocate
0062         npreall = 2e5;
0063         T = LookupTables(opt.dataset); % Do all the bitand's in one go
0064         
0065         S.lat = ones(1,npreall,'single')*-9;
0066         S.lon = ones(1,npreall,'single')*-9;
0067         [S.LNDWTR,S.HITOPO,S.SNOICE,S.MUE,S.IRAD,S.BXICSR,S.GLINT,S.MU0,S.PHI,S.VRAD,...
0068             S.DAYNIT,S.ITHR,S.VTHR,S.SHORE,S.ICSRAD,S.ITMP,S.IPRS,S.ICSTMP,S.ICSPRS,...
0069             S.NREF,S.NTHR,S.NCSREF,S.VCSRAD,S.VALBTA,S.VCSALB,S.VTMP,S.VPRS,S.VTAUIC,...
0070             S.VTMPIC,S.VPRSIC,S.CLOUD,S.CLOUD_PHASE] = deal(ones(1,npreall,'uint8')*255);
0071         S.ARAD = ones(3,npreall,'uint8')*255;
0072         
0073         % Since CLOUD,CLOUD_PHASE is calculated in DXREAD (fortran original), it is
0074         % calculated here too.
0075         
0076         record=2; sNp =1;pix=1;
0077         while ~feof(fid) % while not eof (i.e., loop over records)
0078 %            logtext(atmlab('OUT'),'Processing record #%g. ',record)
0079             WM = fread (fid, 1, 'uint32')/10; % Western-most longitude (0-3600 degrees*10)
0080             if isempty(WM), break, end
0081 %             EM = fread (fid, 1, 'uint32')/10; % Eastern-most longitude (0-3600 degrees*10)
0082 %             NM = fread (fid, 1, 'uint32')/10 -90; % Northern-most latitude (0-1800 degrees*10)
0083 %             SM = fread (fid, 1, 'uint32')/10 -90; % Southern-most latitude (0-1800 degrees*10)
0084 %             fprintf(atmlab('OUT'),'Region: lons between %.1f and %.1f and lats between %.1f and %.1f\n',WM,EM,SM,NM) %#ok<PRTCAL>
0085             [~]=fseek(fid,3*4,0);
0086             n_pixels     = fread (fid, 1, 'uint32'); %Number of pixels reported in data area
0087             iout         = fread (fid, 1, 'uint32'); %Number of bytes in packed data BUFFER
0088             
0089             % ---------------
0090             % GEODATA
0091             % ---------------
0092             eNp = sNp+n_pixels-1; % last index for pixel
0093             
0094             geo = reshape(fread( fid, n_pixels*4, 'uint16'), 4, [])';
0095             S.lon(sNp:eNp) = geo(:,1)/10;
0096             S.lat(sNp:eNp) = geo(:,2)/10 -90;
0097             sNp = eNp+1;
0098             % READ THE REST OF THE DATA AS unit8 and loop over the pixels
0099             buffer = fread(fid, iout, 'uint8=>uint8'); % Packed pixel data for NPIX pixels
0100             beginB=1;
0101             for P = 1:n_pixels
0102                 
0103                 [S1,ADD1,S2,S3,ADD3,S4,beginB]=get_dx_sections(buffer,beginB,S.global_attributes,T);
0104                 
0105                 % -----
0106                 % S1. ALWAYS
0107                 % -----
0108                 % '*' means for internal use only
0109                 %Land/water flag (0-1), 1 = water  pixel
0110                 %S.LNDWTR(pix)  = bitshift(bitand(S1(1),2^5,'uint8'),-5);
0111                 S.LNDWTR(pix)  = T.LNDTAB(S1(1)+1);
0112                 %Topography flag (0-1), 1 = high topography pixel
0113                 S.HITOPO(pix)  = T.TOPTAB(S1(1)+1);
0114                 %Snow/ice code (0-3) (See Table 2.5.10)
0115                 S.SNOICE(pix)  = T.SNITAB(S1(1)+1);
0116                 %Cosine of satellite zenith angle * 100 (0-100)
0117                 S.MUE(pix)     = T.ANGTAB(S1(3)+1);
0118                 %IR radiance (0-254 counts)
0119                 S.IRAD(pix)    = S1(4);
0120                 %First IR clear sky radiance (0-254 counts)
0121                 S.BXICSR(pix)  = S1(5);
0122                 
0123                 
0124                 
0125                 % ----------------------
0126                 % ADD1. depends on number of channels
0127                 for i = 1:3
0128                     %1) First extra wavelength radiance (0-254 counts) (present only when NCHANS > 2)
0129                     %2) Second extra wavelength radiance (0-254 counts)(present only when NCHANS > 3)
0130                     %3) Third extra wavelength radiance (0-254 counts)(present only when NCHANS > 4)
0131                     S.ARAD(i,pix) = ADD1(i);
0132                     
0133                 end
0134                 
0135                 % ----
0136                 % S2
0137                 % ----
0138                 if ~isempty(S2)
0139                     S.GLINT(pix)  = T.GLNTAB(S2(1)+1);%Glint flag (0-1), 1 = glint condition exists
0140                     S.MU0(pix)    = T.ANGTAB(S2(2)+1);%Cosine of solar zenith angle * 100 (0-100)
0141                     S.PHI(pix)    = S2(3);%Relative azimuth angle (0-180 degrees)
0142                     S.VRAD(pix)   = S2(4);%VIS radiance (0-254 counts)
0143                 end
0144                 
0145                 % ----
0146                 % S3
0147                 % ----
0148                 %Day/Night flag (0-1), 1 = night pixel (no VIS)
0149                 S.DAYNIT(pix)   = T.NITTAB(S3(1)+1);
0150                 %Final IR threshold result (0-5), 4,5 = cloudy
0151                 S.ITHR(pix)     = T.ITHRTB(S3(1)+1);
0152                 %Final VIS threshold result (0-5), 4,5 = cloudy
0153                 S.VTHR(pix)     = T.VTHRTB(S3(1)+1);
0154                 %Shore flag (0-1), 1 = near-coastal pixel
0155                 S.SHORE(pix)    = T.NSHTAB(S3(1)+1);
0156                 %IR clear sky composite radiance (0-254 counts)
0157                 S.ICSRAD(pix)   = S3(3);
0158                 %IR-retrieved cloud top or surface temperature     (0-254 counts)
0159                 S.ITMP(pix)     = S3(4);
0160                 %IR-retrieved cloud top or surface pressure (0-254 counts)
0161                 S.IPRS(pix)     = S3(5);
0162                 %IR-retrieved clear sky composite temperature (0-254 counts)
0163                 S.ICSTMP(pix)   = S3(6);
0164                 %IR-retrieved clear sky composite pressure (0-254 counts)
0165                 S.ICSPRS(pix)   = S3(7);
0166                 
0167                 % -------
0168                 % ADD3 (polar-orbiting Satellites)
0169                 % -------
0170                 if ~isempty(ADD3)
0171                     S.NREF(pix)  = ADD3(1); %NIR reflectivity (0-254 counts)
0172                     S.NTHR(pix)  = T.CRETAB(ADD3(2)+1); %NIR threshold result (1-13), > 8 = cloudy
0173                     S.NCSREF(pix)= ADD3(3); %NIR clear sky composite reflectance (0-254 counts)
0174                 end
0175                 
0176                 
0177                 % ----
0178                 % S4
0179                 % ----
0180                 if ~isempty(S4)
0181                     %VIS clear sky composite radiance (0-254 counts)
0182                     S.VCSRAD(pix) = S4(2);
0183                     % VIS-retrieved liquid cloud tau or surface refl. (0-254 counts)
0184                     S.VALBTA(pix) = S4(3);
0185                     % VIS-retrieved clear sky composite reflectance (0-254 counts)
0186                     S.VCSALB(pix) = S4(4);
0187                     % VIS-adjusted cloud top temperature (0-254 counts)
0188                     S.VTMP(pix)   = S4(5);
0189                     % VIS-adjusted cloud top pressure (0-254 counts)
0190                     S.VPRS(pix)   = S4(6);
0191                     % VIS-retrieved ice cloud tau (0-254 counts)
0192                     S.VTAUIC(pix) = S4(7);
0193                     % VIS-adjusted ice cloud top temperature (0-254 counts)
0194                     S.VTMPIC(pix) = S4(8);
0195                     % VIS-adjusted ice cloud top pressure (0-254 counts)
0196                     S.VPRSIC(pix) = S4(9);
0197                 end
0198                 
0199                 % CLOUD=1 /* pixel is cloudy */; CLOUD=0 /* pixel is clear */
0200                 S.CLOUD(pix) = ( S.ITHR(pix) > 3 | S.VTHR(pix) > 3 | S.NTHR(pix) > 8 );
0201                 
0202                 %The actual water cloud - ice cloud decision may be calculated by the user as follows:
0203                 % CLOUD_PHASE = [0 or 1 or 255] (ice,water,not determined)
0204                 S.CLOUD_PHASE(pix) = S.VTMPIC(pix) == 255 | ( S.VTMP(pix) > 74 & S.VTMPIC(pix) > 74 );
0205                 pix = pix+1;
0206             end
0207             
0208             [~]=fseek(fid,record*rlength,'bof');
0209             
0210             record = record+1;
0211         end
0212         
0213         % TRUNCATE
0214         S = truncateS(S,eNp);
0215         
0216         % ------------
0217         % Apply tabulated values
0218         % ------------
0219 
0220         % temperature variables
0221         S.IRAD  = T.TMPTAB(S.IRAD+1);
0222         S.BXICSR  = T.TMPTAB(S.BXICSR+1);
0223         %ARAD    = TMPTAB(ARAD); Don't know if they are visible or IR
0224         S.VRAD    = T.RFLTAB(S.VRAD+1);
0225         S.ICSRAD  = T.TMPTAB(S.ICSRAD+1);
0226         S.ITMP    = T.TMPTAB(S.ITMP+1);
0227         S.IPRS    = T.PRETAB(S.IPRS+1);
0228         S.ICSTMP  = T.TMPTAB(S.ICSTMP+1);
0229         S.ICSPRS  = T.PRETAB(S.ICSPRS+1);
0230         S.NREF    = T.RFLTAB(S.NREF+1);
0231         S.NCSREF  = T.RFLTAB(S.NCSREF+1);
0232         S.VCSRAD  = T.RFLTAB(S.VCSRAD+1);
0233         S.VALBTA(S.CLOUD_PHASE==1) = T.TAUTAB(S.VALBTA(S.CLOUD_PHASE==1)+1);
0234         S.VALBTA(S.CLOUD==0) = T.RFLTAB(S.VALBTA(S.CLOUD==0)+1);
0235         S.VCSALB  = T.RFLTAB(S.VCSALB+1);
0236         S.VTMP    = T.TMPTAB(S.VTMP+1);
0237         S.VPRS    = T.PRETAB(S.VPRS+1);
0238         S.VTAUIC  = T.TAUTAB(S.VTAUIC+1);
0239         S.IWP     = S.VTAUIC*T.PATHIW; %convert tau_v to IWP following D1 and D2 dataset definitions (Sec. 3.1.2 ISCCP documentation new cloud datasets Jan 1996)
0240         S.VTMPIC  = T.TMPTAB(S.VTMPIC+1);
0241         S.VPRSIC  = T.PRETAB(S.VPRSIC+1);
0242         
0243         if npreall<eNp, logtext(atmlab('ERR'),'Need to preallocate to bigger arrays. eNp=%g',eNp), end
0244         toc
0245         
0246     case 'd1'
0247         error(['atmlab:' mfilename],'Nothing set up for the D1 dataset. Feel free to set it up')
0248          % nothing set up yet
0249     case 'd2'
0250         
0251         warning(['atmlab:' mfilename ':legacy'],['The output format doesn''t match the documentation.\n'...
0252             'this is the legacy format from an earlier less generic code.\n'...
0253             'The plan is to fix this part of the reading routine when I (or someone else has time'])
0254         fid = fopen (filename,'rb');
0255         S = uint8(fread (fid));
0256         S = reshape(S,13000,numel(S)/13000)';     % sort to 13000 bytes for each grid
0257         S = S(:,131:end);                         % I don't need prefix
0258         S = reshape(S',numel(S),1);
0259         S = S(1:end-4810);                        % removed bytes not used
0260         S = reshape(S,130,numel(S)/130)';         % number of data fields = 130
0261         S = S (:,cols);
0262         
0263     otherwise
0264         error(['atmlab:' mfilename],'Nothing configured for ISCCP dataset: %s',opt.dataset)
0265 end
0266 logtext(atmlab('OUT'),'DONE\n')
0267 
0268 end
0269 
0270 function attr = DX_attributes()
0271 
0272 % Assign some attributes to the read variables
0273 
0274 attr = struct('lon','Longitudes for NPIX pixels (0-3600 degrees*10)',...
0275     'lat','Latitudes for NPIX pixels (0-1800 degrees*10)',...
0276     'LNDWTR','Land/water flag (0-1), 1 = water  pixel',...
0277     'HITOPO','Topography flag (0-1), 1 = high topography pixel',...
0278     'SNOICE','Snow/ice code (0-3) (See Table 2.5.10)',...
0279     'MUE','Cosine of satellite zenith angle * 100 (0-100)',...
0280     'IRAD','IR radiance',...
0281     'BXICSR','First IR clear sky radiance',...
0282     'ARAD',['1) First extra wavelength radiance (present only when NCHANS > 2)'...
0283     '2) Second extra wavelength radiance (present only when NCHANS > 3)'...
0284     '3) Third extra wavelength radiance (present only when NCHANS > 4)'],...
0285     'GLINT','Glint flag (0-1), 1 = glint condition exists',...
0286     'MU0','Cosine of solar zenith angle * 100 (0-100)',...
0287     'PHI','Relative azimuth angle (0-180 degrees)',...
0288     'VRAD','VIS radiance',...
0289     'DAYNIT','Day/Night flag (0-1), 1 = night pixel (no VIS)',...
0290     'ITHR','Final IR threshold result (0-5), 4,5 = cloudy',...
0291     'VTHR','Final VIS threshold result (0-5), 4,5 = cloudy',...
0292     'SHORE','Shore flag (0-1), 1 = near-coastal pixel',...
0293     'ICSRAD','IR clear sky composite radiance',...
0294     'ITMP','IR-retrieved cloud top or surface temperature',...
0295     'IPRS','IR-retrieved cloud top or surface pressure',...
0296     'ICSTMP','IR-retrieved clear sky composite temperature',...
0297     'ICSPRS','IR-retrieved clear sky composite pressure',...
0298     'NREF','NIR reflectivity',...
0299     'NTHR','NIR threshold result (1-13), > 8 = cloudy',...
0300     'NCSREF','NIR clear sky composite reflectance',...
0301     'VCSRAD','VIS clear sky composite radiance',...
0302     'VALBTA','VIS-retrieved liquid cloud tau or surface reflectance',...
0303     'VCSALB','VIS-retrieved clear sky composite reflectance',...
0304     'VTMP',' VIS-adjusted cloud top temperature',...
0305     'VPRS','VIS-adjusted cloud top pressure',...
0306     'VTAUIC','VIS-retrieved ice cloud tau',...
0307     'VTMPIC','VIS-adjusted ice cloud top temperature',...
0308     'VPRSIC','VIS-adjusted ice cloud top pressure',...
0309     'CLOUD','CLOUD=1 /* pixel is cloudy */; CLOUD=0 /* pixel is clear */',...
0310     'CLOUD_PHASE',' [0 or 1 or 255] (ice,water,not determined)',...
0311     'IWP',['Conversion of tau_v to IWP following D1 and D2 dataset definitions '...
0312     '(Sec. 3.1.2 ISCCP documentation new cloud datasets Jan 1996)']);
0313 
0314 end
0315 
0316 function [S1,ADD1,S2,S3,ADD3,S4,bs]=get_dx_sections(buffer,bs,attr,T)
0317 % GET SECTIONS (as documented in the ISCCP DX documentation)
0318 
0319 % It's one second faster per isccp binary file file to not frivolously
0320 %assign variables for the sake of readibility. To clarify what going on ...
0321 %
0322 % bs = begin index
0323 % be = bs + number of bytes to read
0324 % NODAY = T.NITTAB(S1(1)+1) Night flag
0325 % NCHANS = (attr.NCHANS -2); more than 2 channels
0326 
0327 %S1
0328 S1 = buffer(bs:bs+5-1);
0329 bs=bs+5;
0330 
0331 
0332 %ADD1
0333 ADD1 = buffer(bs:bs+(attr.NCHANS -2)-1);
0334 bs = bs+(attr.NCHANS -2);
0335 
0336 %S2
0337 S2 = buffer(bs:bs+5*~(T.NITTAB(S1(1)+1))-1);
0338 bs = bs+5*~(T.NITTAB(S1(1)+1));
0339 
0340 %S3
0341 S3 = buffer(bs:bs+7-1);
0342 bs = bs+7;
0343 
0344 % ADD3
0345 ADD3=buffer(bs:bs+3*(attr.SATTYP>0)-1);
0346 bs = bs+3*(attr.SATTYP>0);
0347 
0348 % S4
0349 S4=buffer(bs:bs+9*~(T.NITTAB(S1(1)+1))-1);
0350 bs = bs+9*~(T.NITTAB(S1(1)+1));
0351 
0352 end
0353 
0354 function S = truncateS(S,eNp)
0355 
0356 fls = fieldnames(S);fls=fls(~ismember(fls,{'global_attributes','attributes'}));
0357 for F=fls'
0358     S.(F{1}) = S.(F{1})(:,1:eNp);
0359 end
0360 
0361 end
0362 
0363 function T = LookupTables(dataset)
0364 %% LookupTablesBitAnd
0365 % Make look up tables in the same manner as the Fortran code
0366 
0367 if strcmp(dataset,'dx')
0368     % --------------
0369     % For bitands
0370     [T.NITTAB,T.LNDTAB,T.TOPTAB,T.SNITAB,T.ANGTAB,T.ITHRTB,T.VTHRTB,...
0371         T.CRETAB,T.NSHTAB]=deal(zeros(1,255,'uint8'));
0372     
0373     for IVAL=1:256
0374         T.NITTAB(IVAL) = bitshift(bitand(IVAL,2^7),-7);
0375         T.LNDTAB(IVAL) = bitshift(bitand(IVAL,2^5),-5);
0376         T.TOPTAB(IVAL) = bitshift(bitand(IVAL,2^4),-4);
0377         T.SNITAB(IVAL) = bitshift(bitand(IVAL,2^3+2^2),-2);
0378         T.ANGTAB(IVAL) = bitand(IVAL,2^6+2^5+2^4+2^3+2^2+2^1+2^0);
0379         T.ITHRTB(IVAL) = bitshift(bitand(IVAL,2^6+2^5+2^4),-4);
0380         T.VTHRTB(IVAL) = bitshift(bitand(IVAL,2^3+2^2+2^1),-1);
0381         T.CRETAB(IVAL) = bitand(IVAL,2^3+2^2+2^1+2^0);
0382         T.NSHTAB(IVAL) = bitand(IVAL,2^0);
0383     end
0384     T.GLNTAB = T.NITTAB;
0385 end
0386 
0387 % -------------
0388 % CONVERSION tables
0389 
0390 % Temperature conversion tables
0391 T.TMPTAB = [-100.000,165.000,169.000,172.000,175.000,177.800,180.500,183.000, ...
0392     185.500,187.800,190.000,192.000,194.000,195.700,197.500,199.200, ...
0393     201.000,202.700,204.500,206.200,208.000,209.700,211.500,212.800, ...
0394     214.100,215.400,216.700,217.900,219.200,220.500,221.800,223.100, ...
0395     224.400,225.400,226.500,227.500,228.600,229.600,230.600,231.700, ...
0396     232.700,233.800,234.800,235.700,236.600,237.500,238.400,239.200, ...
0397     240.100,241.000,241.900,242.800,243.700,244.500,245.300,246.100, ...
0398     246.900,247.700,248.500,249.300,250.100,250.900,251.700,252.400, ...
0399     253.100,253.900,254.600,255.300,256.000,256.700,257.500,258.200, ...
0400     258.900,259.500,260.200,260.800,261.500,262.100,262.800,263.400, ...
0401     264.100,264.700,265.400,266.000,266.600,267.200,267.800,268.400, ...
0402     269.100,269.700,270.300,270.900,271.500,272.100,272.700,273.200, ...
0403     273.800,274.400,275.000,275.600,276.100,276.700,277.300,277.800, ...
0404     278.400,278.900,279.500,280.000,280.500,281.100,281.600,282.200, ...
0405     282.700,283.200,283.700,284.200,284.700,285.200,285.800,286.300, ...
0406     286.800,287.300,287.800,288.300,288.800,289.300,289.800,290.200, ...
0407     290.700,291.200,291.700,292.200,292.700,293.200,293.600,294.100, ...
0408     294.600,295.000,295.500,296.000,296.500,296.900,297.400,297.800, ...
0409     298.300,298.700,299.200,299.600,300.100,300.500,301.000,301.400, ...
0410     301.900,302.300,302.800,303.200,303.600,304.000,304.500,304.900, ...
0411     305.300,305.800,306.200,306.600,307.000,307.500,307.900,308.300, ...
0412     308.700,309.100,309.600,310.000,310.400,310.800,311.200,311.600, ...
0413     312.000,312.400,312.900,313.300,313.700,314.100,314.500,314.900, ...
0414     315.300,315.700,316.100,316.400,316.800,317.200,317.600,318.000, ...
0415     318.400,318.800,319.200,319.500,319.900,320.300,320.700,321.100, ...
0416     321.400,321.800,322.200,322.600,323.000,323.300,323.700,324.100, ...
0417     324.500,324.900,325.200,325.600,326.000,326.400,326.700,327.100, ...
0418     327.400,327.800,328.200,328.500,328.900,329.200,329.600,329.900, ...
0419     330.300,330.600,331.000,331.300,331.700,332.000,332.400,332.700, ...
0420     333.100,333.400,333.800,334.100,334.500,334.800,335.200,335.500, ...
0421     335.900,336.200,336.600,336.900,337.300,337.600,338.000,338.300, ...
0422     338.600,339.000,339.300,339.700,340.000,345.000,-200.000,-1000.000];
0423 
0424 T.RFLTAB = [-100.000,0.000,0.008,0.012,0.016,0.020,0.024,0.028,0.032,0.036, ...
0425      0.040,0.044,0.048,0.052,0.056,0.060,0.064,0.068,0.072,0.076, ...
0426      0.080,0.084,0.088,0.092,0.096,0.100,0.104,0.108,0.112,0.116, ...
0427      0.120,0.124,0.128,0.132,0.136,0.140,0.144,0.148,0.152,0.156, ...
0428      0.160,0.164,0.168,0.172,0.176,0.180,0.184,0.188,0.192,0.196, ...
0429      0.200,0.204,0.208,0.212,0.216,0.220,0.224,0.228,0.232,0.236, ...
0430      0.240,0.244,0.248,0.252,0.256,0.260,0.264,0.268,0.272,0.276, ...
0431      0.280,0.284,0.288,0.292,0.296,0.300,0.304,0.308,0.312,0.316, ...
0432      0.320,0.324,0.328,0.332,0.336,0.340,0.344,0.348,0.352,0.356, ...
0433      0.360,0.364,0.368,0.372,0.376,0.380,0.384,0.388,0.392,0.396, ...
0434      0.400,0.404,0.408,0.412,0.416,0.420,0.424,0.428,0.432,0.436, ...
0435      0.440,0.444,0.448,0.452,0.456,0.460,0.464,0.468,0.472,0.476, ...
0436      0.480,0.484,0.488,0.492,0.496,0.500,0.504,0.508,0.512,0.516, ...
0437      0.520,0.524,0.528,0.532,0.536,0.540,0.544,0.548,0.552,0.556, ...
0438      0.560,0.564,0.568,0.572,0.576,0.580,0.584,0.588,0.592,0.596, ...
0439      0.600,0.604,0.608,0.612,0.616,0.620,0.624,0.628,0.632,0.636, ...
0440      0.640,0.644,0.648,0.652,0.656,0.660,0.664,0.668,0.672,0.676, ...
0441      0.680,0.684,0.688,0.692,0.696,0.700,0.704,0.708,0.712,0.716, ...
0442      0.720,0.724,0.728,0.732,0.736,0.740,0.744,0.748,0.752,0.756, ...
0443      0.760,0.764,0.768,0.772,0.776,0.780,0.784,0.788,0.792,0.796, ...
0444      0.800,0.804,0.808,0.812,0.816,0.820,0.824,0.828,0.832,0.836, ...
0445      0.840,0.844,0.848,0.852,0.856,0.860,0.864,0.868,0.872,0.876, ...
0446      0.880,0.884,0.888,0.892,0.896,0.900,0.904,0.908,0.912,0.916, ...
0447      0.920,0.924,0.928,0.932,0.936,0.940,0.944,0.948,0.952,0.956, ...
0448      0.960,0.964,0.968,0.972,0.976,0.980,0.984,0.988,0.992,1.000, ...
0449      1.016,1.040,1.072,1.108,-200.000,-1000.000];
0450 
0451 % Pressure table
0452 T.PRETAB = [-100.00,1.00,5.00,10.00,15.00,20.00,25.00,30.00,35.00,40.00, ...
0453     45.00, 50.00, 55.00,60.00,65.00,70.00,75.00,80.00,85.00, ...
0454     90.00,95.00,100.00,105.00,110.00,115.00,120.00,125.00,130.00, ...
0455     135.00,140.00,145.00,150.00,155.00,160.00,165.00,170.00,175.00, ...
0456     180.00,185.00,190.00,195.00,200.00,205.00,210.00,215.00,220.00, ...
0457     225.00,230.00,235.00,240.00,245.00,250.00,255.00,260.00,265.00, ...
0458     270.00,275.00,280.00,285.00,290.00,295.00,300.00,305.00,310.00, ...
0459     315.00,320.00,325.00,330.00,335.00,340.00,345.00,350.00,355.00, ...
0460     360.00,365.00,370.00,375.00,380.00,385.00,390.00,395.00,400.00, ...
0461     405.00,410.00,415.00,420.00,425.00,430.00,435.00,440.00,445.00, ...
0462     450.00,455.00,460.00,465.00,470.00,475.00,480.00,485.00,490.00, ...
0463     495.00,500.00,505.00,510.00,515.00,520.00,525.00,530.00,535.00, ...
0464     540.00,545.00,550.00,555.00,560.00,565.00,570.00,575.00,580.00, ...
0465     585.00,590.00,595.00,600.00,605.00,610.00,615.00,620.00,625.00, ...
0466     630.00,635.00,640.00,645.00,650.00,655.00,660.00,665.00,670.00, ...
0467     675.00,680.00,685.00,690.00,695.00,700.00,705.00,710.00,715.00, ...
0468     720.00,725.00,730.00,735.00,740.00,745.00,750.00,755.00,760.00, ...
0469     765.00,770.00,775.00,780.00,785.00,790.00,795.00,800.00,805.00, ...
0470     810.00,815.00,820.00,825.00,830.00,835.00,840.00,845.00,850.00, ...
0471     855.00,860.00,865.00,870.00,875.00,880.00,885.00,890.00,895.00, ...
0472     900.00,905.00,910.00,915.00,920.00,925.00,930.00,935.00,940.00, ...
0473     945.00,950.00,955.00,960.00,965.00,970.00,975.00,980.00,985.00, ...
0474     990.00,995.00,1000.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0475     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0476     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0477     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0478     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0479     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0480     -200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00,-200.00, ...
0481     -1000.00];
0482 
0483 T.TAUTAB = [-100.000,0.020,0.040,0.060,0.090,0.110,0.140,0.160,0.190, ...
0484     0.220,0.240,0.270,0.300,0.330,0.370,0.400,0.430,0.460,0.500, ...
0485     0.530,0.570,0.600,0.640,0.680,0.720,0.750,0.790,0.830,0.870, ...
0486     0.920,0.960,1.000,1.040,1.090,1.130,1.180,1.220,1.270,1.320, ...
0487     1.370,1.420,1.470,1.520,1.570,1.620,1.670,1.730,1.780,1.830, ...
0488     1.890,1.950,2.000,2.060,2.120,2.180,2.240,2.300,2.360,2.430, ...
0489     2.490,2.550,2.620,2.690,2.750,2.820,2.890,2.960,3.030,3.100, ...
0490     3.180,3.250,3.320,3.400,3.480,3.550,3.630,3.710,3.790,3.880, ...
0491     3.960,4.040,4.130,4.220,4.300,4.390,4.480,4.570,4.670,4.760, ...
0492     4.850,4.950,5.050,5.150,5.250,5.350,5.450,5.560,5.660,5.770, ...
0493     5.880,5.990,6.110,6.220,6.340,6.450,6.570,6.690,6.820,6.940, ...
0494     7.070,7.190,7.330,7.460,7.590,7.730,7.870,8.010,8.150,8.300, ...
0495     8.440,8.590,8.740,8.900,9.060,9.220,9.380,9.540,9.710,9.880, ...
0496     10.050,10.230,10.410,10.590,10.780,10.970,11.160,11.350, ...
0497     11.550,11.760,11.960,12.170,12.390,12.600,12.830,13.050, ...
0498     13.280,13.520,13.760,14.000,14.250,14.510,14.770,15.030, ...
0499     15.300,15.580,15.860,16.150,16.440,16.740,17.050,17.360, ...
0500     17.690,18.020,18.350,18.700,19.050,19.410,19.780,20.160, ...
0501     20.540,20.940,21.350,21.770,22.200,22.630,23.080,23.550, ...
0502     24.030,24.520,25.020,25.540,26.070,26.620,27.190,27.770, ...
0503     28.370,28.990,29.630,30.290,30.970,31.670,32.400,33.160, ...
0504     33.940,34.740,35.580,36.450,37.350,38.290,39.260,40.260, ...
0505     41.320,42.420,43.570,44.760,46.000,47.310,48.680,50.110, ...
0506     51.600,53.170,54.840,56.590,58.430,60.360,62.400,64.590, ...
0507     66.900,69.360,71.960,74.720,77.730,80.940,84.380,88.060, ...
0508     92.020,96.400,101.010,105.510,109.870,114.330,119.590, ...
0509     125.920,133.660,143.120,154.650,169.560,187.490,207.200, ...
0510     228.130,250.440,282.780,323.920,378.650,-200.000,-200.000, ...
0511     -200.000,-200.000,-200.000,-200.000,-200.000,-200.000, ...
0512     -200.000,-200.000,-200.000,-1000.000];
0513 
0514 % CONSTANTS
0515 T.PATHLW = 6.292; %for LWP
0516 T.PATHIW = 10.5; %for the IWP
0517 
0518 if any(strcmp(dataset,{'d1','d2'}))
0519     % LATS & LONS for equal area grid
0520     table.latarray   = -90+2.5/2: 2.5 :90-2.5/2;
0521     table.lonlist   = [  3,   9,  16,  22,  28,  34,  40,  46,  ...
0522         52,  58,  64, 69,  75,  80,  85,  90,  95, 100, 104, ...
0523         108, 112, 116, 120, 123, 126, 129, 132, 134, 136, 138, ...
0524         140, 141, 142, 143, 144, 144, 144, 144, 143, 142, 141, ...
0525         140, 138, 136, 134, 132, 129, 126, 123, 120, 116, 112, ...
0526         108, 104, 100, 95,  90,  85,  80,  75,  69,  64,  58,  ...
0527         52,  46,  40, 34,  28,  22,  16,   9,   3];
0528     
0529     jj = 1; ii = 1;
0530     ngrid = sum(table.lonlist);
0531     table.lon = zeros(1,ngrid);
0532     table.lat = zeros(1,ngrid);
0533     for k = 1:length(table.latarray)
0534         for l = 1:table.lonlist(k)
0535             table.lon(ii) = (l-1)*360.0/table.lonlist(k)+(360.0/table.lonlist(k))/2.0;
0536             table.lat(ii) = table.latarray(jj);
0537             ii = ii+1;
0538         end
0539         jj = jj+1;
0540     end
0541 end
0542 
0543 end
0544 
0545 %% FIELDS NOT READ SINCE THEY ARE INTERNAL
0546 %
0547 % -------
0548 % DX
0549 % -------
0550 
0551 %xpos(sNp:eNp)= geo(:,3); %X-positions for NPIX pixels (1-480)
0552 %ypos(sNp:eNp)= geo(:,4); %Y-positions for NPIX pixels (1-550)
0553 %*BX shore flag (0-1)
0554 %BXSHOR(pix)  = bitshift(bitand(S1(1),2^1,'uint8'),-1);
0555 %*Time/space test result (0-3)
0556 %TIMSPA(pix)  = bitshift(bitand(S1(1),2^6+2^7,'uint8'),-6);
0557 %*IR clear sky composite logic code (0-24)
0558 %ICSLOG(pix)  = bitand(S1(2),2^0+2^1+2^2+2^3+2^4,'uint8');
0559 %*First IR threshold result (0-5)
0560 %BXITHR(pix)  = bitshift(bitand(S1(2),2^5+2^6+2^7,'uint8'),-5);
0561 %* VIS clear sky composite logic code (0-14)
0562 %VCSLOG(pix) = bitshift(bitand(S2(1),sum(2.^(1:4)),'uint8'),-1);
0563 %* First VIS threshold result (0-5)
0564 %BXVTHR(pix) = bitshift(bitand(S2(1),sum(2.^(5:7)),'uint8'),-5);
0565 %* First VIS clear sky radiance (0-254 counts)
0566 %BXVCSR(pix) = S2(5);
0567 %*IR retrieval code (0-12)
0568 %IRET(pix)     = bitand(S3(2),2^0+2^1+2^2+2^3,'uint8');
0569 %*IR clear sky composite retrieval code (0-12)
0570 %ICSRET(pix)   = bitshift(bitand(S3(2),2^4+2^5+2^6+2^7,'uint8'),-4);
0571 %*VIS retrieval code (0-14)
0572 %VRET(pix)   = bitand(S4(1),2^0+2^1+2^2+2^3,'uint8');
0573 %*VIS clear sky composite retrieval code (0-14)
0574 %VCSRET(pix) = bitshift(bitand(S4(1),2^4+2^5+2^6+2^7,'uint8'),-4);

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