Chapter 4. How to use GRASP

4.1. How to run the code

After the installation of GRASP software the grasp command will be available. As it is explained in a previous section (see Section 3.6, “Running the code”) just typing grasp will print some general information about how the software was compiled.

First arguments that grasp executable expects is a path to the settings file. The settings file describes the inversion strategy and general behaviour of the process: where is the input data, in which format is the input data, where to store the output results... Therefore, a deep knowledge of settings parameters is the base to understand GRASP.

4.1.1. Settings file

Settings file is written in YAML format and that brings many benefits: easy to write, clear to read, self-explanatory names, flexible and powerful. The concepts are organized in blocks that are translated to YAML thanks to the fix indentation (we suggest 4 white-spaces). So, for example, first level defines the different modules:

                    

input:
    # Here settings related with input module
    segment: # It define a description of input segment
        x: 2 # It define the maximum size of x dimension of the segment to 2.
output:
    # Settings linked with the output
    # ...
retrieval:
    # Definition of inversion strategy
    # ...
...
            

When grasp is called, the first step is to read settings and the second is to prepare the environment for the settings defined in the main structures. If the settings are not valid, a pretty clean error message will be printed. Please, read the first line of it carefully to understand the error.

The settings parameters can also be defined by the command line. In this case, after the first argument (settings file name) extra settings parameters can be defined with the syntax key=value where key is the parameter name in "dot syntax". For example, for GRASP is equivalent to be called with the argument input.driver=sdata or to have defined in the settings file the following content:

input:
    driver: sdata
                

It is important to define clearly how the relative paths have to be defined in the settings files. Relative paths are always relative to the file that defines it. In the next section the reader will learn about how to include other settings file inside of a settings file, but this rule will keep valid: Relative paths are defined from settings file that define it. In the case of usage of command line, relative paths are relative to the current working directory. In the case of absolute paths all this complexity disappears but the results are less portables.

4.1.1.1. HELP argument

The list of available parameters for GRASP is long. The "help" argument will help you to know the available parameters or to look for something specific. When help argument is present GRASP is not executed normally, but instead, the help information will be printed into screen. Additionally, help can be followed by a search string to filter the results. For example help=input will print only the settings which contain "input" string in its definition and help=input.segment will print only the settings which will be under the block "segment" inside of the block "input" (because "dot" symbol defines block separator).

Table 4.1. The SDATA main structure

