Introduction

The Soil Moisture Ocean Salinity (SMOS) mission is a European Space Agency (ESA) earth observation satellite carrying the Microwave Imaging Radiometer using Aperture Synthetic (MIRAS) instrument. MIRAS is a passive 2-D interferometric radiometer operating at a frequency of 1.413 GHz (21 cm, L-band) and is used to measure the microwave emission from Earth’s surface as brightness temperature (TB), from which global maps of upper-level soil moisture and ocean salinity can be estimated [10].

The SMOS satellite follows a sun synchronous orbit (dawn/dusk at 6:00/18:00), and measurements are performed for a range of incidence angles between 0° and 60° across a swath of approximately 1000 km width with a spatial resolution ranging from 35 to 50 km depending on incidence angle. The MIRAS interferometer measures brightness temperatures in dual and full polarization (4 Stokes vector) mode.

The instrument data product comprises three level-1 sub-products:

• [L1A] calibrated visibilities between individual antenna receivers (i.e. the interferometric measurements from the sensor including the redundant visibilities), combined per integration time of 1.2 s (snapshot) in full polarization and prior the application of image reconstruction, provided in pole-to-pole (half orbit) time-based segments,
• [L1B] brightness temperature product as output of the image reconstruction of the SMOS observations, which consist of the Fourier components of brightness temperatures in the antenna polarization reference frame (XY polarization, hereafter),
• [L1C] multi-incidence angle brightness temperatures at the top of the atmosphere (XY polarization), geo-located in an equal-area grid. The L1C product comprises separate datasets available for ocean and land pixels with information available per pixel and per snapshot. An addition a so-called browse product contains the brightness temperatures averaged for a fixed incidence angle of 42.5°.

This document describes the analysis and processing of the L1C scientific data product, which consists of half-orbit daily archives with two separated physical files for each half-orbit: the product header file (HDR) and the data block binary file (DBL).

The Earth surface is viewed by SMOS in a snapshot-like configuration, which is restructured in the L1C product on an equal area Digital Global Grid (DGG) [13] with pixels based on the Icosahedral Snyder Equal Area (ISEA) projection at a spatial separation of about 15 km [11].

The SMOS L1C data product arranges the brightness temperatures sequentially for the real components of the XX, XY or YY polarizations, along with the same polarization sequence for the imaginary polarization components. The imaginary part is always zero for the single polarization components (e.g. XX or YY polarization). Figure 1 illustrates typical SMOS L1C half-orbit measurements stored in a DBL data file.

In order to use SMOS measurements in further studies for soil moisture or ocean salinity retrievals and for data assimilation, the L1C brightness temperature needs to be transformed from the antenna reference frame (XY) into the Earth’s surface reference frame (HV polarization, hereafter). This in done by considering the amount of rotation suffered by the polarization plane due to geometric and Faraday rotation.

Although the SMOS radiometer operates at L-band, which is a protected frequency band, observations are partly affected by Radio Frequency Interference (RFI) from various sources like radars and TV transmitters. Therefore RFI detection and mitigation for L1C brightness temperatures are important procedures, which need to be applied before brightness temperatures are used for retrievals of soil moisture and ocean salinity. The magnitude at which RFI affect the data, depends on the orbit (ascending or descending), polarization, and observation incidence angle. E.g., at higher incidence angles the pixel footprint is larger, which increases the impact of RFI compared to observations at nadir.

Some tools for polarization reference frame transformation and RFI mitigation have already been developed (e.g. SMOSview and SMOS-Box).

SMOSview is a JAVA based Graphical User Interface (GUI) software, that decodes SMOS products, displays their content graphically or numerically as tables and exports the data to alternative formats. SMOSview is not intended for a detailed analysis, visualization or processing of SMOS data but only for quick and easy look at the raw data products.

Similarly SMOS-Box is a module for the BEAM GUI software dedicated to the decoding, visualization and basic analysis of all operational data products (Levels 1 and 2) from the SMOS mission (see SMOS/Toolbox+BEAM website). SMOS-Box imports SMOS L1C and L2 science data products into the BEAM software and use the existing BEAM functionality to decode, display, manipulate, analyze and export SMOS data products to alternative formats in an interactive way or as batch processing. The BEAM software development has, however, been discontinued.

