Home > atmlab > sensors > cloudsat_read.m

cloudsat_read

PURPOSE ^

CLOUDSAT_READ Reads CloudSat HDF data files.

SYNOPSIS ^

function P = cloudsat_read(filename,fieldlist,choice1,c1arg,choice2,c2arg)

DESCRIPTION ^

 CLOUDSAT_READ   Reads CloudSat HDF data files.

    The function reads data from single files. The data are provided as a
    single structure having field names taken from the HDF data. For
    example, 
       P=cloudsat_read(filename,{'Latitude','Longitude','Sigma_Zero'})
    gives
       P =
            Latitude: [37081x1 double]
           Longitude: [37081x1 double]
          Sigma_Zero: [37081x1 double]

    Some field names are renamed as the corresponding name in HDF can not
    be used in matlab. For example, 'Sigma-Zero' is renamed to
    'Sigma_Zero'.

    Data for field name renaming and unit conversions are hard-coded. The
    following file types are handled so far:
       2B-GEOPROF

    Data are scaled to consider scaling factors and "non-standard"
    units. The data are returned using the following units:
      g, m, dB and deg

    The data are returned as double precision variables, sorted in such
    way that each row corresponds to a position. For example, the cloud 
    mask at position i is: P.CPR_Cloud_mask(i,:)

    The fields existing in a file are obtained as
       fieldnames = cloudsat_read(filename,[],'fields');

    The reading can be restricted in several ways.  To only read data for
    the positions with index *ind*:
       P = cloudsat_read(filename,[],'index',ind);
    To only read data inside latitudes [-30,30]:
       P = cloudsat_read(filename,[],'lat',[-30 30]);
    To only read data inside longitudes [0,180]:
       P = cloudsat_read(filename,[],'lon',[0 180]);
    Note that the range [-180,180] is used for longitudes. To only read 
    data inside leg 2:
       P = cloudsat_read(filename,[],'leg',2);
    The orbit is divided into legs, where the start and end of the orbit
    an latitude turning points are taken as limits of the legs. An orbit
    has accordingly three legs, where e.g. leg 1 extends from the orbit 
    start to the turning point in the South pole area.
       The latitude, longitude and leg options can be combined in pairs,
    such as:
       P = cloudsat_read(filename,[],'leg',2,'lat',[-30 30]);

   The reading can be restricted to a rectangular latitude/longitude area.
   To read "tropical" data for the eastern hemisphere:
        fieldnames = cloudsat_read(filename,[],'lat',[-30 30],'lon',[0 180]);
   Longitude limits can be left out. Default for longitude limits is 
   [-180 180]. 

    The function is based on CPR_L1B_rdr2 (unknown author). Send bug
    reports, suggestions or improvements to patrick.eriksson@chalmers.se

 FORMAT   P = cloudsat_read(filename[,fieldlist,choice1,c1arg,choice2,c2arg])
        
 OUT   P
 IN    filename    Name of file to read.
 OPT   fieldlist   Name of fields to read, as a cell array of strings.
                   Default is [], which results in that all fields are
                   read.
       choice1     Optional choice 1. See above.
       c1arg       Arguments for optional choice 1. See above.
       choice2     Optional choice 2. See above.
       c2arg       Arguments for optional choice 2. See above.

 INFO: If you want to instead read the cloudsat data as it is (i.e. no scaling,
       converting, etc.), use read_cloudsat_hdf directly. The function also
       outputs the variable attributes of the fields you wish to retrieve as a
       second output argument 

 $Id: cloudsat_read.m 7716 2012-07-09 09:36:13Z gerrit $
 2007-10-24   Created by Patrick Eriksson.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

cloudsat_read.m

SOURCE CODE ^