Field NameField Content
help Shows this help information
version Shows GRASP code version information and stops
global. resources Path to framework resources folder
retrieval. general. path_to_internal_files Path to internal data files
retrieval. convergence. stop_before_performing_retrieval Define how the code is run (true=only forward model, false=full inversion)
retrieval. convergence. minimization_convention Minimization in absolute or logarithm space
retrieval. convergence. threshold_for_stopping_Q_iterations Threshold for stopping q - iterations: once the change in the residual smaller than this parameter the iterations are stopped
retrieval. convergence. scale_for_finite_difference Defines DL for calculations of derivatives: (f(x+DL)-f(x))/DL
retrieval. convergence. threshold_for_stopping Threshold for stopping p - iterations: once the change in the residual smaller than this parameter the iterations stops
retrieval. convergence. normal_system_solver Defines the method to solve normal system
retrieval. convergence. maximum_iterations_of_Levenberg-Marquardt The maximum number of ip iterations where Levenberg-Marquardt correction is applied
retrieval. convergence. maximum_iterations_for_stopping Maximum number of iterations performed in the retrieval before stop
retrieval. convergence. shift_for_applying_logarithm_to_negative_values This value (usually 1) is added to negative parameters in order to be able to apply logarithmic transformations
retrieval. regime_of_multipixel_constraints. inversion_regime Flag for single multi-pixel retrieval (VALUES)
retrieval. regime_of_measurement_fitting. polarization 1. -We fit: I, Q, U; 2. - We fit: I, Q/I, U/I; 3. - We fit: I, (Q^2+U^2)^(1/2) (Polarized measurements can be defined in input data as Q and U or P); 4- We fit: I, [(Q^2+U^2)^(1/2)]/I (Polarized measurements can be defined in input data as Q and U or P); 5- We fit: P/I (User has to provide P/I in input data)
retrieval. regime_of_measurement_fitting. radiance 1. -We fit: I; 2. - We fit I/sum(I(:))
retrieval. product_configuration. wavelength_indices_for_angstrom Indices of wavelengths which will be used to calculate Angstrom exponent
retrieval. product_configuration. aerosol_particulate_matter_diameter Diameters of aerosol particles in microns which will be used to calculate Particulate Matter (PM)
retrieval. product_configuration. wavelength_indices_for_ndvi Indices of wavelengths which will be used to calculate NDVI if it is calculated
retrieval. product_configuration. wavelenght_indices_for_aod_error_estimation Indices of wavelengths which will be used to estimate aod error
retrieval. product_configuration. wavelenght_indices_for_ssa_error_estimation Indices of wavelengths which will be used to estimate ssa error
retrieval. product_configuration. wavelenght_indices_for_lidar_error_estimation Indices of wavelengths which will be used to estimate lidar error
retrieval. phase_matrix. size_binning_method_for_triangle_bins Determining the scale used for the binning of size distribution
retrieval. phase_matrix. number_of_elements Number of phase matrix elements used in the calculations and retrieval
retrieval. phase_matrix. use_models Use GRASP 'models' approach which produce results less acurated (based on models) but the process is much quicker
retrieval. phase_matrix. kernels_folder Path to kernels when we retrieve size distribution for triangle bins or lognormal size distribution
retrieval. phase_matrix. radius. mode[]. bins Number of triangle bins set for size distribution representation
retrieval. phase_matrix. radius. mode[]. min Minimum value of the radius used in the triangle bins
retrieval. phase_matrix. radius. mode[]. max Maximum value of the radius used in the triangle bins
retrieval. edges_size. x Size of edges width in pixels (it have to be lower than KIEDGE compilation constant)
retrieval. edges_size. y Size of edges height in pixels (it have to be lower than KIEDGE compilation constant)
retrieval. edges_size. t Size of temporal dimension of edges (it have to be lower than KIEDGE compilation constant)
retrieval. radiative_transfer. number_of_layers Maximum number of vertical layers resolved in radiative transfer (minimum and maximum value. If only one number is assigned maximum is by default KNT-1)
retrieval. radiative_transfer. molecular_profile_vertical_type It defines which model will be used to describe vertical profile of molecular (Rayleigh) scattering. Values: 'Exponential' describes molecular density profile at altitude h as exp(h/8)/8, 'Standard_atmosphere' uses standard atmosphere model (pressure and temperature profiles) to calculate molecular density at each altitude.
retrieval. radiative_transfer. aerosol_profile_vertical_type It defines which model will be used to describe vertical profile of aerosol distribution. Values: 'Exponential' describes aerosol concentration profile at altitude h as exp(h/HM)/HM, 'gaussian' uses normal distribution exp(h-HM)/(sqrt(rpi)*sqrt(2. )*sigma) with sigma and HM parameter taken from parameters.
retrieval. radiative_transfer. phase_matrix_truncation Switch on/off the truncation: the technique to calculate scattering effects from the sharp forward peak of the phase function separately from those of the rest of the phase function. It doesn't effect considerably the accuracy while provides much faster calculations
retrieval. radiative_transfer. absolute_error_rt_calculations Absolute value of truncation threshold of Fourier and order-of-scattering series expansions in radiative transfer calculations
retrieval. radiative_transfer. reference_plane_for_polarization Reference plane for polarization calculations
retrieval. radiative_transfer. look_up_table. mode LUT Mode for radiative transfer (disable by default)
retrieval. radiative_transfer. simulating_observation. order_of_scattering Regime of scattering used for modeling diffuse radiation observations
retrieval. radiative_transfer. simulating_observation. number_of_gaussian_quadratures_for_expansion_coefficients Number of Gaussian quadratures for calculating expansion coefficients used for multiple scattering simulations
retrieval. radiative_transfer. simulating_observation. number_of_guassian_quadratures_for_fourier_expansion_coefficients Number of Gaussian quadratures for calculating Fourier expansion coefficients used for multiple scattering simulations
retrieval. radiative_transfer. simulating_observation. number_of_fourier_expansion_coefficients Number of Fourier expansion coefficients used in multiple scattering simulations
retrieval. radiative_transfer. simulating_derivatives. order_of_scattering Regime of scattering used in calculation of radiance derivatives
retrieval. radiative_transfer. simulating_derivatives. number_of_gaussian_quadratures_for_expansion_coefficients Number of Gaussian quadratures for calculating expansion coefficients used in calculation of radiance derivatives
retrieval. radiative_transfer. simulating_derivatives. number_of_guassian_quadratures_for_fourier_expansion_coefficients Number of Gaussian quadratures used for calculating Fourier expansion coefficients used in calculation of radiance derivatives
retrieval. radiative_transfer. simulating_derivatives. number_of_fourier_expansion_coefficients Number of Fourier expansion coefficients used in calculation of radiance derivatives
retrieval. noises. noise[]. standard_deviation_synthetic Standard deviation of synthetic random noise added to the corresponding inverted data, if value is 0 no synthetic random noise is added
retrieval. noises. noise[]. error_type Error type used for definition of covariance matrices used in measurement fitting and in modeling synthetic noise
retrieval. noises. noise[]. standard_deviation Standard deviation of the random noise expected (this defines the covariance matrix used in fitting)
retrieval. noises. noise[]. measurement_type[]. type Type of relevant measurements for applying the noise assumptions
retrieval. noises. noise[]. measurement_type[]. index_of_wavelength_involved List of indices of relevant wavelengths for applying the noise assumptions
retrieval. products. aerosol. chemistry Retrieve aerosol chemical composition (if retrieved)
retrieval. products. aerosol. lidar Retrieve columnar lidar ratios (e. g. , if lidar data are inverted)
retrieval. products. aerosol. optical_properties Provide aerosol optical properties
retrieval. products. aerosol. phase_matrix Obtain aerosol phase matrix
retrieval. products. aerosol. refractive_index Provide aerosol refractive index
retrieval. products. aerosol. theoretical_bimodal_extinction Provide estimated aerosol extinction for fine and coarse modes
retrieval. products. aerosol. theoretical_bimodal_parameters Provide estimated aerosol microphysical characteristics for fine and coarse modes
retrieval. products. aerosol. particulate_matter Obtain aerosol particulate matter estimation at given particle diameters
retrieval. products. aerosol. type Obtain aerosol type index i. e. 0-Complex mixture,1-Background,2-Maritime,3-Urbn. Polluted,4-Mixed,5-Urbn. clean,6-Smoke flam. ,7-Smoke. sold. ,8-Dust
retrieval. products. error_estimation. aerosol. lidar Implement error estimation for aerosol lidar products
retrieval. products. error_estimation. aerosol. optical_properties Implement error estimation for optical properties of aerosol products
retrieval. products. error_estimation. parameters Implement error estimation for retrieved parameters
retrieval. products. forcing. broadband_flux Provide upward and downward fluxes integrated over the solar spectrum at the user-defined levels for the scenarios with and without aerosols
retrieval. products. forcing. forcing Provide aerosol radiative forcing values at the user-defined levels
retrieval. products. retrieval. fitting Provide obtained measurement fitting
retrieval. products. retrieval. parameters Provide retrieved parameters
retrieval. products. retrieval. residual Provide values of obtained residuals
retrieval. products. surface Provide surface reflectance products
retrieval. debug. verbose Retrieval prints progress information while it is performing the inversion
retrieval. debug. additional_information Print some additional information
retrieval. debug. simulated_sdata_file Filename where simulated observation data are to be printed. This option force retrieval. general. stop_before_performing_retrieval=true
retrieval. debug. path_to_extra_files Path to folder with retrieval extra files: debug information, extra resources like image. dat files. . .
retrieval. debug. use_internal_initial_guess Test option which allows to load different initial guess for each pixel (loading them from image. dat files)
retrieval. constraints. characteristic[]. type Type of characteristic
retrieval. constraints. characteristic[]. retrieved Specify if this characteristic will be retrieved or only used in forward model
retrieval. constraints. characteristic[]. mode[]. initial_guess. value Initial values for a specific (determined by type) characteristic
retrieval. constraints. characteristic[]. mode[]. initial_guess. min Minimum value for the specific characteristic
retrieval. constraints. characteristic[]. mode[]. initial_guess. max Maximum value for the specific characteristic
retrieval. constraints. characteristic[]. mode[]. initial_guess. index_of_wavelength_involved Indices of Wavelengths associated to the specific characteristic
retrieval. constraints. characteristic[]. mode[]. initial_guess. estimate_error Flag to retrieve the error of specific parameter if retrieval. products. error_estimation. parameters is true
retrieval. constraints. characteristic[]. mode[]. single_pixel. a_priori_estimates. lagrange_multiplier Value of the Lagrange multiplier associated to a priori estimate of the retrieved characteristics (applied in each single pixel)
retrieval. constraints. characteristic[]. mode[]. single_pixel. smoothness_constraints. difference_order Order of the derivatives/differences used for applying a priori smoothness constrains for the retrieved characteristics
retrieval. constraints. characteristic[]. mode[]. single_pixel. smoothness_constraints. lagrange_multiplier Value of the Lagrange multiplayer used for applying a priori smoothness constraints for the retrieved characteristics
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. derivative_order_of_X_variability Order of derivatives/differences used for applying a priori smoothness constraints on parameter inter-pixel X-variability in multi-pixel retrieval regime
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. lagrange_multiplier_of_X_variability Value of the Lagrange multiplier used for applying a priori smoothness constraints on parameter inter-pixel X-variability in multi-pixel retrieval regime
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. derivative_order_of_Y_variability Order of derivatives/differences used for applying a priori smoothness constraints on parameter inter-pixel Y-variability in multi-pixel retrieval regime
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. lagrange_multiplier_of_Y_variability Value of the Lagrange multiplier used for applying a priori smoothness constraints on parameter inter-pixel Y-variability in multi-pixel retrieval regime
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. derivative_order_of_T_variability Order of derivatives/differences used for applying a priori smoothness constraints on parameter inter-pixel T-variability in multi-pixel retrieval regime
retrieval. constraints. characteristic[]. mode[]. multi_pixel. smoothness_constraints. lagrange_multiplier_of_T_variability Value of the Lagrange multiplier used for applying a priori smoothness constraints on parameter inter-pixel T-variability in multi-pixel retrieval regime
retrieval. untested. ratio[]. value Values of bins using in retrieved aspect or axis distributions
settings. debug Shows settings read to run this program
settings. strict Force GRASP to continue when there were parse or validation errors in the settings file
settings. dump Stream to dump read settings in short format(experimental)
settings. long_dump Stream to dump read settings in long format (experimental)
input. driver The driver that will be called for inverting data
input. file Name of file(s) which contain input observation
input. center. latitude Latitude of the center of tile to invert
input. center. longitude Latitude of the center of tile to invert
input. corner. row The number of the row for input driver in the native coordinate system of the sensor
input. corner. column The number of the column for input driver in the native coordinate system of the sensor
input. grid_offset. row Information of the first row in the input grid that will be use for normalizing the output. If 0 is used the output coordinate reference will be the same than the input. This offset is for forcing the output row to start at 0 (e. g. if you put 1, output_row == input_row - 1)
input. grid_offset. column Information of the first column in the input grid that will be used for normalizing the output. If 0 is used the output coordinate reference will be the same than the input. This offset is for forcing the output column to start at 0 (e. g. if you put 1, output_column == input_column - 1)
input. area. width The width of the covered area in pixels. It has to be divisible by intput. segment. x value
input. area. height The height of the covered area in pixels. It has to be divisible by intput. segment. y value
input. time. from Initial date and time for data processing
input. time. to Final date and time for data processing
input. segment. x Size of segment width in pixels (it have to be lower than KIX compilation constant)
input. segment. y Size of segment height in pixels (it have to be lower than KIY compilation constant)
input. segment. t Size of segment temporal dimension (it have to be lower than KITIME compilation constant)
input. transformer Name of input data transformer functions to be used after load data
input. debug. raw_segment Stream to print raw segment data loaded
input. debug. clean_segment Stream to print segment information after clean NaN values
input. debug. used_files Stream to print names of the files that have the pixels for inverting
input. sdata. dump Stream where to dump sdata information
input. sdata. dump_original Stream where to dump sdata information just after being generated by the driver. Some transformers can modify sdata information, this setting is thought for debugging purposes where the user is interested in knowing sdata generated by the driver instead of data driving inside the inversion.
input. imagedat. dump Stream where to dump initial guess information (image. dat format)
input. preload_segment. x This parameter specifies how many segments in X dimension will be preloaded in each block
input. preload_segment. y This parameter specifies how many segments in Y dimension will be preloaded in each block
input. preload_segment. t This parameter specifies how many segments in T dimension will be preloaded in each block
input. driver_settings. sdata. debug Print debug information from sdata reader subsystem
input. transformer_settings. segment_imagedat. file File which contains initial guess information in classic input. dat format
output. segment. function Driver to process output for every single retrieval (show information in screen, perform a map, plotting, . . . )
output. segment. stream Stream to dump segment output data
output. iteration. function Driver to process output for every single retrieval (show information in screen, perform a map, plotting, . . . )
output. iteration. stream Stream to dump segment output data
output. tile. function Driver to process output after processing complete tile (show information in screen, perform a map, plotting, . . . )
output. tile. stream Stream to dump tile output data
output. current. function Driver to process output after each retrieval (show information in screen, perform a map, plotting, . . . )
output. current. stream Stream to dump current progress information about the retrieval conducted
output. segment. function_settings. csv. delimiter Separator between fields
output. segment. function_settings. csv. compression If true output is compressed in GZ format (gz extension is automatically added)
output. tile. function_settings. csv. chemical_concentration Calculate and print chemical concentration
output. tile. function_settings. csv. delimiter Separator between fields
output. tile. function_settings. csv. compression If true output is compressed in GZ format (gz extension is automatically added)
controller. segment_range This parameter allows specifying a range of segments that will be inverted/processed. If it is a single number a specific retrieval will be processed. Use -1 as undefined. For example [15, -1] will process all retrievals starting by the segment #15
controller. debug. perform_retrieval Allowing controller to call retrieval. If this parameter is false the framework will work without inverting the data and it will force to use only none output functions (results will not be printed). This parameter is useful for debug the framework or prepare input data
controller. debug. compilation_information Print compilation information at the beginning of the process
controller. debug. tracking_memory_stream Stream where printing all information about memory allocated during the execution (debugging information)
controller. stream Stream to dump controller information
controller. mpi. maximum_job_time Maximum number of seconds that the master node will wait for a working node for obtaining results. After this time (specified in seconds) if the job is not finished the controller will kill the task and the segment will be skipped
controller. mpi. polling_time Number of seconds that the master node wait after checking the workers