The Sentinel Application Platform (SNAP), the successor of BEAM, is a software developed as a common architecture for all Sentinel toolboxes (i.e. Sentinel-1 to -4 and SMOS, see http://step.esa.int/main/toolboxes/snap/). SNAP utilizes the SMOS-Box toolbox and allows to convert any SMOS data levels to NetCDF using its graph processing command line tool or a graphic user interface, also for specific regions of interest. SNAP also provides the Python module snappy, which allows the access to its SNAP Java API from Python; hence it can also work in batch processing mode.

The Centre Aval de Traitement des Donnes provides a Level 3 (L3) brightness temperature product in daily archives, which comprise daily global-coverage products that include all brightness temperatures in full polarization and transformed from the antenna frame (L1B) to the Earth’s surface polarization reference frame. The brightness temperatures (TBs) are binned and averaged for fixed angle classes. L3 TBs are already re-sampled on an Equal-Area Scalable Earth-grid (EASE-grid) [1] with 25 km resolution on a global scale; no detailed processing can be done in regions of interests at higher resolution on the L1C data product on the 15 km ISEA-grid.

Motivation

The discussed software tools mostly work with SMOS data in an interactive way with single data files at a time or via batch jobs. They load the whole data file and the user manually interacts with the software to glance through specific snapshots or regions for a given date. These tools are very helpful for a case-by-case analysis, but when a specific regions of interest need to be analyzed for a longer time period and large amounts of files must be processed limitation in term of reading and processing time becomes obvious. Since BEAM and SNAP are Java based they need to be interpreted by a Java Virtual Machine, which slows down its execution. The snappy module is a wrapper, which calls SNAP from Python and its performance, in terms of speed, might be even lower as expected from a scripting language. Therefore when many SMOS data files need to be fed to observation operators, e.g. for inversion algorithms or for data assimilation, speed-ups for reading and processing are paramount.

We developed the “Process_SMOSxL1C” mex-function as an alternative tool for rapid reading and processing of SMOS L1C data products with Matlab® [9] or GNU/Octave [4] programming languages. A mex-function is a C/C++ source code compiled in order to run within the Matlab workspace as a native Matlab function at a high execution speed. The “Process_SMOSxL1C” mex-function gives the user the ability to either work with a single file and exploit the powerful features of Matlab and Octave for analysis and visualization, or to run batch-jobs for many files and perform the same processing for a selected region and have the results as well as the original data automatically stored as MAT-files for future analysis. In addition it is possible to easily port the function to high performance computer (HPC) systems since most of them have C/C++ as default programming languages.

Implementation and architecture

The “Process_SMOSxL1C” function is implemented in a sequential structure. Since the SMOS L1C datablock (DBL) is a binary file, the data needs to be loaded sequentially and sorted to different variable types according to their specification [11]. The scientific L1C datablock consists of two datasets: one contains a list of snapshots and associated information in the swath (hereafter Swath Snapshot Information or SSI) and a second one contains the list of brightness temperature (TB) samples and associated geophysical information for each grid point sensed by the instrument (hereafter TSF for Temperatures Swath Full). In order to process the TBs, information from both datasets are needed, therefore SSI and TSF needs to be loaded sequentially before a region of interest is selected.

Every piece of information contained in a L1C datablock is sorted as a Matlab structure in their corresponding dataset of origin, either in SSI or TSF. Then the “Process_SMOSxL1C” mex-function returns the TSF structure with the original data contained in the DBL file and also gives the possibility to optionally return the SSI structure (in case of further re-processing).

The processing procedure follows five steps:

1. data filtering according to the RFI flags and other criteria,
2. interpolation of TB_XY values to obtain the four components of the Stokes vector at every time step,
3. transformation from the antenna reference frame (X- & Y-polarization) to the Earth’s surface reference frame (H- & V-polarization),
4. for every grid point the TBs are averaged for equally-spaced classes of incidence angles into 1° bins from 0 to 60°,
5. results are returned either to the Matlab workspace or saved as a MAT-file containing only the selected region of interest.

Step 1 uses the RFI flags contained in the TSF datablock for every grid point [6]. Additionally, other standard RFI detection criteria have been implemented as described by [6] and [15], namely:

with ST₄ the 4^th Stokes parameter, and $\overline{{\mathit{\text{TB}}}_{S1}}=0.5({\mathit{\text{TB}}}_X+{\mathit{\text{TB}}}_Y)$M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ \overline {T{B_{S1}}} \; = \;0.5\;(T{B_X}\; + \;T{B_Y}) \] \end{document} the mean value of the halved 1^st Stokes parameter for all observations in one pixel [6], and DTB_X the radiometric uncertainty for TB_X. Additional criteria (not strictly related to RFI) constrain out-of-range TB values by the requirement M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ 50\;K\; < \;\sqrt {T{B_X}^2\; + \;T{B_Y}^2} \; < \;500\;K \] \end{document} . Here the upper extreme may be related to RFI.

