Home > atmlab > demos > qpack2_demo.m

qpack2_demo

PURPOSE ^

QPACK2_DEMO Demonstration of the Qpack2 retrieval system

SYNOPSIS ^

function L2 = qpack2_demo

DESCRIPTION ^

 QPACK2_DEMO   Demonstration of the Qpack2 retrieval system

    The main features of Qpack2 are demonstrated. The example case is airborne
    measurements of ozone at 110.8 GHz. Synthetic measurement data are
    generated internally. See the code and internal comments for details.

    Everything is here put into a single file. For practical retrievals it
    is probably better to put the definitions of Q (together with O?) in a
    separate function (i.e. [Q,O] = q_mycase). A function to import
    measurement data into the "Y format" (see *qp2_y*) is needed. The
    retrieval result is returned in the L2 format produced by *qp2_l2*.

    This script focuses on giving an introduction, indicating different
    retrieval units and retrieval of other variables. See also *qpack2_demo2*.

 FORMAT   L2 = qpack2_demo

 OUT      L2   L2 data output from *qpack2*.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

qpack2_demo.m

SOURCE CODE ^

0001 % QPACK2_DEMO   Demonstration of the Qpack2 retrieval system
0002 %
0003 %    The main features of Qpack2 are demonstrated. The example case is airborne
0004 %    measurements of ozone at 110.8 GHz. Synthetic measurement data are
0005 %    generated internally. See the code and internal comments for details.
0006 %
0007 %    Everything is here put into a single file. For practical retrievals it
0008 %    is probably better to put the definitions of Q (together with O?) in a
0009 %    separate function (i.e. [Q,O] = q_mycase). A function to import
0010 %    measurement data into the "Y format" (see *qp2_y*) is needed. The
0011 %    retrieval result is returned in the L2 format produced by *qp2_l2*.
0012 %
0013 %    This script focuses on giving an introduction, indicating different
0014 %    retrieval units and retrieval of other variables. See also *qpack2_demo2*.
0015 %
0016 % FORMAT   L2 = qpack2_demo
0017 %
0018 % OUT      L2   L2 data output from *qpack2*.
0019 
0020 % 2010-05-12   Created by Patrick Eriksson.
0021 
0022 function L2 = qpack2_demo
0023   
0024 errid = ['atmlab:' mfilename]; 
0025 
0026 %- Qarts settings
0027 %
0028 Q    = q_demo;           % Local file, found below
0029 
0030 
0031 %- Measurement data
0032 %
0033 Y = y_demo( Q );         % Local file, found below
0034 
0035 
0036 %- Check that all frequencies are OK
0037 %
0038 if ~qp2_check_f( Q, Y, 1e3 );
0039   error( errid, ...
0040             'Some mismatch between Q.F_BACKEND and frequencies of spectra.' );
0041 end
0042 
0043 
0044 %- OEM variables
0045 %
0046 O = qp2_l2( Q );         % This ensures that OEM returns the variables needed
0047 %                          to fill the L2 structure, as defined in Q
0048 O.linear = false;
0049 %
0050 if ~O.linear 
0051   O.itermethod = 'GN';
0052   O.stop_dx    = 0.01;
0053   O.maxiter    = 5;
0054 end
0055 
0056 
0057 %- Make inversion
0058 %
0059 L2 = qpack2( Q, O, Y );
0060 
0061 
0062 %- Plot, if no output argument
0063 %
0064 if ~nargout
0065   
0066   % Profiles
0067   figure(1),clf
0068   plot( L2(1).species1_x*1e6, p2z_simple(L2(1).species1_p)/1e3, 'b', ...
0069         L2(2).species1_x*1e6, p2z_simple(L2(2).species1_p)/1e3, 'r', ...
0070         L2(1).species1_xa*1e6, p2z_simple(L2(1).species1_p)/1e3, 'k-' );
0071   xlabel( 'Ozone [VMR]' );
0072   ylabel( 'Approximate altitude [km]' );
0073   axis( [ 0 8 10 90 ] )
0074   legend( 'Retrieval 1', 'Retrieval 2', 'True and a priori' );
0075 
0076   % Spectra
0077   figure(2),clf
0078   h=plot( L2(1).f/1e9, L2(1).y, 'k.', L2(2).f/1e9, L2(2).y, 'k.', ...
0079           L2(1).f/1e9, L2(1).yf, 'b-', L2(2).f/1e9, L2(2).yf, 'r-', ...
0080           L2(1).f/1e9, L2(1).bl, 'b-.', L2(2).f/1e9, L2(2).bl, 'r-.' );
0081   xlabel( 'Frequency [GHz]' );
0082   ylabel( 'Tb [K]' );
0083   axis( [ min(L2(1).f/1e9) max(L2(1).f/1e9) 0 18 ] )
0084   legend( h(3:end), 'Fitted 1', 'Fitted 2', 'Baseline 1', 'Baseline 2' );
0085 end
0086 
0087 return
0088 %-----------------------------------------------------------------------------
0089 %-----------------------------------------------------------------------------
0090 %-----------------------------------------------------------------------------
0091 
0092 
0093 
0094 
0095 
0096 function Q = q_demo
0097  
0098 errid = ['atmlab:' mfilename];
0099 %- Atmlab settings
0100 %
0101 arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
0102 arts_includes     = atmlab( 'ARTS_INCLUDES' );
0103 if isnan( arts_xmldata_path ) 
0104   error( errid,'You need to set ARTS_XMLDATA_PATH to run this exmaple.' );
0105 end
0106 if isnan( arts_includes )
0107   error( erird,'You need to ARTS_INCLUDES to run this example.' );
0108 end                                                      
0109 %
0110 fascod = fullfile( arts_xmldata_path, 'planets', 'Earth', 'Fascod' );
0111 
0112 
0113 %- Init Q
0114 %
0115 Q = qarts;
0116 %
0117 Q.INCLUDES            = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
0118                           fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
0119                           fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
0120                           fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
0121 Q.ATMOSPHERE_DIM      = 1;
0122 Q.STOKES_DIM          = 1;
0123 Q.J_DO                = true;
0124 Q.CLOUDBOX_DO         = false;
0125 
0126 
0127 %= Define agendas
0128 %
0129 % Here we do it by using the predefined agenda templates
0130 %   (found in arts/controlfiles/general/agendas.arts)
0131 % This works only if the pre-defined agenda is names following the pattern:
0132 %   name_of_agenda__(Something)
0133 %
0134 Q.PPATH_AGENDA               = { 'ppath_agenda__FollowSensorLosPath'   };
0135 Q.PPATH_STEP_AGENDA          = { 'ppath_step_agenda__GeometricPath'    };
0136 Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck'  };
0137 Q.IY_SPACE_AGENDA            = { 'iy_space_agenda__CosmicBackground'   };
0138 Q.IY_SURFACE_AGENDA          = { 'iy_surface_agenda__UseSurfaceRtprop' };
0139 Q.IY_MAIN_AGENDA             = { 'iy_main_agenda__Emission'            };
0140 
0141 
0142 %- Radiative transfer
0143 %
0144 Q.IY_UNIT             = 'RJBT'; 
0145 Q.YCALC_WSMS          = { 'yCalc' };
0146 %
0147 Q.PPATH_LMAX          = 250;
0148 
0149 
0150 %- Surface
0151 %
0152 Q.Z_SURFACE           = 1e3;          % Just a dummy value. A 10 km
0153                                       % observation altitude is assumed here
0154 
0155 
0156 %- Absorption
0157 %
0158 Q.ABS_LINES           = fullfile( atmlab_example_data, 'o3line111ghz' );
0159 Q.ABS_LINES_FORMAT    = 'Arts';
0160 %
0161 Q.ABSORPTION          = 'OnTheFly';
0162 Q.ABS_NLS             = [];
0163 
0164 %- Pressure grid
0165 %
0166 z_toa                 = 95e3;
0167 %
0168 Q.P_GRID              = z2p_simple( Q.Z_SURFACE-1e3 : 2e3 : z_toa )';
0169 
0170 
0171 
0172 %- Frequency, spectrometer and pencil beam antenna
0173 %
0174 % The hypothetical spectrometer has rectangular response functions
0175 %
0176 Q.F_GRID              = qarts_get( fullfile( atmlab_example_data , ...
0177                                                       'f_grid_111ghz.xml' ) );
0178 %
0179 H                     = qartsSensor;
0180 %
0181 H.SENSOR_NORM         = true;
0182 %
0183 H.BACKEND_DO          = true;
0184 df                    = 0.5e6;
0185 H.F_BACKEND           = ( min(Q.F_GRID)+df : df : max(Q.F_GRID)-df )';
0186 %
0187 B.name                = 'Spectrometer channel response function';
0188 B.gridnames           = { 'Frequency' };
0189 B.grids               = { [-df/2 df/2] };
0190 B.dataname            = 'Response';
0191 B.data                = [1 1];
0192 %
0193 H.BACKEND_CHANNEL_RESPONSE{1} = B;
0194 clear B
0195 %
0196 Q.SENSOR_DO           = true;
0197 Q.SENSOR_RESPONSE     = H;
0198 %
0199 Q.ANTENNA_DIM         = 1;
0200 Q.MBLOCK_ZA_GRID      = 0;
0201 
0202 
0203 
0204 %- Correlation of thermal noise
0205 %
0206 f              = H.F_BACKEND;
0207 cl             = 1.4 * ( f(2) - f(1) );
0208 cfun           = 'gau';
0209 cco            = 0.05;
0210 %
0211 Q.TNOISE_C     = covmat1d_from_cfun( f, [], cfun, cl, cco );
0212 %
0213 clear H f
0214 
0215 
0216 %- Define L2 structure (beside retrieval quantities below)
0217 %
0218 Q.L2_EXTRA     = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', ...
0219                    'mresp', 'A', 'e', 'eo', 'es', 'date' };
0220 
0221 
0222 %- Temperature
0223 %
0224 Q.T.RETRIEVE   = false;
0225 Q.T.ATMDATA    = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
0226                         'msis90', 'msis90.t.xml' ), 'Temperature', 't_field' );
0227 
0228 %- Determine altitudes through HSE
0229 %
0230 Q.HSE.ON       = true;
0231 Q.HSE.P        = Q.P_GRID(1);
0232 Q.HSE.ACCURACY = 0.1;
0233 
0234 
0235 %- Species
0236 
0237 % Ozone, only species is retrieved here
0238 Q.ABS_SPECIES(1).TAG      = { 'O3' };
0239 Q.ABS_SPECIES(1).RETRIEVE = true;
0240 Q.ABS_SPECIES(1).L2       = true;
0241 Q.ABS_SPECIES(1).GRIDS    = { Q.P_GRID, [], [] };
0242 Q.ABS_SPECIES(1).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0243       'midlatitude-winter', 'midlatitude-winter.O3.xml' ), 'O3', 'vmr_field' );
0244 %
0245 % If you don't apply a min value (by MINMAX), you could need to active this:
0246 %Q.VMR_NEGATIVE_OK = true;
0247 %
0248 % For demonstration, setting for several units are provided:
0249 switch 1
0250  case 1  % Constant VMR
0251   Q.ABS_SPECIES(1).UNIT     = 'vmr';
0252   Q.ABS_SPECIES(1).SX       = ...
0253                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 1.5e-6, ...
0254                                                'lin', 0.2, 0.00, @log10 ) + ...
0255                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.3e-6, ...
0256                                                'lin', 0.5, 0.00, @log10 );
0257    Q.ABS_SPECIES(1).MINMAX  = 1e-12;
0258  case 2 % Constant rel
0259   Q.ABS_SPECIES(1).UNIT     = 'rel';
0260   Q.ABS_SPECIES(1).SX       = ...
0261                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0262                                                'lin', 0.2, 0.00, @log10 ) + ...
0263                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0264                                                'lin', 0.5, 0.00, @log10 );
0265   Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0266  case 3 % Mimic case 2 in vmr
0267   Q.ABS_SPECIES(1).UNIT     = 'vmr';
0268   Q.ABS_SPECIES(1).SX       = ...
0269                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, ...
0270                                        [ Q.ABS_SPECIES(1).ATMDATA.GRID1,...
0271                                    0.5 * Q.ABS_SPECIES(1).ATMDATA.DATA ],...
0272                                                'lin', 0.2, 0.00, @log10 ) + ...
0273                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, ...
0274                                        [ Q.ABS_SPECIES(1).ATMDATA.GRID1, ...
0275                                    0.1 * Q.ABS_SPECIES(1).ATMDATA.DATA ],...
0276                                                'lin', 0.5, 0.00, @log10 );
0277   Q.ABS_SPECIES(1).MINMAX  = 1e-12;
0278  case 4 % Constant logrel
0279   Q.ABS_SPECIES(1).UNIT     = 'logrel';
0280   Q.ABS_SPECIES(1).SX       = ...
0281                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.5, ...
0282                                                'lin', 0.2, 0.00, @log10 ) + ...
0283                      covmat1d_from_cfun( Q.ABS_SPECIES(1).GRIDS{1}, 0.1, ...
0284                                                'lin', 0.5, 0.00, @log10 );
0285   Q.ABS_SPECIES(1).MINMAX  = 1e-6;
0286 end  
0287   
0288   
0289   
0290 %- Water
0291 %
0292 % This generates no absorption, as linefile has no H2O lines
0293 %
0294 Q.ABS_SPECIES(2).TAG      = { 'H2O' };
0295 Q.ABS_SPECIES(2).RETRIEVE = false;
0296 Q.ABS_SPECIES(2).ATMDATA  = gf_artsxml( fullfile( fascod, ...
0297     'midlatitude-winter', 'midlatitude-winter.H2O.xml' ), 'H2O', 'vmr_field' );
0298 
0299 
0300 
0301 %- Wind
0302 %
0303 % Here demonstrated for v-component, that should the component of main
0304 % concern for ground-based instruments. This component can be seen as the
0305 % part of the horizontal wind going along the viewing direction, where
0306 % positive values mean a movement away from the sensor.
0307 %
0308 Q.WIND_V_FIELD    = [];      % This is a short-cut for zero wind.
0309 %
0310 Q.WIND_V.RETRIEVE = false;   % Set to true to activate retrieval
0311 Q.WIND_V.L2       = true;
0312 Q.WIND_V.GRIDS    = { Q.P_GRID(1:2:end), [], [] };
0313 Q.WIND_V.SX       = covmat1d_from_cfun( Q.WIND_V(1).GRIDS{1}, 40, ...
0314                                                     'lin', 0.5, 0.00, @log10 );
0315 if Q.WIND_V.RETRIEVE
0316   Q.WSMS_AT_START{end+1} = 'IndexSet(abs_f_interp_order,3)';
0317 end
0318 
0319 
0320 
0321 %- Pointing
0322 %
0323 % Here just included for testing purposes, of little interest for ground-based
0324 % spectrometers.
0325 %
0326 Q.POINTING.RETRIEVE       = false;
0327 Q.POINTING.DZA            = 0.1; 
0328 Q.POINTING.POLY_ORDER     = 0;
0329 Q.POINTING.CALCMODE       = 'recalc';
0330 Q.POINTING.SX             = 1;  
0331 Q.POINTING.L2             = true;
0332 
0333 
0334 %- Frequency shift
0335 %
0336 Q.FSHIFTFIT.RETRIEVE      = true;
0337 Q.FSHIFTFIT.DF            = 25e3; 
0338 Q.FSHIFTFIT.SX            = 50e3^2;  
0339 Q.FSHIFTFIT.L2            = true;
0340 
0341 %- Frequency stretch
0342 %
0343 Q.FSTRETCHFIT.RETRIEVE    = false;   % Set to true to activate retrieval
0344 Q.FSTRETCHFIT.DF          = 25e3; 
0345 Q.FSTRETCHFIT.SX          = 50e3^2;  
0346 Q.FSTRETCHFIT.L2          = true;
0347 
0348 
0349 %- Polyfit
0350 %
0351 % A polynomial of order 3 is used for "baseline fit".
0352 %
0353 Q.POLYFIT.RETRIEVE        = true;
0354 Q.POLYFIT.ORDER           = 3;
0355 Q.POLYFIT.L2              = true;
0356 Q.POLYFIT.SX0             = 1^2; 
0357 Q.POLYFIT.SX1             = 0.5^2; 
0358 Q.POLYFIT.SX2             = 0.2^2;
0359 Q.POLYFIT.SX3             = 0.1^2; 
0360 
0361 
0362 %- Sinefit
0363 %
0364 % Here demonstrated with two period lengths
0365 %
0366 Q.SINEFIT.RETRIEVE        = false;   % Set to true to activate retrieval
0367 Q.SINEFIT.PERIODS         = [ 75e3 200e3 ]';
0368 Q.SINEFIT.L2              = true;
0369 Q.SINEFIT.SX1             = 0.2^2; 
0370 Q.SINEFIT.SX2             = 0.4^2; 
0371 
0372 return
0373 %-----------------------------------------------------------------------------
0374 %-----------------------------------------------------------------------------
0375 %-----------------------------------------------------------------------------
0376 
0377 
0378 
0379 
0380 
0381 function Y = y_demo(Q)
0382 
0383 % The data should be loaded from one or several files, but are here genereted
0384 % by a forward model call to show how qpack2 can also be used to generate
0385 % simulaled measurements (matching a priori assumptions).
0386 
0387 % The simulated data model airborn measurements at two different zenith
0388 % angles, from two nearby positions.
0389   
0390 % Init Y
0391 %
0392 Y = qp2_y;
0393 
0394 % Set a date
0395 %
0396 Y.YEAR  = 2008;
0397 Y.MONTH = 2;
0398 Y.DAY   = 25;
0399 
0400 % Lat / lon
0401 %
0402 Y.LATITUDE  = 45;
0403 Y.LONGITUDE = exp(1);
0404 
0405 % An airborn measurement assumed here
0406 %
0407 Y.Z_PLATFORM = 10.5e3;
0408 Y.ZA         = 50;
0409 
0410 % Reference point for hydrostatic equilibrium
0411 %
0412 Y.HSE_P = 100e2;
0413 Y.HSE_Z = 16e3;
0414 
0415 % Set backend frequencies
0416 %
0417 Y.F = Q.SENSOR_RESPONSE.F_BACKEND;
0418 
0419 % Thermal noise standard deviation
0420 %
0421 Y.TNOISE = 0.1;
0422 % To test varying noise
0423 %Y.TNOISE = linspace( 0.03, 0.07, length(Y.F) )';
0424 
0425 % Simulate a measurement
0426 %
0427 Y.Y = [];                           % A flag to tell qpack2 to calculate the
0428 %                                     spectrum ({} signifies undefined!).
0429 
0430 % Add a second measurement
0431 %
0432 Y(2) = Y(1);
0433 %
0434 Y(2).LONGITUDE = pi;
0435 Y(2).ZA        = 45;
0436 
0437 
0438 % Calculate simulated spectra
0439 %
0440 Y = qpack2( Q, oem, Y );            % Dummy oem structure OK here
0441 
0442 % Add thermal noise
0443 %
0444 % The correlation specified in Q is included
0445 %
0446 for i = 1 : length(Y)
0447   Y(i).Y = Y(i).Y + Y(i).TNOISE .* make_noise(1,Q.TNOISE_C);
0448 end
0449 
0450 % Add a constant "baseline shift" for measurement 2
0451 %
0452 Y(2).Y = Y(2).Y + 1;
0453

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