4.1.1.2. Extending settings: command line, import and template statements

The first argument of GRASP has to be the settings file but this file can be modified by another mechanism proposed by the settings module. The main way is by the command line, which allows to replace every settings parameter with "dot" syntax (replacing indentation by "dot" symbol and colon symbol by equal). All parameters that have been defined before and being replaced in the command line will cause a "note" information during the execution of GRASP. That sentence is just to inform the user that command line arguments always have higher priority than parameters in settings files. The value from command line will be the value that will be used to run the code. The command line is a very powerful feature to be used in production scripts.

But the command line is not the only way to modify GRASP settings files. Settings files accept "import" and "template" statement. This statements could look similar but theirs behaviour is a bit different. Both of them allow defining other settings files that are read before the current one, but in case of import the settings can not be overwritten. Template statement allows loading other settings files and then, modifying some settings to customize the loaded file. Take into account that these statements can be used in cascade, creating problems to debug the code. Please, use these statements carefully.

4.1.1.3. Streams

Some of the settings parameters are defined as "streams". GRASP output streams allow users to create dynamic names avoiding overwriting files or having to change the filenames each time they execute GRASP. When a description of a parameter is defined as "output stream" the user can set up as a regular output path, for example ./folder/file.extension or use the "magic" behind the output streams by using a wildcard that will be replaced by dynamic values. For instance:

                    

output:
    segment:
        function: hdf
        stream: "GRASP_Banizoumbou_20080101_20080331_2x2+3286+1376.hdf"
    tile:
        function: [ ascii, hdf ] 
        stream: [ "GRASP_Banizoumbou_20080101_20080331.txt",
                "GRASP_Banizoumbou_20080101_20080331.hdf" ]
                

That definition is ok for many cases but if many tiles or segments are going to be processed the fixed names will produce name collisions (the content of some files will be overwritten during the process). It would be tedious (and not always possible) for the user to change oneself the dates or other numeric substrings in the file names. For this reason, the configuration system provides some wildcards that will be automatically replaced with the given values depending on the state of the processing. These wildcards are marked with curly braces and their names are quite self-explanatory. The previous example can be rewritten in a more generic way using the stream wildcards:

                    

output:
    segment:
        function: hdf
        stream: "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}_{segment_nx}x{segment_ny}+{segment_corner_column(4)}+{segment_corner_row(4)}.hdf"
    tile:
        function: [ ascii, hdf ] 
        stream: [ "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}.txt",
                "GRASP_Banizoumbou_{tile_from(%Y%m%d)}_{tile_to(%Y%m%d)}.hdf" ]
                

In addition wildcards will provide the user the capability to set some system streams. If you use the values "true", "screen", "t" or "1" the information will be printed in the terminal (stdout). If the stream is set to "false", "none", "null", "f" or 0 nothing will be printed (like redirect to /dev/null). Finally, using the value "stderr" or "-1" the output will be redirected to the standard error output.