For the step 2, the TB values come in sequence for X-pol, Y-pol or XY, YX-pol at every snapshot time step. In order to have TB values with four Stokes parameters, an interpolation is needed in such a way that at every time step the TBs contain the [X, Y, Re{XY}, Im{XY}] Stokes parameters.

Step 3 transforms the polarization reference frame via a transformation matrix given by [6] as:

where the leftmost vector is A = [Re{TB_X}, Re{TB_Y}, 2Re{TB_XY}, –2Im{TB_XY}]^T for full polarization and A = [TB_X, TB_Y]^T for dual-polarization. The angle α = ω_Fa + α_ge is the sum of ω_Fa, the Faraday angle, and α_ge, the geometric rotation angle. Every matrix and vector element depends on the incidence angle of observation.

In order to obtain the Earth’s polarization reference frame vector, the above matrix system is solved for [TB_H, TB_V, ST₃, ST₄]^T at every incidence angle and for every pixel-point by using the GNU Scientific Library [7].

Step 4 re-samples the vector [TB_H, TB_V, ST₃, ST₄]^T at fixed incidence angle classes of 1° between 0 and 60° in the current version. This can be changed to other ranges and bin sizes, but requires re-compilation of the source code. The new re-sampled Stokes (or dual-pol) TBs are sorted in the following elements of the TSF structure:

> TSF.TB_Fixed_IncAngle,
> TSF.Fixed_IncAngle

where the first is an 3D array with the number of DGG pixel points in the first dimension (depending on the geographic region selected), the number of incidence angles (61 elements) in the second dimension and the polarization (either 4 or 2 elements) in the third dimension.

All information contained in the L1C datablock DBL is kept inside the structures TSF and SSI, with TSF having additional information related to the specific processing; therefore any step can be repeated or re-calculated within Matlab/Octave workspace. For a detailed description of the complete list of elements in the TSF and SSI structures, see the mex-function documentation at https://github.com/pablosaa/Process_SMOS_L1C/wiki.

The general usage of the “Process_SMOSxL1C” mex-function for a desired geographic region needs to be called from the Matlab or Octave workspace as:

> [TSF, SSI] = Process_SMOSxL1C
               (‘DBL_filename’,…
[lat_min, lon_min, lat_max, lon_max],…
  ‘PATH_TO_OUTPUT’);

with ‘DBL_filename’ a string with the absolute location of the L1C DBL data file, the second argument the limits of a latitude and longitude box for a geographic region of interest, and the last argument an optional string which indicates where the MAT-file containing the TSF and SSI variables will be saved. At least one output argument is mandatory, i.e. TSF, while SSI is optional.

Figure 2 shows the TB_H extracted and processed from the raw data file shown in Figure 1 for the Neckar catchment in south-west Germany, with the limits of the geographic region given by [47.5, 7.5, 50, 10] and assigned to the second argument of “Process_SMOSxL1C”.

See Table 1 for a reference benchmark regarding the processing time for a full archive with half-orbits as in Figure 1 and for a selected region as in Figure 2.

Quality control

The “Process_SMOSxL1C.cpp” code has been successfully compiled with GNU C++ version 4.7.1 in a Linux 64bit system. The compiled mex-function has been tested in three different Matlab versions: v7.8.0.347 (R2009a), v7.11.0 (R2010b), and v9.1.0.441655 (R2016b). Similarly, the code has been compiled and tested for Octave version 3.6.2. and 4.0.3.

The processing speed is listed in Table 1 for two case studies: a whole SMOS L1C data product containing one half-orbit (see Figure 1), and for a region of interest limited by its bottom-left corner [47.5°N, 7.5°E] and upper-right corner [50°N, 10°E] latitude and longitude, respectively (see Figure 2). This region contains the Neckar river catchment and is used for data assimilation studies by [12, 14] using of the “Process_SMOSxL1C” mex-function tool.

