Home > atmlab > arts > scenegen > asg_cloudsat_dBZe.m

asg_cloudsat_dBZe

PURPOSE ^

ASG_CLOUDSAT_DBZE reads in cloudsat data on gformat

SYNOPSIS ^

function [G]=asg_cloudsat_dBZe(G,data,leg,latitude,varargin)

DESCRIPTION ^

 ASG_CLOUDSAT_DBZE reads in cloudsat data on gformat

          G(1) will holds radar reflectivity data.
          If it is desired G(1+i) will holds
          auxilirary data

 FORMAT   G=asg_cloudsat_dbze( G, data, leg, latitude, z_grid, aux ,dBZemin)
 
 IN      
          G        original gformat data
                   can be empty             
          data     name of a cloudsat file
          leg      options 1,2, or 3 see cloudsat_read
                   for definitions
          latitude latitude region to be read in [lat1 lat2] 
 OPT      z_grid   the data will be binned on this grid
                   default is the mean of the cloudsat
                   array of grids
          aux      auxilirary data,options can for instance be
                   {'DEM_elevation','Longitude','TAI_start','Profile_time'}
                   default are none
          dBZemin  Lowest value of dBZe to be considered
                   default is -50  

example usage on how to read in data from
a cloudsat 2B-GEOPROF file

datadir='/home/bengt/CIWSIR/Dataset_gen/CloudSAT_data/Data/';
data=[datadir,'2006166131201_00702_CS_2B-GEOPROF_GRANULE_P_R04_E00.hdf'];
z_grid=[0:1e3:20e3]'; 
aux={'DEM_elevation','Longitude','TAI_start','Profile_time'};
leg=2;
latitude=[0 5];
[G,D]=gfin_cloudsat_dBZe([],data,aux,[],leg,latitude,[]);

2007-11-05 created by Bengt Rydberg

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

asg_cloudsat_dBZe.m

SOURCE CODE ^