The following list show all available wildcards that can be used for creating dynamic output filenames:

  • auto(N): itime x icol x irow with N zeros at the left
  • icol(N): current column number with N zeros at the left
  • irow(N): current column number with N zeros at the left
  • itime(N): current column number with N zeros at the left
  • iinversion(N): current inversion id with N zeros at the left
  • segment_nx(N): number of X elements per segment with N zeros at the left
  • segment_ny(N): number of Y elements per segment with N zeros at the left
  • segment_nt(N): number of T elements per segment with N zeros at the left
  • tile_from(FORMAT): start tile date in FORMAT. By default FORMAT is %FT%H_%M_%SZ
  • tile_to(FORMAT): final tile date in FORMAT. By default FORMAT is %FT%H_%M_%SZ
  • tile_corner_column(N): number of the corner (column) of the segment defined in settings file. Requirement: Input data have to be defined using input.corner instead of input.center
  • tile_corner_row(N): number of the corner of (row) the segment defined in settings file. Requirement: Input data have to be defined using input.corner instead of input.center
  • tile_center_longitude(FORMAT): longitude of the center of the tile defined in settings file. Requirement: Input data have to be defined using input.center instead of input.corner
  • tile_center_latitude(FORMAT): latitude of the center of the tile defined in settings file. Requirement: Input data have to be defined using input.center instead of input.corner
  • tile_coordinate_x(I): x input reference of center of the tile defined in settings file. It can be defined by corner or latitude. If is N in case it was defined by corner or 0.I in case it was defined like center
  • tile_coordinate_y(I): y input reference of center of the tile defined in settings file. It can be defined by corner or latitude. If is N in case it was defined by corner or 0.I in case it was defined like center
  • tile_width(N): Number of X elements in tile with N zeros at the left
  • tile_height(N): Number of Y elements in tile with N zeros at the left
  • segment_corner_column(N): number of column of the segment corner with N zeros at the left. Requirement: Input data have to be defined using input.corner instead of input.center
  • segment_corner_row(N): number of row of the segment corner with N zeros at the left Requirement: Input data have to be defined using input.corner instead of input.center
  • segment_first_date(FORMAT): date of first pixel inside the segment in FORMAT. By default FORMAT is %FT%H_%M_%SZ
  • segment_last_date(FORMAT): date of last pixel inside the segment in FORMAT. By default FORMAT is %FT%H_%M_%SZ
  • iteration(N): Number of iterations with N zeros at the left. (Note. In case of single pixel it returns the number of iterations of first pixel)
  • settings_filename: the name of settings file used to run the retrieval
  • version: version of grasp if it is compiled with saving this information
  • branch: git branch of grasp if it is compiled with saving this information
  • commit: reference of git commit of grasp if it is compiled with saving this information
  • constants_set: constants set used in compilation time
  • pwd: this is replaced by current folder and is only valid at the beginning of the stream definition
  • yml: this is replaced by current folder of main configuration file and it is only valid at the beginning of the stream definition

4.1.2. Retrieved characteristics

Understanding GRASP inversion procedure means understanding how the algorithm is starting from an initial guess and obtains a results array. This is an iterative procedure explained in the literature and introduced in Section 2.3, “GRASP scientific core algorithm”. The purpose of this section is to explain how initial guess is represented inside the code as an array which evolve in each iteration until get the result array. Following diagram shows how initial guess is read from settings file and translated to an internal array in the code. This detail could look very technical and related with the development but understanding of some internal concepts of the code helps to understand how it works. Following diagram shows this transformation:

Figure 4.1. Translation of settings file into initial guess array

Translation of settings file into initial guess array


Once that the initial guess is loaded, it is set as first array of characteristics to be retrieved. Retrieval process iterates over it until it gets the results. The results of GRASP retrieval is the array with the same shape as the initial guess but containing the results of retrieval of these parameters to match them with SDATA file. Then, the GRASP forward model is called using this results array, GRASP obtains the rest of the derived results it provides. So, we can talk about two kinds of results: basic results and derived products. Following diagram shows this process:

Figure 4.2. Evolution of retrieved characteristics during GRASP processing

Evolution of retrieved characteristics during GRASP processing


To know all details about the products obtained by GRASP please, see Section 4.3.1, “The list of GRASP output parameters”.

4.2. Input module

The input module is responsible for reading the data and for setting the internal input data structures. GRASP provides a very flexible way to inject the input data offering developers the capability to create the input drivers. An input driver is an extension of GRASP that is added during compilation and is selected in the settings file. These generic drivers allow the developers to create a custom way to read a specific database and load it into the GRASP scientific algorithm. Therefore, GRASP can read infinite kinds of input databases, as much as drivers exist. Additionally, drivers can take the responsibility of performing some pre-processing actions (such as calibration corrections) before retrieving the data. Finally, when GRASP is used for massive data processing, drivers are extremely important since they help to prepare input data and load it into scientific module without the use of any intermediate files, everything is performed in memory. For general use GRASP proposes a generic driver called the SData driver, which reads files in SData format (Sensor-Data format). These files are not standard, they have a specific format proposed by GRASP to start working with the code. This format is described in Section 4.2.1, “The SDATA format”.

Transformers are the second kind of extensions that the input module offers to users (and developers). They are called after getting input data for a segment, it means after calling the driver. The purpose of transformers is to add the capability to modify the segment after obtaining the data. An example of transformation is to load a climatology database and to modify the initial guess of each pixel to optimize the number of iterations needed to retrieve the data. Performing this action after getting the data allows reusing it between different drivers. Default installation of GRASP does not offer any transformers. All of them are considered as optional extensions and have to be externally added manually or using the grasp-manager.

4.2.1. The SDATA format

The SDATA (sensor data) format is the original input data format of the GRASP code. It is a simple text format designed by the science team at the early stages of the development of the scientific code and it's the easiest way to create test data.

In the context of the GRASP project, the SDATA format has a number of pros:

  • Designed by the scientific team at the origin of the GRASP project, it is well adapted to its needs.
  • It is very simple to describe.
  • It is a text format, and therefore portable, quite easy to check (for the accustomed eye!) and to edit and to make some quick experiments.
  • It is a piece of cake to read in Fortran :-)

It has also a number of cons:

  • It is not standard (a by-product of the "not designed by a committee" approach). No off-the-shelf library is available to parse it and validate it (meaning outside of the GRASP project).
  • It lacks of flexibility: the order and the number of values are fixed for one version and any change in the format to meet new requirements is likely to break the compatibility with the former versions.
  • It is fragile: a malformed file may easily make the code crash or produce mysterious bugs. The design of the format, while simple, makes it hard to develop a really reliable validator.
  • Comments are expected only after values, on the same line (after a colon sign)
  • While the format uses a text representation for the data, it contains lots of numeric values, with limited accuracy and with no comment. Large files are tedious to read and it's easy to make shift errors while reading and editing even for the experienced user.
  • Being a text format, it becomes very inefficient for large data volumes. While it can be compressed for archiving, it must be uncompressed for processing, and only sequential access is possible. It is still possible to perform regional processings (several dozens of thousands of pixels that cover more than a few hundreds of kilometres in both directions) with this format (it was actually done for the sake of necessity), but it stresses the computing system a lot and can't be scaled up to the global processing.

Whatever the number and seriousness of the cons, one of the design objectives is to keep the code simple and flexible allowing scientific community to play with the code. It does not make sense to implement a driver for a single user who wants to do some tests with GRASP. That's why this easy format is maintained by the developer team.

In the following description elements in fixed-width font are the snippets of content. The numeric values in these snippets (e.g. in 2 2 2 : NX NY NT) are given only as examples.

Figure 4.3. An example of SDATA file

        

SDATA version 2.0
2   2   2  : NX NY NT

4       2008-01-04T13:15:00Z    70000.0   0   0   : NPIXELS  TIMESTAMP  HOBS_km  NSURF  IFGAS
1	    1		1	 3286	     1377   2.599	13.528       252.0	 100.0      6  0.443   0.490   0.565 ...
2	    1		1	 3287	     1377   2.657	13.528       242.0	 100.0      6  0.443   0.490   0.565 ...
1	    2		1	 3286	     1376   2.601	13.583       241.0	 100.0      6  0.443   0.490   0.565 ...
2	    2		1	 3287	     1376   2.658	13.583       239.0	 100.0      6  0.443   0.490   0.565 ...


4       2008-01-06T13:02:41Z    70000.0   0   0   : NPIXELS  TIMESTAMP  HOBS_km  NSURF  IFGAS
1	    1		1	 3286	     1377   2.599	13.528       252.0	 100.0      6  0.443   0.490   0.565 ...
2	    1		1	 3287	     1377   2.657	13.528       242.0	 100.0      6  0.443   0.490   0.565 ...
1	    2		1	 3286	     1376   2.601	13.583       241.0	 100.0      6  0.443   0.490   0.565 ...
2	    2		1	 3287	     1376   2.658	13.583       239.0	 100.0      6  0.443   0.490   0.565 ...
                