Table 1 includes in its 2^nd column the number of Digital Global Grid pixel points contained. The 3^rd column shows the result of only loading the raw L1C data product into the workspace by a Matlab m-script without using the mex-function, and the 4^th and 5^th columns indicate the total processing time by using the “Process_SMOSxL1C” mex-function. The last row in Table 1 indicates the percentage of minimum time needed in addition, when the option to save the processing outputs as a MAT-file is used.

Table 1 clearly shows the efficiency in terms of processing time by the mex-function, with Octave being slightly faster than Matlab. When the processed data is saved as MAT-file, Octave requires more time because all variables need to be translated to the MATIO library variable types [8].

For validation purposes the processed SMOS observations are “compared” against measurements by NASA’s Soil Moisture Active Passive (SMAP) mission launched in March 2015. SMAP carries a radiometer sensing microwave radiation at 1.4 GHz at one fixed incidence angle of 40° [5, 2]. SMAP has an antenna configuration different from SMOS, and SMAP’s L1C data products (36 km resolution) are delivered with the Earth’s polarization reference frame, i.e. no further processing is needed.

For this validation exercise the region corresponding to the Neckar catchment has been selected. Observations by SMOS and SMAP from April to December 2015 have been taken into account when both satellites have a morning overpass at the same day. Figure 3 shows the statistics of TBs at H-, and V-polarization after SMOS TB has been processed with the “Process_SMOSxL1C” mex-function. Note the good agreement despite the fact that the data come from two completely different instruments and methodologies.

Figure 4 shows the excellent agreement in dynamic range and daily variability between SMOS and SMAP TBs after processing the SMOS L1C product, which lends strong confidence to the capabilities of the software to adequately yield results free of artifacts which could jeopardize further applications of the data.

Compilation and testing

The best way to ensure the that “Process_SMOSxL1C” works for a particular system is by compiling the C++ source code on the used system. To that goal, the user needs first to assure, that Matlab’s mex compiler is properly configured and tested. Compilation procedures are indicated in the https://github.com/pablosaa/Process_SMOS_L1C/wiki. The user can also try the pre-compiled version provided and test whether or not the Process_SMOSxL1C mex-function is working. The easiest way to test is by trying to process a single file interactively from the Matlab or Octave workspace, following the simple command:

> TSF = Process_SMOSxL1C;

Immediately, a file browser pops-up for navigating thought directories to select a SMOS L1C DBL file. Once selected another dialogue-box pops-up (for the Matlab version) for entering the latitude and longitude limits for the processing region (these are done provided as command line inputs for the Octave version). If the Process_SMOSxL1C mex-function shows any message or warnings, it could mean that the input arguments are wrong, but indicates already that the mex-function is working in the system. When no warning appears the process should be working and the user can look at the content of the TSF structure.

Note that the “Process_SMOSxL1C” mex-function has as file extension either “.mex” or “.mexa64” for Octave and Matlab, respectively. “Process_SMOSxL1C” requires that the corresponding L1C DBL and HDR files are located in the same directory. Further description and examples on how to run the function and plot results can be found at https://github.com/pablosaa/Process_SMOS_L1C/wiki.

Operating system

GNU/Linux kernel 3.4.63-2.44, tested on OpenSuse version 12.2 or higher.

Programming language

Source code is written in GNU C++ language.

Additional system requirements

Minimum Memory RAM: 2 GB, disk space depending on number of half orbit L1C binary data files (each file is ~435 MB).

Dependencies

To compile the software the MEX library is needed, with two possibilities for compilation:

• as a MATLAB® function via MEX compiler, revision: 1.1.6.3 [9],
• as a GNU/Octave [4] mex function via MKOCTFILE compiler version 3.6.2 or higher.

The program requires:

• The GNU Scientific Library (GSL) version 2.2.1 for C and C++ language,
• The GNU C GSL Basic Linear Algebra Subprograms (GSLCBLAS) version 3.6.0. [7],
• MATIO library version 1.5.9 [8] is required for Octave to support MAT-file archiving.

List of contributors

Pablo Saavedra Garfias: project scientist, Meteorological Institute, University of Bonn.

Clemens Simmer: project Principal Investigator (PI), Meteorological Institute, University of Bonn.

Software location

Language

English