0001 % ASG_CLOUDSAT_DBZE reads in cloudsat data on gformat
0002 %
0003 %          G(1) will holds radar reflectivity data.
0004 %          If it is desired G(1+i) will holds
0005 %          auxilirary data
0006 %
0007 % FORMAT   G=asg_cloudsat_dbze( G, data, leg, latitude, z_grid, aux ,dBZemin)
0008 %
0009 % IN
0010 %          G        original gformat data
0011 %                   can be empty
0012 %          data     name of a cloudsat file
0013 %          leg      options 1,2, or 3 see cloudsat_read
0014 %                   for definitions
0015 %          latitude latitude region to be read in [lat1 lat2]
0016 % OPT      z_grid   the data will be binned on this grid
0017 %                   default is the mean of the cloudsat
0018 %                   array of grids
0019 %          aux      auxilirary data,options can for instance be
0020 %                   {'DEM_elevation','Longitude','TAI_start','Profile_time'}
0021 %                   default are none
0022 %          dBZemin  Lowest value of dBZe to be considered
0023 %                   default is -50
0024 %
0025 %example usage on how to read in data from
0026 %a cloudsat 2B-GEOPROF file
0027 %
0028 %datadir='/home/bengt/CIWSIR/Dataset_gen/CloudSAT_data/Data/';
0029 %data=[datadir,'2006166131201_00702_CS_2B-GEOPROF_GRANULE_P_R04_E00.hdf'];
0030 %z_grid=[0:1e3:20e3]';
0031 %aux={'DEM_elevation','Longitude','TAI_start','Profile_time'};
0032 %leg=2;
0033 %latitude=[0 5];
0034 %[G,D]=gfin_cloudsat_dBZe([],data,aux,[],leg,latitude,[]);
0035 %
0036 %2007-11-05 created by Bengt Rydberg
0037 
0038 function [G]=asg_cloudsat_dBZe(G,data,leg,latitude,varargin)
0039 
0040 [z_grid,aux,dBZemin]=optargs(varargin,{[],[],-50});
0041 
0042 fields= {'Latitude','Height','Radar_Reflectivity',...
0043                     'CPR_Cloud_mask','Longitude'};
0044 
0045 name = 'Radar_Reflectivity';
0046 type = 'dBZe';
0047 dims = [1 2];
0048 
0049 D.DIM=3;
0050 D.GRID1_NAME='altitude';
0051 D.GRID2_NAME='latitude';
0052 D.GRID3_NAME='longitude';
0053 
0054 
0055 for id = 1 : D.DIM
0056     gname = lower( D.(sprintf('GRID%d_NAME',id)) );
0057     if strncmp( gname, 'altitude', 8 )
0058       dims(1) = id;
0059     elseif strncmp( gname, 'latitude', 8 )
0060       dims(2) = id;
0061     elseif strncmp( gname, 'longitude', 8 )
0062       dims(3) = id;
0063     end
0064     
0065 end
0066 
0067 if any( dims == 0 )
0068     error( [ 'Setting of *dims* failed. Altitude, latitude or ',...
0069              'longitude was nof found among the grid names.' ] );
0070 end
0071 
0072 if ~isempty(aux)
0073    for i=1:length(aux)
0074        if isempty(find(strcmp(aux(i),fields)))
0075           fields{end+1}=aux{i};
0076        end
0077    end
0078 end
0079 
0080 %read data
0081 P = cloudsat_read(data,fields,'leg',leg,'lat',latitude);
0082 
0083 if ischar(P)
0084    error('the selected cloudsat file was not found')
0085 end
0086 
0087 %The data in P is a function of (lat,alt) and the altitude is decreasing,
0088 %but we want is as function of (alt,lat) and increasing altitude,
0089 %therefore we transpose the data and flip it upside down
0090 %note that nothing happens with the latitude vector
0091 %since the flipud will be performed on a row vector for latitude
0092 
0093 for i=1:length(fields)
0094     P.(fields{i}) = flipud(P.(fields{i})');
0095 end
0096 
0097 
0098 % The code was: else exist('z_grid','var')
0099 % Correct that it shall be a "pure" else?yes
0100 if ~exist('z_grid','var') 
0101    P.z_grid=mean(P.Height')';
0102 else 
0103    if ~isempty(z_grid)
0104       P.z_grid=z_grid;
0105    else
0106       P.z_grid=mean(P.Height')';
0107    end
0108 end
0109 
0110 %remove noisy data not to be identified as clouds
0111 cind=zeros(size(P.CPR_Cloud_mask));
0112 cind(find(P.CPR_Cloud_mask>=20))=1;%cloudy data
0113 P.Radar_Reflectivity=P.Radar_Reflectivity.*cind-(cind-1)*dBZemin;
0114 P.Radar_Reflectivity(find(P.Radar_Reflectivity<dBZemin))=dBZemin;
0115 
0116 
0117 
0118 if length(fields)>4
0119    kvec=[1 5:length(fields)];
0120 else
0121   kvec=1;
0122 end
0123 
0124 
0125 for ik=kvec
0126     if ik==1
0127        ind=3;
0128        H=G;
0129        H.DATA_UNIT='dBZe';
0130        i=1;
0131     else
0132       ind=ik;
0133       H(end+1)=H(end);
0134       i=length(H);
0135       H(i).DATA_UNIT='?';
0136     end
0137     H(i).NAME=fields{ind};
0138     H(i).DATA_NAME=fields{ind};
0139     y=P.(fields{ind});
0140     %find out dimensionality of loaded data
0141     S=[size(y)==size(P.Latitude)];
0142     if S(1)==1 & S(2)==1 
0143        %data is a vector
0144        for id = 1 : length(dims)
0145            gname = sprintf( 'GRID%d', dims(id) );
0146            if id == 1
0147           [H(i).(gname)] = [];
0148            elseif id == 2
0149           [H(i).(gname)] = vec2col(P.Latitude)';
0150            elseif id == 3
0151               [H(i).(gname)] = [];
0152            end
0153        end
0154        H(i).DIMS=dims(2);
0155        H(i).DATA=y;
0156     elseif  S(1)==1 & S(2)==0
0157        %data is a scalar
0158        H(i).GRID1=[];
0159        H(i).GRID2=[];
0160        H(i).GRID3=[];
0161        H(i).DIMS=[];
0162        H(i).DATA=y;
0163     elseif  S(1)==0 & S(2)==1
0164        %data is a matrix
0165        y1=zeros(length(P.z_grid),size(y,2));
0166        warning off;
0167        %bin the data on altitude grid
0168        for iy=1:size(y,2)
0169            y1(:,iy)=binning(P.z_grid,P.Height(:,iy),y(:,iy));
0170        end
0171        y=y1;
0172        warning on;
0173        if ik~=1
0174       warning(['the data in fields ',fields{ind},' has been binned',...
0175                   ' to fit on a z_grid']);
0176        end
0177        for id = 1 : length(dims)
0178            gname = sprintf( 'GRID%d', dims(id) );
0179            if id == 1
0180               [H(i).(gname)] = P.z_grid;
0181            elseif id == 2
0182           [H(i).(gname)] = vec2col(P.Latitude)';
0183            elseif id == 3
0184               [H(i).(gname)] = [];
0185            end
0186        end
0187        if dims(1)~=1 | dims(2)~=2 | dims(3)~=3
0188       %data is y(z,lat) change order if necessary
0189       grid_ind=[dims(1) dims(2)];
0190       if grid_ind(1)>grid_ind(2)
0191              y=y';
0192           end
0193        end
0194        H(i).DATA=y;
0195        H(i).DIMS=sort([dims(1) dims(2)]);
0196    else
0197        error('unknown dimensionality of data');
0198    end
0199 
0200    H(i).SOURCE=data;
0201 end
0202  
0203 G=H;
0204 
0205 
0206 G0 = gf_empty( 3 );
0207 field={G(:).NAME};
0208 ind=find(strcmp(field,{'Radar_Reflectivity'}));
0209 H1=G0;
0210 H1=gf_set_fields(H1,'TYPE','atm_data','NAME',G(ind).NAME,...
0211                'SOURCE',G(ind).SOURCE,'DATA_NAME',G(ind).DATA_NAME,...
0212                'DATA_UNIT',G(ind).DATA_UNIT);
0213 H1=gf_set(H1,G(ind).DATA,[{G(ind).GRID1},{G(ind).GRID2'},{}]);     
0214 H1=gf_set_fields(H1,'GRID1_NAME','altitude',...
0215                'GRID1_UNIT','m','GRID2_NAME','latitude',...
0216                'GRID2_UNIT','degree');
0217 
0218 ind=find(strcmp(field,{'Longitude'}));
0219 H1(2)=G0;
0220 H1(2)=gf_set_fields(H1(2),'TYPE','atm_data','NAME',G(ind).NAME,...
0221                'SOURCE',G(ind).SOURCE,'DATA_NAME',G(ind).DATA_NAME,...
0222                'DATA_UNIT','degree');
0223 H1(2)=gf_set(H1(2),G(ind).DATA',[{G(ind).GRID2'}]);     
0224 H1(2)=gf_set_fields(H1(2),'GRID1_NAME','latitude',...
0225                'GRID1_UNIT','degree');
0226 
0227 ind=find(strcmp(field,{'DEM_elevation'}));
0228 H1(3)=G0;
0229 H1(3)=gf_set_fields(H1(3),'TYPE','atm_data','NAME',G(ind).NAME,...
0230                'SOURCE',G(ind).SOURCE,'DATA_NAME',G(ind).DATA_NAME,...
0231                'DATA_UNIT','m');
0232 H1(3)=gf_set(H1(3),G(ind).DATA',[{G(ind).GRID2'}]);     
0233 H1(3)=gf_set_fields(H1(3),'GRID1_NAME','latitude',...
0234                'GRID1_UNIT','degree');
0235 
0236 ind1=find(strcmp(field,{'TAI_start'}));
0237 ind2=find(strcmp(field,{'Profile_time'}));
0238 sec_of_day=60*60*24;
0239 mjd1=date2mjd(1993,1,1,0,0,0)+(G(ind1).DATA+G(ind2).DATA)/sec_of_day;
0240 H1(4)=G0;
0241 H1(4)=gf_set_fields(H1(3),'TYPE','mjd','NAME','mjd',...
0242                'SOURCE',G(ind1).SOURCE,'DATA_NAME','mjd',...
0243                'DATA_UNIT','days');
0244 H1(4)=gf_set(H1(4),mjd1',[{G(ind2).GRID2'}]);     
0245 H1(4)=gf_set_fields(H1(4),'GRID1_NAME','latitude',...
0246                'GRID1_UNIT','degree');
0247 
0248 G=H1;

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