Home > atmlab > datasets > cpr_postprocessor.m

cpr_postprocessor

PURPOSE ^

% CPR_POSTPROCESSOR

SYNOPSIS ^

function S = cpr_postprocessor(self,S,flds)

DESCRIPTION ^

% CPR_POSTPROCESSOR

 PURPOSE: To create "pseudo" fields that are not in the original cloudsat
          data by post-processing fields that are in the original data


 IN
     self = dataset
     S = Data structure to add pseudo field to (from self.reader)
     flds = {'names','of','pseudo','fields'} (fields)

 OUT
     S = Data structure + pseudo field/s

 NOTE: See also help Satdataset/pseudo_fields

 $Id: cpr_postprocessor.m 8433 2013-05-22 13:29:28Z seliasson $
 Salomon Eliasson

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

cpr_postprocessor.m

SOURCE CODE ^

0001 function S = cpr_postprocessor(self,S,flds)
0002 %% CPR_POSTPROCESSOR
0003 %
0004 % PURPOSE: To create "pseudo" fields that are not in the original cloudsat
0005 %          data by post-processing fields that are in the original data
0006 %
0007 %
0008 % IN
0009 %     self = dataset
0010 %     S = Data structure to add pseudo field to (from self.reader)
0011 %     flds = {'names','of','pseudo','fields'} (fields)
0012 %
0013 % OUT
0014 %     S = Data structure + pseudo field/s
0015 %
0016 % NOTE: See also help Satdataset/pseudo_fields
0017 %
0018 % $Id: cpr_postprocessor.m 8433 2013-05-22 13:29:28Z seliasson $
0019 % Salomon Eliasson
0020 
0021 % --------------
0022 % Do operations
0023 % --------------
0024 
0025 
0026 for F =flds
0027     field = F{1};
0028     % get the missing value
0029     mv = self.pseudo_fields.(field).atts.missing_value;
0030     switch field
0031         case {'Cloud_Types','Cloud_Types_multiLayer'}
0032 
0033             % preallocate
0034             x.cloudy = false(size(S.cloud_scenario));
0035             x.Ci = x.cloudy; x.As = x.cloudy; x.Ac = x.cloudy;
0036             x.St = x.cloudy; x.Sc = x.cloudy; x.Cu = x.cloudy;
0037             x.Ns = x.cloudy; x.DC = x.cloudy;
0038             
0039             % -----------------------
0040             % extract clouds
0041             % -----------------------
0042             CS = sum(2.^(1:4)); %scenario mask
0043             Nprof = size(S.cloud_scenario,1);
0044             ph=mv*ones(1,10);%placeholder
0045             for i = 1:Nprof
0046                 determined = bitand(typecast(S.cloud_scenario(i,:),'uint16'),2^0)==2^0;
0047                 cloudScenario = bitand(typecast(S.cloud_scenario(i,:),'uint16'),CS);
0048                 x.cloudy(i,:) = determined & cloudScenario > 0;
0049                 x.Ci(i,:) = determined & cloudScenario == 2;
0050                 x.As(i,:) = determined & cloudScenario == 2^2;
0051                 x.Ac(i,:) = determined & cloudScenario == 2^1+2^2;
0052                 x.St(i,:) = determined & cloudScenario == 2^3;
0053                 x.Sc(i,:) = determined & cloudScenario == 2^1+2^3;
0054                 x.Cu(i,:) = determined & cloudScenario == 2^2+2^3;
0055                 x.Ns(i,:) = determined & cloudScenario == 2^1+2^2+2^3;
0056                 x.DC(i,:) = determined & cloudScenario == 2^4;
0057                 if strcmp(field,'Cloud_Types_multiLayer')
0058                     singleThings = (1:8).*(any([x.Ci(i,:);x.As(i,:);x.Ac(i,:);x.St(i,:);x.Sc(i,:);x.Cu(i,:);x.Ns(i,:);x.DC(i,:)],2))';
0059                     singleThings=singleThings(singleThings~=0); if isempty(singleThings), singleThings=0;end
0060                     S.Cloud_Types_multiLayer(i,:) = ph;
0061                     S.Cloud_Types_multiLayer(i,1:length(singleThings)) = singleThings;
0062                 end
0063 
0064                 
0065             end            
0066             % -------------------------
0067             % Make logical true if all
0068             % clouds in profile are the same type
0069             % ------------------------
0070             
0071             cloudfree = ~any(x.cloudy,2);
0072             cloudTypes = {'Ci','As','Ac','St','Sc','Cu','Ns','DC'};
0073             for i = 1:length(cloudTypes)
0074                 CT = cloudTypes{i};
0075                 % in profile: sum(cloudtype) = sum(cloudy vertical bins)
0076                 singleType.(CT) = sum(x.(CT),2)==sum(x.cloudy,2);
0077                 singleType.(CT)(cloudfree) = false;
0078             end
0079             
0080             % mixed clouds is where the profile is not cloudfree nor does it contain only one cloud type (=9)
0081             S.Cloud_Types = (length(cloudTypes)+1)*ones(size(singleType.Ci),'int8');
0082             for i = 1:length(cloudTypes)
0083                 S.Cloud_Types(singleType.(cloudTypes{i})) = i;
0084             end
0085             S.Cloud_Types(cloudfree) = 0;
0086         case {'Cloud_Types_Lidar','Cloud_Types_Lidar_multiLayer'}
0087             
0088             % The main data
0089             clouds = S.CloudLayerType;
0090             
0091             % cloud types come from CloudLayerType, but the pseudo_field: Cloud_Types_Lidar is a
0092             % column integrated quantity. Initialize.
0093             if strcmp(field,'Cloud_Types_Lidar_multiLayer')
0094                 S.(field)=ones(size(clouds),'int8')*mv;
0095                 logtext(atmlab('OUT'),'Filtering cloudtypes according to quality and assigning CT=0 if cloud free\n')
0096             elseif strcmp(field,'Cloud_Types_Lidar')
0097                 S.(field)=ones(size(clouds,1),1,'int8')*mv;
0098                 logtext(atmlab('OUT'),'Assigning cloudtype to each profile\n')
0099             end
0100             
0101             % I need to loop over every point in order to use the index of the cloud
0102             % fraction and the clouds together
0103             
0104             for i = 1:size(clouds,1)
0105                 
0106                 % ---------------
0107                 % Assign clouds
0108                 if strcmp(field,'Cloud_Types_Lidar')
0109                     cloudColumn = clouds(i,clouds(i,:)~=0);
0110                     if isempty(cloudColumn)
0111                         % undetermined cloud type (is = 0 in the original [clear, or
0112                         % undetermined]). I set "clear" explicitly using the
0113                         % CloudFraction field
0114                         S.(field)(i) = 10;
0115                     elseif length(unique(cloudColumn))==1 % this only works since, if it is clear, all values are 0
0116                         % column with a single cloud type
0117                         S.(field)(i) = unique(cloudColumn);
0118                     else
0119                         % mixed cloud column
0120                         S.(field)(i) = 9;
0121                     end
0122                 elseif strcmp(field,'Cloud_Types_Lidar_multiLayer')
0123                     S.(field)(i,:) = clouds(i,:);
0124                     S.(field)(i,clouds(i,:)==0)=10; % reassigning undetermined to 10. Will put zeros if CloudFraction says zero (see below)
0125                 end
0126                 
0127                 
0128                 % -------------------
0129                 % MASK special cases
0130                 % -------------------
0131                 
0132                 % ------------
0133                 % Bad quality
0134                 % Assign a missing value to a profile if any of the cloud boxes have CloudTypeQuality < 0.5
0135                 bad = S.CloudTypeQuality(i,:)<.5;
0136                 if strcmp(field,'Cloud_Types_Lidar')
0137                     if any(bad)
0138                         S.(field)(i) = mv;
0139                     end
0140                 elseif strcmp(field,'Cloud_Types_Lidar_multiLayer')
0141                     S.(field)(i,bad) = mv;
0142                 end
0143                 
0144                 % -----------------
0145                 % Cloud Fraction
0146                 % Use the lidar cloud fraction to find the truly CLOUD FREE
0147                 % measurements (-99), and mask partly cloudy cloudsat
0148                 % footprints/ volumes
0149                 cloudy = S.CloudFraction(i,:);
0150                 if strcmp(field,'Cloud_Types_Lidar')
0151                     if any(cloudy>0&cloudy<1)
0152                         % not used for single profile
0153                         S.(field)(i) = mv;
0154                     elseif all(cloudy==-99)
0155                         % cloud free
0156                         S.(field)(i) = 0;
0157                     end
0158                 elseif strcmp(field,'Cloud_Types_Lidar_multiLayer')
0159                     S.(field)(i,cloudy>0&cloudy<1) = mv;
0160                     S.(field)(i,cloudy==-99)=0;
0161                 end
0162                 
0163                 % ------------------
0164                 % Final sweep for multi clouds
0165                 % ------------------
0166                 % only count clouds covering consecutive bins once.
0167                 
0168                 if strcmp(field,'Cloud_Types_Lidar_multiLayer')
0169                     singleThings = S.(field)(i,diff(S.(field)(i,:))~=0);
0170                     if isempty(singleThings) % if not every box has a cloud in it
0171                         singleThings = 0; % i.e, cloud free
0172                     end
0173                     S.(field)(i,:) = mv*ones(1,10); % put missing values as fillers
0174                     S.(field)(i,1:length(singleThings)) = singleThings;
0175                 end
0176             end
0177         otherwise
0178             error(['atmlab:' mfilename],'field: "%s" is not listed here',field)
0179             
0180     end
0181 end
0182 
0183 end

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