A SDATA file has a simple structure:

  1. The first line is the FILE HEADER. It contains a magic identifier SDATA followed by a version number.
  2. The second line is the SEGMENT HEADER. It contains three numbers, NX, NY and NT, the spatial and temporal dimensions of the segment that this SDATA file represents. You can notice a colon and names of fields after the values. This is the way how comments are written in the SDATA files. Everything starting from the colon will be ignored by the SDATA parser.
  3. An empty line follows the SEGMENT HEADER.
  4. Then comes the first CELL (group of neighbouring pixels) of the SEGMENT. Each CELL has at most NX*NY pixels (but it may have less, for various reasons: cloudy pixels that have been filtered, missing pixels, etc.). The number of CELLs in the SEGMENT is given by the NT number provided in the SEGMENT HEADER.

Table 4.2. The SDATA main structure

Field NameField Content
FILE HEADER SDATA version 1.0
SEGMENT HEADER

2 2 2 : NX NY NT

empty line 
CELL 1cell content, look for CELL structure
empty line 
CELL 2cell content, look for CELL structure
empty line 
... 
CELL itcell content, look for CELL structure
empty line 
... 
CELL NTcell content, look for CELL structure
empty line 


A CELL is a set of neighbouring pixels, that form the base of a SEGMENT. Each CELL has a HEADER and a number of PIXELs, supposedly acquired at the same time.

The CELL HEADER contains:

  1. The number of pixels in the CELL (NPIXELS). It may not be larger than NX*NY
  2. The timestamp of acquisition of the pixels, in the ISO8601 time format.
  3. A "height" of observation, in metres. The value here is a bit weird (70000), and doesn't correspond to the satellite altitude (that is at least 10 times larger). Actually the value doesn't really matter as long as it is large. Historically, the scientific team has used this value of 70000 in many SDATA files.
  4. Two values for the number of surface and gas parameters. These two values are currently not documented and can be set to 0 for the moment.
  5. Comments starting with a colon.

Table 4.3. The CELL structure

Field NameField Content
CELL HEADER
4 2008-01-04T13:15:00Z 70000.0 0 0 : NPIXELS TIMESTAMP HOBS NSURF IFGAS
PIXEL 1a line of values, look for PIXEL structure
PIXEL 2a line of values, look for PIXEL structure
... 
PIXEL NPIXELSa line of values, look for PIXEL structure


Each line of data after the CELL HEADER represents exactly one pixel, with all its fields. The Table 4.4, “The PIXEL structure” describes the order and type of these fields. For the types, the Fortran notation is used: array types are described with the dimensions of arrays between parentheses, and the ordering is such that the first index increases faster. Indices start from 1, not from 0 like in C. For instance, when one reads real(nwl) for wavelengths, that means that one has to read a list of nwl real values that represent wavelengths.

Table 4.4. The PIXEL structure

Field TypeVariable Name (in source code)Field Content
integerpixel[ipix].ixcoordinate x in the current cell, starting at 1 (in the direction EW)
integerpixel[ipix].iycoordinate y in the current cell, starting at 1 (in the direction NS)
integerpixel[ipix].icloudycloud flag: 0 = cloud, 1 = clear [a]
integerpixel[ipix].icolline of the pixel in its original grid or database (can be set to 0 when not relevant)[b]
integerpixel[ipix].rowcolumn of the pixel in its original grid or database (can be set to 0 when not relevant)[b]
realpixel[ipix].x

longitude of the pixel, in decimal degrees, in the range [-180..180]

0: Greenwich meridian, east of Greenwich: positive, west of Greenwich: negative

realpixel[ipix].y

latitude of the pixel, in decimal degrees, in the range [-90..90]

realpixel[ipix].MASLaltitude of the ground, in metres (MASL: metres above sea level)
realpixel[ipix].land_percentpercentage of land, in the range [0 (sea) .. 100 (land)]. Intermediate values correspond to coastal pixels
integerpixel[ipix].nwlnumber of available wavelengths (nwl)
real(nwl)pixel[ipix].meas[nwl].wllist of wavelengths, in micrometers
integer(nwl)pixel[ipix].meas[nwl].nipnumber of types of measurements for each wavelength (nip)
integer(nip, nwl)pixel[ipix].meas[nwl].meas_type[nip]

list of types of measurements meas_type (see Table 4.5, “TYPES OF MEASUREMENTS”)

The ordering is as follows: meas_type(1, wl1) meas_type(2, wl1) ... meas_type(nip, wl1) meas_type(1, wl2) ... meas_type(nip, wln)

integer(nip, nwl)pixel[ipix].meas[nwl].nbvm[nip]

number of valid measurements (nbvm), for each type of measurement and for each wavelength

The ordering is as follows: nbvm(1, wl1) nbvm(2, wl1) ... nbvm(nip, wl1) nbvm(1, wl2) ... meas_type(nip, wln)

real(nwl)pixel[ipix].meas[nwl].szasolar/sounding zenith angle (sza or θs) in decimal degrees ([0..90]), for each wavelength
real(nbvm, nip, nwl)pixel[ipix].meas[nwl].thetav[nip][nbvm]

viewing zenith angle (vza or θv) in decimal degrees ([0..90] for PARASOL, TBD for PHOTOMETERS)

The ordering is as follows: θv(1, 1, wl1) vza(2, 1, wl1) ... θv(nbvm, 1, wl1) θv(1, 2, wl1) ... θv(nbvm, nip, wln)

In case of lidar or vertical observations this field contains altitudes or ranges of the observations in meters.

real(nbvm, nip, nwl)pixel[ipix].meas[nwl].phi[nip][nbvm]

relative azimuth angle (raa or Δϕ) in decimal degrees ([-180..180] or [0..360])

The ordering is as follows: Δϕ(1, 1, wl1) Δϕ(2, 1, wl1) ... Δϕ(nbvm, 1, wl1) Δϕ(1, 2, wl1) ... Δϕ(nbvm, nip, wln)

real(nbvm, nip, nwl)

pixel[ipix].meas[nwl]

.tau[nbvm];...;.P[nbvm]

measurements (depending on meas_type), for each wavelength

The ordering is as follows: meas(1, 1, wl1) meas(2, 1, wl1) ... meas(nbvm, 1, wl1) meas(1, 2, wl1) ... meas(nbvm, nip, wln)

Where meas can be: tau ... P

real(nsurf, nwl)pixel[ipix].meas[nwl].groundpar[nsurf]

ground parameters

This part can be ignored for now. The ordering is as follows: groundpar(1, wl1) groundpar(2, wl1) ... groundpar(nsurf, wl1) groundpar(1, wl2) ... groundpar(nsurf, wln)

real(nwl)pixel[ipix].meas[nwl].gaspar

gas absorption (tau gases) or molecular depolarization ratio. This parameter has to be provided only if the setting IFGAS (in the CELL HEADER) is set to 1.

The ordering is as follows: gaspar(wl1) gaspar(wl2) ... gaspar(wln)

integer(nip, nwl)pixel[ipix].meas[nwl].ifcov[nip]

ifcov (1 if a covariance matrix is available, 0 otherwise)

The ordering is as follows: ifcov(1, wl1) ifcov(2, wl1) ... ifcov(nip, wl1) ifcov(1, wl2) ... ifcov(nip, wln)

real(nbvm, nip, nwl)pixel[ipix].meas[nwl].cmtrx[nip][nbvm]

cmtrx (diagonal of covariance matrix, also known as Ω). These values have to be skipped if ifcov=0

The ordering is as follows: cmtrx(1, 1, wl1) cmtrx(2, 1, wl1) ... cmtrx(nbvm[c], 1, wl1) cmtrx(1, 2, wl2) ... cmtrx(nbvm[c], nip, wln)

integer(nip, nwl)pixel[ipix].meas[nwl].ifmp[nip]

ifmp (1 if a vertical profile (mprof) is available, 0 otherwise)