0001 % CLOUDSAT_READ   Reads CloudSat HDF data files.
0002 %
0003 %    The function reads data from single files. The data are provided as a
0004 %    single structure having field names taken from the HDF data. For
0005 %    example,
0006 %       P=cloudsat_read(filename,{'Latitude','Longitude','Sigma_Zero'})
0007 %    gives
0008 %       P =
0009 %            Latitude: [37081x1 double]
0010 %           Longitude: [37081x1 double]
0011 %          Sigma_Zero: [37081x1 double]
0012 %
0013 %    Some field names are renamed as the corresponding name in HDF can not
0014 %    be used in matlab. For example, 'Sigma-Zero' is renamed to
0015 %    'Sigma_Zero'.
0016 %
0017 %    Data for field name renaming and unit conversions are hard-coded. The
0018 %    following file types are handled so far:
0019 %       2B-GEOPROF
0020 %
0021 %    Data are scaled to consider scaling factors and "non-standard"
0022 %    units. The data are returned using the following units:
0023 %      g, m, dB and deg
0024 %
0025 %    The data are returned as double precision variables, sorted in such
0026 %    way that each row corresponds to a position. For example, the cloud
0027 %    mask at position i is: P.CPR_Cloud_mask(i,:)
0028 %
0029 %    The fields existing in a file are obtained as
0030 %       fieldnames = cloudsat_read(filename,[],'fields');
0031 %
0032 %    The reading can be restricted in several ways.  To only read data for
0033 %    the positions with index *ind*:
0034 %       P = cloudsat_read(filename,[],'index',ind);
0035 %    To only read data inside latitudes [-30,30]:
0036 %       P = cloudsat_read(filename,[],'lat',[-30 30]);
0037 %    To only read data inside longitudes [0,180]:
0038 %       P = cloudsat_read(filename,[],'lon',[0 180]);
0039 %    Note that the range [-180,180] is used for longitudes. To only read
0040 %    data inside leg 2:
0041 %       P = cloudsat_read(filename,[],'leg',2);
0042 %    The orbit is divided into legs, where the start and end of the orbit
0043 %    an latitude turning points are taken as limits of the legs. An orbit
0044 %    has accordingly three legs, where e.g. leg 1 extends from the orbit
0045 %    start to the turning point in the South pole area.
0046 %       The latitude, longitude and leg options can be combined in pairs,
0047 %    such as:
0048 %       P = cloudsat_read(filename,[],'leg',2,'lat',[-30 30]);
0049 %
0050 %   The reading can be restricted to a rectangular latitude/longitude area.
0051 %   To read "tropical" data for the eastern hemisphere:
0052 %        fieldnames = cloudsat_read(filename,[],'lat',[-30 30],'lon',[0 180]);
0053 %   Longitude limits can be left out. Default for longitude limits is
0054 %   [-180 180].
0055 %
0056 %    The function is based on CPR_L1B_rdr2 (unknown author). Send bug
0057 %    reports, suggestions or improvements to patrick.eriksson@chalmers.se
0058 %
0059 % FORMAT   P = cloudsat_read(filename[,fieldlist,choice1,c1arg,choice2,c2arg])
0060 %
0061 % OUT   P
0062 % IN    filename    Name of file to read.
0063 % OPT   fieldlist   Name of fields to read, as a cell array of strings.
0064 %                   Default is [], which results in that all fields are
0065 %                   read.
0066 %       choice1     Optional choice 1. See above.
0067 %       c1arg       Arguments for optional choice 1. See above.
0068 %       choice2     Optional choice 2. See above.
0069 %       c2arg       Arguments for optional choice 2. See above.
0070 %
0071 % INFO: If you want to instead read the cloudsat data as it is (i.e. no scaling,
0072 %       converting, etc.), use read_cloudsat_hdf directly. The function also
0073 %       outputs the variable attributes of the fields you wish to retrieve as a
0074 %       second output argument
0075 %
0076 % $Id: cloudsat_read.m 7716 2012-07-09 09:36:13Z gerrit $
0077 % 2007-10-24   Created by Patrick Eriksson.
0078 
0079 
0080 function P = cloudsat_read(filename,fieldlist,choice1,c1arg,choice2,c2arg)
0081   
0082 %- Check input
0083 %
0084 if ~ischar( filename )
0085   error( 'Input argument *filename* must be a string' );
0086 end
0087 %
0088 if ~exist( filename, 'file' );
0089   P = 'File not found.';
0090   return
0091 end
0092 %
0093 %
0094 if nargin < 2
0095   fieldlist = [];
0096 end
0097 %
0098 if ~( isempty( fieldlist )  |  iscellstr( fieldlist ) )
0099   error( ...
0100      'Input argument *fieldlist* must be empty or a cell array of strings.' );
0101 end
0102 %
0103 if nargin < 3
0104   choice1 = [];
0105 elseif ~ischar( choice1 )
0106   error( 'Input argument *choice1* must be a string' );
0107 else
0108   choice1 = lower( choice1 );
0109   %
0110   if ~( strcmp(choice1,'fields')  |  strcmp(choice1,'index')  |  ...
0111         strcmp(choice1,'lat')     |  strcmp(choice1,'lon')  |  ...
0112         strcmp(choice1,'leg') )
0113     error( ['Valid options for *choice1* are: ''fields'', ''index'', ',...
0114             '''lat'', ''lon'' and ''leg''.'] );
0115   end
0116 end
0117 %
0118 if nargin < 5
0119   choice2 = [];
0120 elseif ~ischar( choice2 )
0121   error( 'Input argument *choice2* must be a string' );
0122 else
0123   choice2 = lower( choice2 );
0124   %
0125   if ~( strcmp(choice2,'lat')     |  strcmp(choice2,'lon')  |  ...
0126         strcmp(choice2,'leg') )
0127     error( ['Valid options for *choice1* are: ''lat'', ''lon'', ',...
0128             'and ''leg''.'] );
0129   end
0130 end
0131 
0132 
0133 %- Get information about the file
0134 %
0135 fileinfo = hdfinfo( filename, 'eos' );
0136 %
0137 allfields = { fileinfo.Swath.GeolocationFields.Name, ...
0138              fileinfo.Swath.DataFields.Name };
0139 
0140 %- Handle empty *fieldlist* and choice1=='fields* option
0141 %
0142 if isempty( fieldlist )  |  strcmp( choice1, 'fields' )
0143 
0144   if strcmp( choice1, 'fields' )
0145     P = rename_fields( allfields, 1 );
0146     return
0147   end
0148 
0149   hdfnames  = allfields;
0150   fieldlist = rename_fields( hdfnames, 1 );
0151 
0152 else
0153   hdfnames = rename_fields( fieldlist );
0154 end
0155 
0156 
0157 %- Determine index for data points to keep
0158 %
0159 iout = [];
0160 %
0161 if ~isempty( choice1 )
0162   
0163   %- Index pre-defined
0164   if strcmp( choice1, 'index' )
0165     iout = c1arg;
0166     clear c1arg;
0167   
0168   %- Filtering according to leg, lat and lon
0169   else
0170 
0171     P = cloudsat_read( filename, {'Latitude','Longitude'} );
0172     
0173     %- leg
0174     if strcmp( choice1, 'leg' )  |  strcmp( choice2, 'leg' )
0175       if strcmp( choice1, 'leg' )
0176         if ~( isnumeric(c1arg) & isvector(c1arg) & length(c1arg)==1 )
0177           error( '*c1arg* must for ''leg'' option be a vector of length 1.' );
0178         end
0179         leg = c1arg;
0180       else
0181         if ~( isnumeric(c2arg) & isvector(c2arg) & length(c2arg)==1 )
0182           error( '*c2arg* must for ''leg'' option be a vector of length 1.' );
0183         end
0184         leg = c2arg;
0185       end
0186       %
0187       if leg == 1
0188         i1     = 1;
0189         [u,i2] = min( P.Latitude );
0190       elseif leg == 2
0191         [u,i1] = min( P.Latitude ); i1 = i1 + 1;
0192         [u,i2] = max( P.Latitude );
0193       elseif leg == 3
0194         [u,i1] = max( P.Latitude ); i1 = i1 + 1;
0195         i2     = length( P.Latitude );
0196       else
0197         error( 'Possible choices for leg index are 1, 2 and 3.' );
0198       end
0199       inleg        = zeros( length(P.Latitude), 1 );
0200       inleg(i1:i2) = 1;
0201     else
0202       inleg = ones( length(P.Latitude), 1 );
0203     end   
0204   
0205     %- latitude
0206     if strcmp( choice1, 'lat' )  |  strcmp( choice2, 'lat' )
0207       if strcmp( choice1, 'lat' )
0208         if ~( isnumeric(c1arg) & isvector(c1arg) & length(c1arg)==2 )
0209           error( '*c1arg* must for ''lat'' option be a vector of length 2.' );
0210         end
0211         latlims = c1arg;
0212       else  
0213         if ~( isnumeric(c2arg) & isvector(c2arg) & length(c2arg)==2 )
0214           error( '*c2arg* must for ''lat'' option be a vector of length 2.' );
0215         end
0216         latlims = c2arg;
0217       end
0218     else
0219       latlims = [-90 90 ];
0220     end   
0221 
0222     %- longitude
0223     if strcmp( choice1, 'lon' )  |  strcmp( choice2, 'lon' )
0224       if strcmp( choice1, 'lon' )
0225         if ~( isnumeric(c1arg) & isvector(c1arg) & length(c1arg)==2 )
0226           error( '*c1arg* must for ''lon'' option be a vector of length 2.' );
0227         end
0228         lonlims = c1arg;
0229       else  
0230         if ~( isnumeric(c2arg) & isvector(c2arg) & length(c2arg)==2 )
0231           error( '*c2arg* must for ''lon'' option be a vector of length 2.' );
0232         end
0233         lonlims = c2arg;
0234       end
0235     else
0236       lonlims = [-180 180 ];
0237     end   
0238 
0239     iout = inleg & ...
0240         P.Latitude  >= latlims(1) & P.Latitude  <= latlims(2) & ...
0241         P.Longitude >= lonlims(1) & P.Longitude <= lonlims(2);
0242     
0243     if isempty(iout)
0244       P = 'No data matching reading restrictions.';
0245       return
0246     end
0247   end
0248 end
0249 
0250 
0251 
0252 %- Read data
0253 %
0254 % The reading routines cause a lot of warnings (at least in Matlab 7). Turn
0255 % off these warnings temporariliy.
0256 %
0257 % Some files appear to be corrupted and cause error. Use try/catch.
0258 %
0259 warning( 'off', 'MATLAB:hdfwarn:generic' );
0260 P = [];
0261 
0262 %
0263 try
0264     P = read_cloudsat_hdf(fileinfo.Filename,fieldlist);
0265     for F = fieldlist
0266         % convention to make everything double
0267         if ~ischar(P.(F{1}))
0268             P.(F{1}) = double(P.(F{1}));
0269         end
0270         % convention to make all vectors column vectors
0271         if isvector(P.(F{1}))
0272             P.(F{1}) = P.(F{1})(:);
0273         end
0274         % apply conditions
0275         if ~isempty(iout) && ( isvector(P.(F{1})) || ismatrix(P.(F{1})) )
0276             P.(F{1}) = P.(F{1})(iout,:);
0277         end
0278     end
0279 catch ME
0280     P = sprintf('Error while reading file.: %s\n',ME.message);
0281     warning('MATLAB:hdfwarn:generic','%s',P);
0282 end
0283 
0284 if ischar( P )
0285   return
0286 end
0287 
0288 
0289 %- Factor and offset
0290 %
0291 attributes = strvcat(fileinfo.Swath.Attributes.Name);
0292 for it=1:length(fieldlist)
0293   ifactor = strmatch([hdfnames{it} '.factor'],attributes);
0294   ioffset = strmatch([hdfnames{it} '.offset'],attributes);
0295   iunit   = strmatch([hdfnames{it} '.units'],attributes);
0296   if ~isempty(iunit)
0297     unit=fileinfo.Swath.Attributes(iunit).Value;
0298     unit=unitconv(unit,hdfnames{it});
0299   else
0300     unit=1;
0301   end
0302 
0303   if ~isempty(ifactor) & ~isempty(ioffset) 
0304      factor = double(fileinfo.Swath.Attributes(ifactor).Value(1));
0305      offset = double(fileinfo.Swath.Attributes(ioffset).Value(1));
0306      P.(fieldlist{it}) = (P.(fieldlist{it})/factor - offset)*unit;
0307   end
0308 end
0309 
0310   
0311 return
0312 
0313 
0314 
0315 
0316 %--- Sub functions ---------------------------------------------------------
0317 
0318 %--- Conversion of field names
0319 %
0320 % Renaming information stored in internal variable *table*. Default is
0321 % conversion from matlab to HDF names.
0322 %
0323 % OUT   names       Renamed field names
0324 % IN    names       Original field names
0325 %       backwards   Set to true for conversion to matlab names. Default
0326 %                   is false.
0327 %
0328 function names = rename_fields( names, backwards )
0329 
0330   if nargin < 2, backwards = 0; end  
0331 
0332   %- Renaming table
0333   %
0334   % HDF name in column 1, matlab name in column 2
0335   %
0336   table{1} = { 'Sigma-Zero', 'Sigma_Zero' };
0337   
0338   if backwards
0339     iout = 2;
0340     iin  = 1;
0341   else
0342     iout = 1;
0343     iin  = 2;
0344   end
0345 
0346   for it = 1 : length(table)
0347     ind = find( strcmp( names, table{it}{iin} ) );
0348     if ~isempty( ind )
0349       names{ind} = table{it}{iout};
0350     end
0351   end
0352 
0353 return
0354 
0355 
0356 %--- Unit conversion
0357 %
0358 % Conversion information stored in internal variable *table*.
0359 %
0360 function u=unitconv(unit,field)
0361 
0362 % Make sure that *unit* is a row vector (can be column for old HDF versions)
0363 %
0364 unit = vec2row( unit );
0365   
0366 if strcmp(unit,['mg m^{-3}']) |  strcmp(unit,['mg/m^3'])
0367    u=1e-3;
0368 elseif strcmp(unit,['um'])
0369    u=1e-6;
0370 elseif strcmp(unit,['degrees'])
0371    u=1;
0372 elseif strcmp(unit,['m'])  |  strcmp(unit,['/m'])
0373    u=1;
0374 elseif strcmp(unit,['meters'])
0375    u=1;
0376 elseif strcmp(unit,['dBZe'])  |  strcmp(unit,['dbz'])  |  strcmp(unit,['dBZ'])
0377    u=1;
0378 elseif strcmp(unit,['dB*100'])
0379    u=1e-2;
0380 elseif strcmp(unit,['km'])
0381    u=1e3;
0382 elseif strcmp(unit,['1/km'])  |  strcmp(unit,['/km'])
0383    u=1e-3;
0384 elseif strcmp(unit,['seconds'])
0385    u=1;
0386 elseif strcmp(unit,['--'])  |  strcmp(unit,['none'])  |  strcmp(unit,['None'])
0387    u=1;
0388 elseif strcmp(unit,['%'])
0389    u=1;
0390 elseif strcmp(unit,['g m^{-2}'])  |  strcmp(unit,['g/m^2'])
0391    u=1;
0392 elseif strcmp(unit,['K'])
0393    u=1;
0394 elseif strcmp(unit,['L^{-1}'])
0395    u=1e3;
0396 elseif strcmp(unit,['cm^{-3}'])
0397    u=1e6;
0398 elseif strcmp(unit,['Pa'])
0399    u=1;
0400 elseif strcmp(unit,['kg/kg']) |  strcmp(unit,['kg kg**-1']) 
0401    u=1;
0402 else 
0403    u=1;
0404    warning( 'cloudsat_read:unknownUnit', ...
0405         ['Could not identify field type unit in subfunction unitconv. ',...
0406          'No unit conversion performed for ',field,'. ',...
0407          'Unit is in ',unit,'.'] );
0408 end
0409 
0410 return

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