The ordering is as follows: ifmp(1, wl1) ifmp(2, wl1) ... ifmp(nip, wl1) ifmp(1, wl2) ... ifmp(nip, wln)

real(nbvm, nip, nwl)pixel[ipix].meas[nwl].mprof[nip][nbvm]

mprof (vertical profile of Rayleigh backscattering). These values have to be skipped if ifmp=0

The ordering is as follows: mprof(1, 1, wl1) mprof(2, 1, wl1) ... mprof(nbvm[d], 1, wl1) mprof(1, 2, wl2) ... mprof(nbvm[d], nip, wln)

[a] This fairly counter-intuitive coding has a reason: the cloud flag was at first intended to be a general processing flag (0 = pixel not to be processed, 1 = to be processed), cloud contamination is only one particular case. Now the flag is limited to cloud screening, but unfortunately the coding couldn't be changed right away. Since the framework is still in development, it is planned to correct this unnatural feature in the near future.

[b] These fields are actually not used by the processing and therefore the SDATA implementer is free to put whatever he or she likes here (e.g. 0 for non-gridded data). They are intended mainly for documentation and debugging. For satellite data, they make it possible to retrieve the pixel original information in the original database.

[c] nbvm is actually to be multiplied by ifcov(ip, iwl). If this last number equals 0, the array reduces to an empty set and no value is to be read.

[d] nbvm is actually to be multiplied by ifmp(ip, iwl). If this last number equals 0, the array reduces to an empty set and no value is to be read.


The field pixel[ipix].meas[nwl].meas_type[nip] of pixel structure is a special code which define the type of measure. Following table describes the valid codes and its interpretation:

Table 4.5. TYPES OF MEASUREMENTS

CONSTANT NAME (used in source code)VALUE (SDATA 2.0)MEANING
MEAS_TYPE_UNKNOWN 0 The measurement type is invalid or not yet implemented
MEAS_TYPE_TOD 11 Total Optical Depth
MEAS_TYPE_AOD 12 Aerosol Optical Depth
MEAS_TYPE_ABS 13 Aerosol absorption optical depth
MEAS_TYPE_P11 21 Phase Matrix Element P11
MEAS_TYPE_P12 22 Phase Matrix Element P12
MEAS_TYPE_P22 23 Phase Matrix Element P22
MEAS_TYPE_P33 24 Phase Matrix Element P33
MEAS_TYPE_P34 25 Phase Matrix Element P34
MEAS_TYPE_P44 26 Phase Matrix Element P44
MEAS_TYPE_LS 31 Lidar Signal
MEAS_TYPE_RL 32 Raman Lidar Signal
MEAS_TYPE_DP 35 Volume Depolarization Ratio
MEAS_TYPE_VEXT 36 Vertical Extinction profile
MEAS_TYPE_VBS 39 Vertical Backscatter profile
MEAS_TYPE_I 41 Normalized Radiance I [a]
MEAS_TYPE_Q 42 Polarized radiance Q[a]
MEAS_TYPE_U 43 Polarized Radiance U[a]
MEAS_TYPE_P 44 Polarization Rate: P = sqrt(Q*Q + U*U)/I

[a] All the Stokes Parameters are to be expressed as reduced quantities, without dimension

Equation 4.1. Conversion from absolute radiances to normalized, reduced radiances

I = radiance * π / E0


where radiance is the radiance of the instrument, and E0 the solar spectral flux, that may be both in mW / (m^2 * sr * nm) or equivalent units


4.2.2. Angle definition

This chapter main goal is to describe how the angles should be defined to be used inside of GRASP code. The universal spirit of GRASP, where many different instruments coexist (from satellite to ground based measurements) creates challenges to define a homogeneous way to define the angles keeping an unique geometry. As it is shown in the figure Figure 4.4, “Definition of GRASP geometry” GRASP angles are defined to be considered as "normal" for satellite reference.

Figure 4.4. Definition of GRASP geometry

Definition of GRASP geometry


Since GRASP angles are defined using a satellite reference, it provokes some problems to define what we could call as an intuitive "ground based" reference system. It is why we are going to put special emphasis on definition of the angles to these less intuitive applications. The intuitive reference for a "ground based" measurements, in spherical geometry, is given as follows:

  • θgb zenith angle: with the zero established in the zenith
  • ϕgb azimuth angle: with the zero considered in the sun position

where the sub index "gb" makes reference to "ground based".

The conversion to the GRASP geometry is done as follows:

Equation 4.2. Conversion from θgb (ground based) to θG (GRASP)

θG = 180° - θgb


Equation 4.3. Conversion from ϕgb (ground based) to ϕG (GRASP)

ϕG = 180° + ϕgb


Here we propose some examples for better understanding of the process. Before defining the measurement angles introduced in the code, both "intuitive" and "GRASP", we need first to consider the instrument viewing angle for each scenario (θv, ϕv). Following figure and table will provide some examples of angles defined for the ground based applications.

Figure 4.5. Ground based angles definition example

Ground based angles definition example


Table 4.6. Specific examples in ground based angle definition example

ExampleθsθvϕvθgbϕgbθGϕG
125°25°155°180°
225°12.5°12.5°167.5°180°
325°25°180°180°
425°37.5°12.5°180°167.5°
525°50°25°180°155°
625°90°65°180°115°
725°30°25°30°155°210°
825°90°25°90°155°270°

Since sunphptometers are widely used with GRASP following table provides information specifically about these instruments. Considering the instrument viewing angle for each scenario (θv, ϕv). They can be understood as the "movements of the motors". The process will be as follows:

instrument viewing angle -> angle in (intuitive) ground based -> angle in GRASP

Table 4.7. Sunphotometer angle description

Measure typeAngleInst. View.RangeGr. BasedRangeGRASPRange
Direct sunθθv = 0°[0°]θgb = θs[0° -- 90°][a]θG = 180° - θs[180° -- 90°][a]
ϕϕv = 0°[0°]ϕgb = 0°[0°]ϕG = 180°[180°]
Almucantarθθv = 0°[0°]θgb = θss]θG = 180° - θs[180° -- θs]
ϕϕv = 3°, 3.5°, 4°, 5°, ..., 90°, ... 180°[0° - 180°]ϕgb = 3°, 3.5°, 4°, 5°, ..., 90°, ... 180°[0° -- 180°]ϕG = 183°, 183.5°, 184°, 185°, ... 270°, ... 360°[180° -- 360°]
Principal plane measurement: Before the zenith θθv = -6°, ... , -3°, 3°, θmax < θs [-6° -- (θmax < θs)]θgb = θs + 6° ... θs + 3°, θs - 3°, ... θs - 6°, θs - θmax[(θs + 6°) -- 0][b]θG = 180° - θs - 6°, ... 180° - θs - 3°, 180° - θs + 3°, ... 180° - θs + θmax [(180° - θs - 6°) -- 180°]
ϕϕv = 0°[0°]ϕgb = 0°[0°]ϕG = 180°[180°]
Principal plane measurement: After the zenith θθv = θmin < θs , ... 140° or θmax - θs > 90°[(θmin < θs) -- 140°]θgb = θmin - θs ... 140° - θs or θmax - θs[0° -- 90°][b]θG = 180° + θs - θmin, ... 180° + θs - 140° or 180° + θs - θmax [180° -- 90°][b]
ϕϕv = 0°[180°]ϕgb = 180°[180°]ϕG = 0°[0°]
Polarized principal plane measurement: before the zenith[d]θ  θgb = 85°, 80°, 75°, ... 10°, 5°, 0° θG = 95°, 100°, 105°, ... 170°, 175°, 180° 
ϕ  ϕgb = 0° ϕG = 180° 
Polarized principal plane measurement: after the zenith[d]θ  θgb = 0°, 5°, 10°, ... 75°, 80°, 85° θG = 180°, 175°, 170°, ... 105°, 100°, 95° 
ϕ  ϕgb = 180° ϕG = 0° 

[a] θs refers to the solar zenith angle (for different measurements)

[b] decreasing values

[c] decreasing values

[d] The data of the polarized principal plane correspond always to fixed points in the sky and it is given for the instrument in so-called ground based coordinates.


In the case of nephelometer data angle definition is a bit different since only θ angle has to be defined, rest of the angles will be ignored. To provide nephelometer data (phase matrix) the conversion to the GRASP geometry is done as follows:

Equation 4.4. Conversion from θn (nephelometer scattering angle) to θG (GRASP)

θG = 180° - θn


4.2.3. How to prepare the photometer data

Sunphotometers are widely used with GRASP. They take measurements of sky radiance and direct sun. Many inversion strategies can be used to retrieve sunphotometer data but this section will explain how to define input data. The information that runs inside of GRASP has to be pre-pocessed in order to remove clouds an to calibrate and to normalize the data.

Following Section 4.2.1, “The SDATA format” the direct sun measurements can be described as AOD or TOD defined in Table 4.5, “TYPES OF MEASUREMENTS” as measurements of type 11 or 12. At this point it is needed to take into account that if "ifgas" field is defined as 1 in the case of AOD no gaseous absorption optical depth will be accounted but in the case of TOD they will be subtracted. The gases also affect the radiance measurements but in lower magnitude. In the case of TOD + radiances with ifgas=1, the same model will be applied to all measurements. If AOD is used some (minor) incongruences could come from the use of different models to calculate gases for AOD and for radiances.

Radiance measurements are defined with the constant MEAS_TYPE_I(41). Polarized measurements can be defined as Q,U (42, 43) or as polarization rate (44). It is also important to check how polarized data is going to be manipulated in the retrieval code based on inversion strategy defined in the settings file.

4.2.4. How to prepare the lidar data

Note that all processing will be considering range corrected profile for one wavelength. Procedure for other profiles from different wavelengths are exactly the same. Range corrected profile implies that at least background noise was subtracted and altitude correction were applied to the raw signal but if it's possible to consider all other corrections (electrical noise and overlap correction, dead time correction, gluing analog and photon-counting signals and all other your system may have) you should apply them.

Step 1. Background noise subtraction and range correction.

Let B' be the estimation of the background noise. Usually B' is estimated as P(ZB), accumulated and averaged around selected altitude ZBZ, much higher than the maximum altitude of lidar extraction in step 2 (Zmax). For example with maximum altitude Zmax=15km, it is averaged around 30 km and accumulated for the whole period of lidar observation. The noise and range corrected signal will be:

S(Zi) = (P(Zi)-B')*Zi2

Step 2. Altitude range selection and signal cropping.

The minimum Zmin and maximum Zmax altitudes are selected and the signal is cropped, so Zmin<Zi<Zmax. The minimum altitude should be selected as low as possible, preferably in the region where the overlap correction could be correctly applied. The Zmax should be selected from the following considerations: maximum altitude where the noise levels of lidar measurement are acceptable and the amount of atmospheric aerosols is still noticeable.

In case if the lidar used is inclined The sounding zenith angle should be provided for the corresponding wavelenghts. The angle value should be placed in the position of the SZA corresponding to this wavelength. All SZA's for all vertical profile measurements should be the same. Vertically pointed lidar should have this value set to 0. Note: Keep in mind that retrieved aerosol vertical profiles will be retrieved for vertically projected altitudes, i.e. if your lidar is inclined, the maximum altitude of the profile will be equivalent to Zmax*cos(SZA).

Step 3. Backscatter of molecular profile

Backscatter vertical profile is calculated inside of GRASP based on standard atmosphere model for each lidar wavelength[1]:

βmol(Zi,λ) = N * σ(λ) * (P(Zi)/PSA) * (TSA/T(Zi)/)

        where:

  • N molecular number density
  • σ(λ) total Rayleigh cross section per molecule which analytical formula can be written like σ(λ) = A * λ -B-C*λ-D/λ
  • Psa pressure of standard atmosphere model
  • Tsa temperature of standard atmosphere model
  • P(Zi) pressure profile of atmosphere
  • T(Zi) temperature profile of atmosphere

Step 4. Reducing the number of points in profiles

GRASP can use arbitrary altitude/range scale. However, to keep number of retrieved parameters reasonable and to fight higher noise contamination of lidar signals at higher altitudes it is recommended to use a logarithmical altitude/range scale with NZ points to represent aerosol profiles in the atmosphere. For that we have to present all vectors (altitude/range vector and profiles of lidar signals) in logarithmically equidistant manner.

1. Move to logarithmic scale and find altitude/range step:

logarithmic scale: Zi lg = lg(Zi)

step in log scale: ΔZ = Zmax lg - Zminlg/NZ

logarithmic altitude ranges (from hk to hk+1 ) for averaging data in logarithmically equidistant manner, k = 1 .. Nz: hk = Z0 lg + (k - 1) * ΔZ

2. Average the data profiles

Ak = (Σj=1n Aj(hk, hk+1) ) / n

        where:

  • A altitude vector or profile of lidar signal or molecular backscatter
  • n number of points inside logarithmic altitude ranges

After such procedure number of points in three main vectors reduced to NZ points in logarithmically equidistant manner. These values should be placed in places corresponding to the zenith wieving (vza or θv) angles for the vavelenths corresponding to lidar or vertical measurements.

Hint: the altitude/range vectors for measurements at all wavelengths have to be the same.

Step 5. Profile normalization

Values of lidar signals (except for volume depolarization profiles) vary from instrument to instrument, from detector to detector that is why GRASP requires normalized lidar signal. To be consistent to molecular optical depth inside code profile of molecular backscatter have to be normalized too. Normalized lidar and backscatter profiles:

A'k = Ak / ∫ZminZmax Ak dZ

        where A represents profile of lidar signal or molecular backscatter.

Caution: integration have to be done using meters in altitudes.

At the end for each wavelength you have to have normalized lidar profiles and altitude vector. Volume depolarization profiles don't need to be normalized, the only requirement for such observations is to be presented in the percentage range i.e. [1.0e-9, 100]

References: Anthony Bucholtz, Rayleigh-scattering calculations for the terrestrial atmosphere, Optical Society of America, 1995.

4.3. Output module

The output module is responsible for managing results for each segment or entire tile and driving them to the correspond destination (file or screen). There are two kinds of main output structure into GRASP: tile and segment. Segment output structure represents the output results obtained from the retrieval library. Then, core unit compact it and store it in a tile output structure, which contains the results of entire process.

Output module can be extended in the same way as input module. In the case of the output there are three kinds of extensions:

  • output segment function: it is called after retrieving a segment and is called with the results of that single segment.
  • output current function: it is called after retrieving each segment but it is called with the partial tile processed until that moment. In each call to this function it gets a partial tile nearest to completion. Last call to that function will send the entire tile results.
  • output tile function: at the end of the process a function is called sending to it the entire tile results. It can print a complete map of the process.

As the user can see, GRASP is very flexible in the way to work with the output. In a high optimized process it can be adapted directly to the format of the output database. By default, GRASP come with some ASCII output functions, which allow to print the results in a readable way (ASCII or specific GRASP format) into a file (using stream library) or on the screen. Additional extensions can be added to get the output in different formats such us HDF, NetCDF, png plots...

4.3.1. The list of GRASP output parameters

The output from GRASP is quite complex and strongly dependent on the settings used. As it was discussed in Section 4.1.2, “Retrieved characteristics”, there are two kinds of output products: direct and derived. Direct results are the values directly inverted in the initial guess array, then after obtaining them, forward model is called one more time to obtain the derived products. The list of direct and derived products obtained depends on the data and the inversion strategy selected (defined in the settings file). This chapter contain a complete list of products that can be obtained with GRASP but it does not mean that all of then can be obtained with all inversions strategies.

In the internal GRASP output structures there are some products that are repeated between direct products structure (array with the same shape as initial guess) and derived products. It is because that for some applications they are direct information, for others they are derived results. This depends on input data and inversion strategy defined in the settings file.

Before defining the output we are going to define the size of some arrays. This information is needed to understand the list of the products. For example, when a product such as AOD is defined as many times as wavelengths it is because the output will be wavelength dependent (a value for each wavelength):

  • NW: Number of wavelengths
  • NSD: Number of aerosol components
  • NKNOISE: Number of noises defined
  • NPARS: Number of parameters to be retrieved

The following list includes all products that can be obtained using GRASP:

  • Number of iterations (niter)
  • Total final measurement fitting residual for multi-pixel retrieval (rest)
  • Detailed (i.e., separated by the type of observation) final absolute measurement fitting residuals for segment (resat)
  • Detailed (i.e., separated by the type of observation) final relative measurement residuals for segment (resrt)
  • If SD (size distribution) retrieved in form of binned SD, the number of grid radii for SD (used to print output) (radius)
  • If SD retrieved in form of binned SD, grid radii
  • If the pre-calculated lognormal bins where used, the function describing each lognormal SD bin (used to print output) (SDL)
  • Main output for each single pixel (for both single- and multiple-pixel retrievals):
    • Single pixel total residual (meas. + smoothness constraints) (res)
    • Detailed absolute measurement residuals (resa[NKNOISE])
    • Detailed relative measurement residuals (resr[NKNOISE])
    • Retrieved aerosol and surface reflectance parameters (par[NPARS])
    • Angstrom exponent (Aexp)
    • For each wavelength:
      • Spectral total aerosol extinction (extt)
      • Spectral total aerosol single scattering albedo (ssat)
      • Spectral total aerosol absorption extinction (aext)
      If retrieved aerosol consist of several components:
      • Spectral extinction for each component (ext[NSD])
      • Spectral single scattering albedo for each component (ssa[NSD])
      • Real part of refractive index for each aerosol component (mreal[NSD])
      • Imaginary part of refractive index for each component(mimag[NSD])
  • If SD retrieved in form of binned SD, the optical properties can be calculated for fine and coarse modes, separated using chosen inflection radii:
    • Volume median radius (rv)
    • Standard deviation (std)
    • Concentration (cv)
    • Effective radius (reff)
    • Spectral extinction for each wavelength (ext[NW])
  • Phase matrix parameters:
    • Number of scattering angles (nangle)
    • Values of scattering angles (angles)
    • For each pixel and for all the wavelengths:
      • Phase matrix elements for each aerosol component ph11, ph12, ph22, ph33, ph34, ph44
      • Total aerosol phase matrix elements pht11, pht12, pht22, pht33, pht34, pht44
      • Lidar and depolarization ratios for each aerosol component (lr, dlpr)
      • Total aerosol lidar and depolarization ratios (lrt, dlprt)
  • If chemical composition is retrieved, then for each pixel and each aerosol component:
    • Relative humidity (rh(NSD))
    • Fraction of soluble (fslbl(NSD))
    • Insoluble fractions of soot (fsoot(NSD))
    • Insoluble fractions of iron (firon(NSD))
    • Insoluble fractions of “quartz” (fslbl(NSD))
    • Water fraction (fwtr(NSD))
  • Surface reflectance parameters for each pixel and for each wavelength:
    • All parameters of BRDF
    • All parameters of BPRDF
    • Surface Albedo
  • Lidar characteristics for each pixel:
    • Levels for vertical profiles
    • Vertical profiles for retrieved aerosol components (avp(NSD))
  • Lidar optical characteristics for each pixel and for each wavelength:
    • Extinction profiles (extv)
    • SSA profiles (ssav)
    • Lidar ratio profiles (lrv)
    • Lidar depolarization profiles (ldprv)
    • Retrieved lidar calibration coefficients (for lidar wavelength only)
  • Fit of every measured characteristic for each pixel and for each wavelength
  • Error estimation for each pixel:
    • Standard deviations of the random errors of the retrieved parameter logarithms (~relative errors) (ERRP)
    • Standard deviation of systematic errors of the retrieved parameter logarithms (BIASP)
    • Standard deviations of the random errors of the retrieved extinction for each aerosol component (~relative errors) (ERR_ext)
    • Standard deviations of systematic errors of the retrieved extinction for each aerosol component (BIAS_ext)
    • Standard deviations of the random errors of the retrieved total extinction (~relative errors) (ERR_extt)
    • Standard deviations of systematic errors of retrieved total extinction (BIAS_extt)
    • Standard deviations of the random errors of the retrieved single scattering albedo for each aerosol component (~relative errors) (ERR_ssa)
    • Standard deviations of systematic errors of retrieved single scattering albedo for each aerosol component (BIAS_ssa)
    • Standard deviations of the random errors of the of retrieved total single scattering albedo (~relative errors) (ERR_ssat)
    • Standard deviations of systematic errors of the retrieved total single scattering albedo (BIAS_ssat)
    • Standard deviations of the random errors of the of retrieved lidar ratio for each aerosol component (ERR_lr)
    • Standard deviations of systematic errors of the lidar ratio for each aerosol component (BIAS_lr)
    • Standard deviations of the random errors of the of retrieved depolarization ratio for each aerosol component (ERR_dr)
    • Standard deviations of systematic errors of the depolarization ratio for each aerosol component (BIAS_dr)
  • Radiative forcing for each pixel:
    • The heights for forcing output (HLV)
    • Broad band up-ward flux without aerosol at each height (BBUFX0)
    • Broad band down-ward flux without aerosol at each height (BBDFX0)
    • Broad band up-ward flux with aerosol at each height (BBUFXA)
    • Broad band down-ward flux with aerosol at each height (BBDFXA)
  • Estimations of aerosol particulate matter at the ground level (PM)
  • Aerosol type for each pixel (requires that the optical properties for fine and coarse modes are included in the calculated output)

4.4. Forward model

GRASP has several forward models and each of them are used (or not) depending on the application. For example, to retrieve nephelometer data just single scattering (particle properties) will be used. For other applications GRASP has also a multiple scattering module (radiative transfer) and a lidar signal module.

4.4.1. How to use the forward model: Derived products and reprocesing data

As it was explained in Section 4.1.2, “Retrieved characteristics” section, the retrieval algorithm works iteratively over an array of parameters (in its first definition is called the initial guess) until it represents the best solution. This solution array has same shape as initial guess (the same parameters and defined in the same position). Once it is obtained a final call of the forward model with the resulting array provides a complete list of output products.

For some applications it can be useful to use the forward model without inverting any data. It can be done easily in GRASP with the use of the setting parameter retrieval.convergence.stop_before_performing_retrieval=true. When no retrieval is performed, just one call of forward model is performed. If in the initial guess array the user has set an aerosol model, it will be used inside of the forward model obtaining therefore an entire output structure, with information in all fields.

This procedure can be used also to reprocess some data. If output parameters of a retrieval are stored then, they can be set as initial guess and then, running GRASP with the same settings except retrieval.convergence.stop_before_performing_retrieval=true the entire output can be obtained again. This procedure can be used to reprocess data with many objectives such us saving storage space (just save the output array of grasp and reprocess to obtain the rest if it is needed) or obtaining extra products in the future.

4.4.2. Synthetic data

The previous procedure can also help to simulate input data. In this case an aerosol model is set as initial guess and the code works just in forward model. As input data is necessary to set an SDATA file with a valid geometry. Then, using retrieval.debug.simulated_sdata_file parameter the user can set the path to the simulated filed. A SDATA file will be dump to that path where the geometry is the same and the measurements are filled with the output of forward model. Then, this SDATA file can be used to self-consistency tests, where synthetic data is retrieved.