Software Tool for Soil Surface Parameters Retrieval from Fully Polarimetric Remotely Sensed SAR Data

The retrieval of soil surface parameters, in particular soil moisture and roughness, based on Synthetic Aperture Radar (SAR) data, has been the subject of a large number of studies, of which results are available in the scientific literature. However, although refined methods based on theoretical/analytical scattering models have been proposed and successfully applied in experimental studies, at the operative level very simple, empirical models with a number of adjustable parameters are usually employed. One of the reasons for this situation is that retrieval methods based on analytical scattering models are not easy to implement and to be employed by non-expert users. Related to this, commercially and freely available software tools for the processing of SAR data, although including routines for basic manipulation of polarimetric SAR data (e.g., coherency and covariance matrix calculation, Pauli decomposition, etc.), do not implement easy-to-use methods for surface parameter retrieval. In order to try to fill this gap, in this paper we present a user-friendly computer program for the retrieval of soil surface parameters from Polarimetric Synthetic Aperture Radar (PolSAR) imageries. The program evaluates soil permittivity, soil moisture and soil roughness based on the theoretical predictions of the electromagnetic scattering provided by the Polarimetric Two-Scale Model (PTSM) and the Polarimetric Two-Scale Two-Component Model (PTSTCM). In particular, nine different retrieval methodologies, whose applicability depends on both the used polarimetric data (dual- or full-pol) and the characteristics of the observed scene (e.g., on its topography and on its vegetation cover), as well as their implementation in the Interactive Data Language (IDL) platform, are discussed. One specific example from Germany’s Demmin test-site is presented in detail, in order to provide a first guide to the use of the tool. Obtained retrieval results are in agreement with what was expected according to the available literature.


Introduction
Electromagnetic (EM) scattering from the Earth surface depends on a number of parameters (i.e., permittivity, roughness, soil moisture content, vegetation biomass-index, etc.) that characterize the observed scene. Accordingly, large-scale and frequent monitoring of these quantities can be carried out exploiting remote sensing data, on condition that reliable estimation procedures are used [1].
Parameters estimation from remote sensing measurements is, however, everything but an easy task. It indeed requires, on one hand, to recognize and isolate the impact that each parameter has on the measured data and, on the other hand, to build up well-posed estimation procedures, where the number of unknowns does not exceed the available independent measurements (a wider discussion on this issue is reported for instance in [2]). The facet random tilt causes a random variation of the local incidence angle with respect to the radar look-angle , and a random rotation of the local incidence plane around the Line-Of-Sight (LOS) (see Figure 2), as stated by the following equations [19]: cos = cos + sin √1 + + tan = sin − cos (1) where and represent, respectively, the facets' slopes along the azimuth and range directions, which are modeled as independent Gaussian random processes. The facet random tilt causes a random variation of the local incidence angle ϑ l with respect to the radar look-angle ϑ, and a random rotation β of the local incidence plane around the Line-Of-Sight (LOS) (see Figure 2), as stated by the following equations [19]: where a and b represent, respectively, the facets' slopes along the azimuth and range directions, which are modeled as independent Gaussian random processes. Therefore, the incoherent first order backscattered electric field in the H (horizontal) and V (vertical) polarizations, namely = ( , ) , that arises from a generic randomly tilted facet can be expressed in terms of the incident field = , as: via the following scattering matrix [19,23]: where is the wave number, r is the radar-to-target distance, ( , ) is a polarizationindependent function depending on the small-scale roughness , while ( , ) and ( , ) represent the Bragg (SPM) or Fresnel (PO) coefficients for the H and V polarizations, with being the soil relative permittivity [19]. As Equation (3) shows, reciprocity (i.e., = ) is assumed for the PTSM data analysis.
According to Equations (2) and (3), the second order statistics of the EM field backscattered from the generic tilted facet can be obtained by averaging the quantities * with respect to the smallscale roughness, where each of the indexes p, q, u, w can stand for the polarization H or V.
Moreover, assuming that the facets' sizes are greater than both the EM wavelength and the correlation length of the small-scale roughness, but smaller than both the sensor geometric resolution and the correlation length of the large-scale roughness, it follows that the EM returns from different facets are uncorrelated. Accordingly, the second order statistics of the EM field backscattered from the whole surface can be obtained by averaging those of a single facet over and or, equivalently, due to Equation (1), over the slopes and .
In particular, the elements of the polarimetric covariance matrix, that is, the powers and correlations of the polarimetric channels for the whole scattering surface, are: where the operator 〈•〉 | represents the statistical mean with respect to the random variable .
Unfortunately, the statistical averages over and in Equation (4) cannot be computed in closed form [19,24]. However, assuming small values for the standard deviation of the facet slopes, a second order Taylor expansion of the quantities 〈 * 〉 | around the mean slopes can be considered, thus leading to the PTSM solution for the EM field. In this regard, it is worth recalling as: via the following scattering matrix [19,23]: S = S HH S HV S HV S VV = k 2 cos 2 ϑ l I ζ (ϑ l , β) where k is the wave number, r is the radar-to-target distance, I ζ (ϑ l , β) is a polarization-independent function depending on the small-scale roughness ζ, while F H (ε, ϑ l ) and F V (ε, ϑ l ) represent the Bragg (SPM) or Fresnel (PO) coefficients for the H and V polarizations, with ε being the soil relative permittivity [19]. As Equation (3) shows, reciprocity (i.e., S HV = S VH ) is assumed for the PTSM data analysis. According to Equations (2) and (3), the second order statistics of the EM field backscattered from the generic tilted facet can be obtained by averaging the quantities S pq S * uw with respect to the small-scale roughness, where each of the indexes p, q, u, w can stand for the polarization H or V.
Moreover, assuming that the facets' sizes are greater than both the EM wavelength and the correlation length of the small-scale roughness, but smaller than both the sensor geometric resolution and the correlation length of the large-scale roughness, it follows that the EM returns from different facets are uncorrelated. Accordingly, the second order statistics of the EM field backscattered from the whole surface can be obtained by averaging those of a single facet over ϑ l and β or, equivalently, due to Equation (1), over the slopes a and b.
In particular, the elements of the polarimetric covariance matrix, that is, the powers σ pq and correlations Γ pquw of the polarimetric channels for the whole scattering surface, are: where the operator · | ψ represents the statistical mean with respect to the random variable ψ. Unfortunately, the statistical averages over a and b in Equation (4) cannot be computed in closed form [19,24]. However, assuming small values for the standard deviation σ of the facet slopes, a second order Taylor expansion of the quantities S pq S * uw | ζ around the mean slopes can be considered, thus leading to the PTSM solution for the EM field. In this regard, it is worth recalling that the mean slopes are tied to the topography of the considered surface, which in turn can be obtained by exploiting Digital Elevation Models.
The full expressions of the PTSM solution for the quantities in Equation (4), as well as their validity limits and the direct model validation, can be found in [2,19,23,24]. Here, we just want to stress that the PTSM is the only analytical model employed in the literature for soil moisture retrieval from PolSAR data, for which the statistics of the scattered field and, in particular, its depolarization and cross-polarization, come straightly from the stochastic properties of the scattering surface.

Retrieval Method Based on Co-Pol, Cross-Pol Ratios and Co-Polar Correlation Coefficient
The PTSM expressions for the covariance matrix elements of previous section depend on a number of parameters, whether characterizing the employed radar system (e.g., the operating frequency and the radar look-angle), or describing the geometrical (e.g., the large-and small-scale roughness) and physical (e.g., the soil permittivity) properties of the scattering surface. Among these, the latter are typically unknown, but can be retrieved from polarimetric radar data by adopting a well-posed inversion procedure.
To this aim, in [19,23,24] it has been shown that the co-polar (CP) and cross-polar (XP) ratios, as well as the co-polar correlation coefficient (γ), namely: are substantially insensitive to the small-scale roughness since, for fixed values of the radar look-angle and operating frequency, they only depend on the large-scale roughness σ and on the soil relative permittivity ε. Accordingly, it is possible to build up numerical charts where co-pol and cross-pol loci or co-pol and correlation coefficient loci are drowned for different values of σ and ε. Therefore, these charts, an example of which is depicted in Figure 3, can be used to retrieve the soil permittivity (and then, the soil moisture content [23]) and the large-scale roughness from measured co-and cross-pol ratios, or from measured co-pol ratios and co-polar correlation coefficients, which can be obtained, for instance, using polarimetric SAR images. In this regard, it has been shown that for bare soil surfaces exhibiting an anisotropic roughness (i.e., different roughness characteristics along the azimuth and range directions, as is the case, for instance, for tilled soils), the co-pol/correlation method outperforms to co-pol/cross-pol one for the estimation of ε (using a mixing model [26] for soil moisture) and σ [27]. PTSM-based retrieval methods have been applied and/or validated for P-, L-and X-band SAR data sets acquired over bare soils under various radar look-angles, and for different surface roughness, and soil moisture values (see for instance [2,16,19,[23][24][25]). The best results are obtained at lower frequencies (P and L bands), which are less influenced by vegetation [16][17][18][19][20][21][22][23][24][25][26]. Sensors 2020, 20, x 6 of 20 Figure 3. Co-pol/Cross-pol chart (radar look-angle equal to 45°) based on the PTSM approach and its usage for the soil permittivity and large-scale roughness retrieval: measured co-pol and cross-pol ratios equal to 4 dB and −24 dB, respectively, lead to soil permittivity and soil roughness estimates equal to 4 and 0.09, respectively.

EM Scattering from a Slightly Vegetated Soil: The PTSTCM
What shown in Sections 2.1 and 2.2 allows us to analytically predict the EM scattering from a bare soil surface and to make use of these predictions to retrieve soil surface parameters from polarimetric radar data. However, in many situations of practical interest, radar data are acquired over vegetated soils. In these contexts, the surface backscattering contribution is mixed with the radar returns arising from the vegetation cover, which in turn can be modeled as the superposition of double-bounce (i.e., ground-trunks interactions) and volume (i.e., canopy) scattering contributions (see Figure 4) [20]. Therefore, in such scenarios, a single-component model for the prediction of the EM backscattering, and its usage for parameters estimation, is shown to be ill-suited [19]. On the other hand, accounting for all the above-mentioned scattering components for parameters inversion purposes would require some quantitative a-priori knowledge about the characteristics of the observed scene (e.g., the value of the standard deviation of the soil roughness), as to devise wellposed inversion procedures [18,25].  Co-pol/Cross-pol chart (radar look-angle equal to 45 • ) based on the PTSM approach and its usage for the soil permittivity and large-scale roughness retrieval: measured co-pol and cross-pol ratios equal to 4 dB and −24 dB, respectively, lead to soil permittivity and soil roughness estimates equal to 4 and 0.09, respectively.

EM Scattering from a Slightly Vegetated Soil: The PTSTCM
What shown in Sections 2.1 and 2.2 allows us to analytically predict the EM scattering from a bare soil surface and to make use of these predictions to retrieve soil surface parameters from polarimetric radar data. However, in many situations of practical interest, radar data are acquired over vegetated soils. In these contexts, the surface backscattering contribution is mixed with the radar returns arising from the vegetation cover, which in turn can be modeled as the superposition of double-bounce (i.e., ground-trunks interactions) and volume (i.e., canopy) scattering contributions (see Figure 4) [20]. Therefore, in such scenarios, a single-component model for the prediction of the EM backscattering, and its usage for parameters estimation, is shown to be ill-suited [19]. On the other hand, accounting for all the above-mentioned scattering components for parameters inversion purposes would require some quantitative a-priori knowledge about the characteristics of the observed scene (e.g., the value of the standard deviation of the soil roughness), as to devise well-posed inversion procedures [18,25].
Sensors 2020, 20, x 6 of 20 Figure 3. Co-pol/Cross-pol chart (radar look-angle equal to 45°) based on the PTSM approach and its usage for the soil permittivity and large-scale roughness retrieval: measured co-pol and cross-pol ratios equal to 4 dB and −24 dB, respectively, lead to soil permittivity and soil roughness estimates equal to 4 and 0.09, respectively.

EM Scattering from a Slightly Vegetated Soil: The PTSTCM
What shown in Sections 2.1 and 2.2 allows us to analytically predict the EM scattering from a bare soil surface and to make use of these predictions to retrieve soil surface parameters from polarimetric radar data. However, in many situations of practical interest, radar data are acquired over vegetated soils. In these contexts, the surface backscattering contribution is mixed with the radar returns arising from the vegetation cover, which in turn can be modeled as the superposition of double-bounce (i.e., ground-trunks interactions) and volume (i.e., canopy) scattering contributions (see Figure 4) [20]. Therefore, in such scenarios, a single-component model for the prediction of the EM backscattering, and its usage for parameters estimation, is shown to be ill-suited [19]. On the other hand, accounting for all the above-mentioned scattering components for parameters inversion purposes would require some quantitative a-priori knowledge about the characteristics of the observed scene (e.g., the value of the standard deviation of the soil roughness), as to devise wellposed inversion procedures [18,25].  However, in absence of any a-priori quantitative information, and assuming that the polarimetric radar data are acquired over slightly vegetated soils, where the ground-trunks interactions can be considered negligible, it is possible to resort to the PTSTCM [2]. This model is indeed founded on a two-component approach, where the scattered field is modeled as the superposition of independent surface and volume scattering contributions. In particular, the former is modeled by using the PTSM, via the scattering matrix and field statistics provided in Section 2.1. As for the latter, instead, it is modeled by a cloud of randomly oriented cylinder-like dipoles, whose scattering matrix is [20]: where ϕ is the dipoles' random orientation angle with respect to the vertical polarization direction, and α v is a function of parameters describing the dipole cloud. Accordingly, the second order statistics for the backscattering from the canopy layer can be expressed as: where f v = |α v | 2 , while the δ pq and ρ pquw coefficients only depend on the probability density function p(ϕ) of the dipoles' orientation ϕ. In particular, the value of these coefficients for three different probability density functions, corresponding to uniform dipoles orientation, prevalently horizontally distributed dipoles and prevalently vertically distributed dipoles, can be found in [2]. Finally, due to the independence between the surface and volume scattering mechanisms, the overall polarimetric powers σ tot pq and correlations Γ tot pquw relevant to the slightly vegetated soils can be expressed as the sum of the quantities in Equations (4) and (7).

Retrieval Method Based on the Modified Co-Pol Ratio and Co-Polar Correlation Coefficient
In Section 2.2, it has been shown that a proper combination of the elements of the covariance matrix in Equation (4) allows us to build up a well-posed retrieval procedure to get both the soil roughness and permittivity from polarimetric radar data acquired over a bare soil.
Similarly, in this section, we resume the methodology proposed in [2] to retrieve the soil roughness and permittivity from radar data acquired over a slightly vegetated soil. In particular, the following modified co-pol ratio and co-polar correlation coefficient: are only related to the surface parameters ε and σ, since the volumetric contribution f v cancels out. Accordingly, also in this case, it is possible to build up numerical charts where modified co-pol ratio and correlation coefficient loci are drowned for different values of σ and ε. These charts can then be used to retrieve the soil permittivity (and then the soil moisture content) and roughness from polarimetric radar data. The method has been tested with L-band SAR data acquired over natural fields characterized by different phenology, and during a period covering all the phases of the vegetation growth. From such an analysis, it turned out that, as already mentioned, in presence of low vegetation (vegetation height lower than 50 cm, or cross-polarized ratio smaller than 0.1, and negligible double-bounce component), the PTSTCM soil moisture estimates were in better agreement with in situ measurements than the ones provided by the retrieval method based on the IGHD [2].

The Hallikainen's Mixing Model of Volumetric Soil Moisture Retrieval from Permittivity
The Hallikainen's mixing model has been employed to retrieve volumetric soil moisture from permittivity [26]. This mixing model was developed exploiting data of five soil types with wide textural (geotechnical) compositions, and employing nominal frequencies of 1.4,4,6,8,10,12,14,16 and 18 GHz [26]. In the Hallikainen's mixing model, the real and imaginary parts of the (complex) permittivity are second order polynomial functions of the volumetric soil moisture, whose coefficients depend on the sand and clay textural components of the soil. The general form of the Hallikainen's mixing equation reads: where SAND and CLAY are sand and clay textural components in percent by weight, and a i , b i , and c i are the function of the (given) frequency, with their values are given in [26].

Software Design and Implementation
PTSM-OOP from the University of Naples Federico II is a software tool developed in the IDL environment, aimed at retrieving soil surface parameters from polarimetric SAR data. PTSM-OOP is an open source project available for public usage via [28] under course material.
We chose IDL for several reasons. First of all, it allows to easily manage large images, such as remote sensing SAR ones. In addition, with almost no required change in the source code, it can be run both on Windows and Linux systems. Finally, IDL is the basic language behind ENVI, so that our software can be easily employed in conjunction with ENVI.
A freely available IDL virtual machine is needed to run the software. The overall PTSM-OOP workflow is depicted in the block scheme of Figure 5 and it will be described in the following.
Sensors 2020, 20, x 9 of 20 3), and the best fit pair of above-mentioned parameters is selected. Finally, the retrieved dielectric constant is converted into volumetric soil moisture using the Hallikainen model (Section 2.5). This approach can be selected if the analyzed scenario is a bare soil in a flat area.  The user provides the PolSAR data to the software and gets as output the maps of soil moisture, soil permittivity and large-scale roughness. The IDL compiler checks the source codes that are selected by the user, and then applies the proper algorithm. PTSM-OOP Graphical User Interface (GUI) has menu entries in the main bar (see top of Figure 6) to access all application's functionality. The user can read PolSAR data channels (i.e., HH, VV and HV complex data), ancillary file (including line and column numbers and spacings), and incidence-angle maps (that are particularly useful for scenarios with significant topography and/or for airborne SAR acquisitions, where the radar look angle significantly changes within the radar swath) in the "file" menu. Further details on input data will be provided in Section 4.   The "exit" menu item of the "file" menu allows to exit from the application. After reading the PolSAR data and the ancillary information in the "file" menu, the next step is the multi-look analysis. Here, the user can select between two options, namely "Whole image multi-looking mode" and "ROI image multi-looking mode," depending on whether the entire image or a selected Region of Interest (ROI) has to be considered. "Multi-looking" operation computes powers of the three polarimetric channels by averaging the square modulus of each channel over rectangular windows of size 10 by n, where n is such that the final pixel is approximately square. In addition, the copol correlation is computed by averaging in the same way the product of HH by the complex conjugate of VV. These operations also allow to reduce the effect of speckle. After this, it follows the core of the software tool, namely the "Processing" step. Here, the user can select among nine possible modes from the menu, according to the desired approach for the soil moisture and roughness retrieval. These nine approaches are described in the following section.

Co-and Cross-Polarized Ratios (Co-Pol/Cross-Pol) Method
If this option is selected, the retrieval method described in [19], and briefly recalled in Section 2.2, is employed. From the input measured SAR polarimetric data, σ HH , σ HV , σ VV are evaluated and the co-polarized and cross-polarized ratios (i.e., Equation (5)) are computed. For each (multi-look) image pixel, they are then compared to the charts of co-polarized and cross-polarized ratios computed from the model as a function of dielectric constant ε and large-scale rms slope σ (Figure 3), and the best fit pair of above-mentioned parameters is selected. Finally, the retrieved dielectric constant is converted into volumetric soil moisture m V using the Hallikainen model (Section 2.5). This approach can be selected if the analyzed scenario is a bare soil in a flat area.

Co-Polarized Ratio and Correlation Coefficient (Co-Pol/Correlation) Method
In this option, the initial settings and retrieval rationale are the same as those of Section 3.1, but the co-polarized ratio and the co-polarized correlation coefficient are exploited for the inversion, according to the approach described in [24]. Just as the previous one, this approach can be safely used only if the analyzed scenario is a bare soil in a flat area. However, it has the advantage that only co-polarized channels HH and VV are needed (dual-pol). In addition, recent studies [27] show that this approach is more robust to surface roughness anisotropy than the previous one. Therefore, it must be preferred if a tilled bare soil is considered.

Modified Co-Polarized Ratio and Correlation Coefficient (Modified Co-Pol/Correlation) Method
This option is founded on PTSTCM approach summarized in Sections 2.3 and 2.4, in which the volumetric scattering contribution from the (moderate) vegetation cover is described via a dipole-cloud model with a uniform distribution of the dipole orientation angle [2]. This option can be selected if the considered scenario is a moderately vegetated soil in a flat area, and no specific a priori information on vegetation type is available.

Co-Pol/Correlation Method with DEM
This approach is basically the same as the one of Section 3.2, extended to account for the local surface non-zero mean slopes related to macroscopic topography, according to the approach described in [23,24]. In this case, a Digital Elevation Model (DEM) must also be provided as input. The latter is used for the calculation of the local (i.e., at the pixel scale) mean slopes. This approach can be selected if the analyzed scenario is a bare soil in the presence of non-negligible topography.

Modified Co-Pol/Correlation Method with DEM
This approach is basically the same as the one of Section 3.3, extended to account for the local surface non-zero mean slopes related to macroscopic topography, according to the approach described in [23,24]. In this case, also a DEM must be provided as input, see Section 3.4. This approach can be selected if the analyzed scenario is a moderately vegetated soil in the presence of non-negligible topography, and no specific a priori information on vegetation type is available.

Modified Co-Pol/Correlation Method-Horizontal
This approach is basically the same as the one of Section 3.3, but in this case the volumetric scattering contribution from (moderate) vegetation is described via a dipole-cloud model with prevalently horizontal dipole orientation angle [2]. This option can be selected if the considered scenario is a moderately vegetated soil in a flat area, and available a priori information on vegetation type leads to assume a prevalently horizontal distribution of leaves and/or branches.

Modified Co-Pol/Correlation Method-Vertical
This approach is the same as the one in Section 3.6, but with a dipole-cloud model with prevalently vertical dipole orientation angle.

Modified Co-Pol/Correlation Method with DEM-Horizontal
This option is the same as the one of Section 3.6, but it requires an input DEM to account for macroscopic topography. This approach can be selected if the analyzed scenario is a moderately vegetated soil in the presence of non-negligible topography, and available a priori information on vegetation type leads to assume a prevalently horizontal distribution of leaves and/or branches.

Program Functionality
GUI is a point-and-click way of running an application in a windowed environment. Complex and repetitive tasks can be simplified by an interactive GUI, which hides the function and implementation details and allows the user to focus on the data, rather than on the mechanism of entering IDL commands to the program itself. PTSM-OOP GUI environment is provided in IDL programming language. Another advantage of PTSM-OOP GUI toolkit is the ability to be run on different platforms, without changing the appearance of the mapping GUI.
PTSM-OOP GUI application contains a number of components known as widgets, each of which having a specific appearance. In particular, PTSM-OOP GUI has the following widgets: 1. A menu bar (main bar), which allows the configuration of application properties. 2. Scrolling lists, which enable the selection of input data, processing method and output map to display (Figure 7). 3. IDL graphical windows displaying the output imageries.

Program Functionality
GUI is a point-and-click way of running an application in a windowed environment. Complex and repetitive tasks can be simplified by an interactive GUI, which hides the function and implementation details and allows the user to focus on the data, rather than on the mechanism of entering IDL commands to the program itself. PTSM-OOP GUI environment is provided in IDL programming language. Another advantage of PTSM-OOP GUI toolkit is the ability to be run on different platforms, without changing the appearance of the mapping GUI.
PTSM-OOP GUI application contains a number of components known as widgets, each of which having a specific appearance. In particular, PTSM-OOP GUI has the following widgets: 1. A menu bar (main bar), which allows the configuration of application properties. 2. Scrolling lists, which enable the selection of input data, processing method and output map to display (Figure 7). 3. IDL graphical windows displaying the output imageries.
PTSM-OOP application responds to any interaction with each widget in a unique way. When user's interaction with the PTSM-OOP GUI is detected, the PTSM-OOP GUI takes particular course of action depending on the selected widgets. The main steps of the input procedure are guided, and the catch error function provides the control on the possible user errors. The execution time for each PTSM-OOP run depends on the size of the input images (i.e., HH, VV and HV layers).  PTSM-OOP application responds to any interaction with each widget in a unique way. When user's interaction with the PTSM-OOP GUI is detected, the PTSM-OOP GUI takes particular course of action depending on the selected widgets. The main steps of the input procedure are guided, and the catch error function provides the control on the possible user errors. The execution time for each PTSM-OOP run depends on the size of the input images (i.e., HH, VV and HV layers).
The PTSM-OOP's functions are listed below: 1. A base widget named PTSM which has the form of widget_base(). This widget contains all the widgets subsequently created inside the GUI. This type of the widget is known as top-level base widget, which is a parent to all the child widgets contained in the PTSM-OOP GUI. 2. A widget_label(), creates a label widget and returns the widget identifier label. 3. Number of widget_button(), have control of the GUIs buttons. 4. Finally, a widget_control(), which is called to render on screen the widget hierarchy controlled by base_ widget().
Widget buttons needed to read input data are under the "file" button. Before mentioning them, it is useful to describe the required format of input data for the subsequent "multi-looking" step. Each polarimetric channel (which is a complex image) must be a two-dimensional array of complex numbers, stored in an *.dat file. The number of lines and columns, and line and column spacings, must be provided as a simple text file (ancillary data). Finally, a two-dimensional array of real numbers, with the same size of SAR data, containing the incidence angle map, stored in an *.dat file (incidence-angle map file), can be optionally used. Therefore, a preprocessing step aimed at converting the input PolSAR data in this format is required. Of course, this step strongly depends on the storage format of the considered dataset: the PTSM-OOP software can perform this preprocessing step, with no need for external tools, when the original SAR data are in *.h5 (HDF) format or if they are E-SAR data of the AGRISAR mission [29]. Otherwise, the user can use the freely available ESA SNAP package to generate HDF format data from original SAR data of the main polarimetric sensors, such as ALOS, RADARSAT, etc. Alternatively, preprocessing can be carried out with most of the licensed (e.g., ENVI, IDL, Matlab) or free (e.g., SNAP, PolSAR PRO) packages that are usually employed to process SAR data. For instance, since IDL is connected to ENVI software, all the formats that can be read by ENVI can be transformed to the required *.dat format via Export to IDL Variable, which is a built-in function in ENVI.
In conclusion, three widget_buttons read the input data under the "file" -> reading data button. They are named: Read_the_HH_data, Read_the_HV_data, Read_the_VV_data, Read_the_incidence angle map, Read_the_ancillary data and PTSM_exit.
In the "Multilook analysis" button, two widget_buttons, which are named Whole_image _multilooking_mode and ROI_image_multilooking_mode, are handling the Multilooking step.
In the "Processing" button, nine widget_buttons allow the selection among the processing methods listed in Figure 5 and described in Section 3.
In each of the remaining buttons (Iplot, Icontour, Imap, Iimage, ENVI, Help, About), only one widget_button is acting. They are designed to call functions from ENVI software, and to display HELP/About information.

Program Run Example
In this section we illustrate a step by step example procedure of PTSM-OOP, to obtain meaningful soil parameters from fully polarimetric SAR images. The employed SAR data were acquired in the framework of the 2006 AgriSAR campaign (Figure 8) [29]. The campaign was funded by the European Space Agency (ESA) with the main goal of collecting in-situ measurements, microwave and optical remote sensing data over an intensive farming region in Demmin (Germany), during an entire vegetation growth period. DEMMIN (Durable Environmental Multidisciplinary Monitoring Network) is a well-known test site located in Mecklenburg-Western Pomerania North-East Germany, approximately 60 km north of Neustrelitz and 150 km north of Berlin (Figure 8). The general information about the E-SAR data are given in Table 1  The sensor operated at three different bands/channels (X, C and L), and it was mounted onboard a Dornier Do-228 aircraft, operating in a range of flight altitudes between 300 m and 6000 m above ground. For each flight the sensor was able to collect data at a single center frequency (i.e., no  The sensor operated at three different bands/channels (X, C and L), and it was mounted onboard a Dornier Do-228 aircraft, operating in a range of flight altitudes between 300 m and 6000 m above ground.
For each flight the sensor was able to collect data at a single center frequency (i.e., no simultaneous multi-bands acquisitions were possible), either with one, two, or four polarization channels. However, the system was enabled to get fully polarimetric acquisitions only at L-band. Data were freely delivered to the scientific community at two processing levels, namely the Radar Geometry Image product (RGI) and the Geo-coded and Terrain-Corrected product (GTC). RGI (SLC) products exhibits a slant range resolution equal to 2 m and an azimuth resolution equal to {0.9 m, 0.9 m, 1 m} (for X-, C-and L-band, respectively), while a 2 m × 2 m grid in WGS-84 UTM projection, zone 33, has been used for GTC products [29].
In particular, in this study retrieval of soil surface parameters via PTSM-OOP algorithm was applied on the available Single Look Complex (SLC) geo-coded L-band quad-polarimetric GTC SAR images.
The study area is characterized by several crop types, such as: sugar beet, wheat, winter barley, winter rape, grassland and corn. Optical and SAR images of the area are shown in Figure 9. A wide range of ground truth data collected at the same time of the radar acquisitions are available: vegetation phenology, terrain conditions, precipitations, and volumetric soil moisture. In particular, the last was measured with different techniques (i.e., time-domain reflectometry, gravimetric and capacitive measurements) and different time-sampling scenarios (intensive campaigns over many fields, see Figure 9a, weekly measures on a limited set of fields and via continuous measurements stations over few fields) [29]. All details on the in-situ measurement campaign are reported in [29].
Sensors 2020, 20, x 14 of 20 simultaneous multi-bands acquisitions were possible), either with one, two, or four polarization channels. However, the system was enabled to get fully polarimetric acquisitions only at L-band. Data were freely delivered to the scientific community at two processing levels, namely the Radar Geometry Image product (RGI) and the Geo-coded and Terrain-Corrected product (GTC). RGI (SLC) products exhibits a slant range resolution equal to 2 m and an azimuth resolution equal to {0.9 m, 0.9 m, 1 m} (for X-, C-and L-band, respectively), while a 2 m × 2 m grid in WGS-84 UTM projection, zone 33, has been used for GTC products [29].
In particular, in this study retrieval of soil surface parameters via PTSM-OOP algorithm was applied on the available Single Look Complex (SLC) geo-coded L-band quad-polarimetric GTC SAR images.
The study area is characterized by several crop types, such as: sugar beet, wheat, winter barley, winter rape, grassland and corn. Optical and SAR images of the area are shown in Figure 9. A wide range of ground truth data collected at the same time of the radar acquisitions are available: vegetation phenology, terrain conditions, precipitations, and volumetric soil moisture. In particular, the last was measured with different techniques (i.e., time-domain reflectometry, gravimetric and capacitive measurements) and different time-sampling scenarios (intensive campaigns over many fields, see Figure 9a, weekly measures on a limited set of fields and via continuous measurements stations over few fields) [29]. All details on the in-situ measurement campaign are reported in [29].  AGRISAR data (including both SAR data and ground truth) can be received free of charge from ESA by following the procedure described in [30].

Read the Input Data
The PTSM-OOP procedure is started by reading the input PolSAR data including HH, HV and VV SAR channels/bands, incidence angle map of the study area and ancillary data of the available PolSAR images, via the "file" button of the PTSM-OOP program (Figure 7, top-left). The used data is named "06agrsar0609 × 1_t01" [29], and is part of 2006 AgriSAR campaign. "i06agrsar0609 × AGRISAR data (including both SAR data and ground truth) can be received free of charge from ESA by following the procedure described in [30].

Read the Input Data
The PTSM-OOP procedure is started by reading the input PolSAR data including HH, HV and VV SAR channels/bands, incidence angle map of the study area and ancillary data of he available PolSAR images, via the "file" button of the PTSM-OOP program (Figure 7, top-left). The used data is named "06agrsar0609 × 1_t01" [29], and is part of 2006 AgriSAR campaign. "i06agrsar0609 × 1_ch1_t01_int_slc_geo.dat" (HH), "i06agrsar0609x1_ch2_t01_int_slc_geo.dat" (HV) and "i06agrsar0609x1_ch3_t01_int_slc_geo.dat" (VV) are the three PolSAR bands of the Demmin study area. "incmap06agrsar0609x1_t01.dat" is the given incidence angle map; "06agrsar0609x1_t01_README_GTC" is the ancillary data, which provides some information about the E-SAR sensor such as description of the campaign, SAR processing parameters file, repeat pass master identifier, DEM slant range, range pixel positions, pixel spacing north, minimum/maximum easting/northing etc. The E-SAR data are in the complex format. The processing parameters of E-SAR data are given in Table 1, including original pixel spacing; we considered a further multi-looking to the entire image in the software run, as described below.

Multi-Looking the Input Data
After reading the data to PTSM-OOP software, we multi-look the given PolSAR imageries as has been shown in Figure 7, top-right. In this step, powers and correlation are also computed according to Equation (4) by averaging over rectangular windows, and, accordingly, measured co-polar and cross-polar ratios, and co-polar correlation coefficients are computed via Equation (5) [19,24]. Multi-looking can be carried out in two different sections, depending on the necessity: (a) whole image multi-looking mode and (b) ROI multi-looking mode. In the given example, we use whole image multi-looking mode, and the employed window size is 10 by 10 pixels, so that the final map resolution is 20 m by 20 m. The incidence angle data also is resampled accordingly for further analysis. Depending on the necessities, ROI multi-looking mode also is available, in which the operator can select the regions which need to be processed.

Processing
After multi-looking analysis, we proceed to the processing step. The nine available processing modules are presented in the Section 3 and Figure 5. In this example, we use the "Co-and cross-polarized ratios (co-pol/cross-pol)" method (the first choice in the scrolling list) to retrieve the soil moisture parameters.
Theoretical values of Normal Radar Cross Sections (NRCS) for HH, VV and HV polarizations and of the HHVV correlation are computed for a set of values of ε and σ, as described briefly in Section 2.1, see Equation (4), and in more detail in [19] and [24]. Co-polar and cross-polar ratios are then computed via Equation (5). This results in a look-up table, which associates a co-polar and cross-polar ratios pair to each ε, σ pair (a graphical example of such a look-up table is provided by Figure 3).
The final step is the comparison of the measured co-polar and cross-polar ratios (i.e., those obtained from SAR data and computed in the multi-looking step) with the just-produced lookup tables to retrieve soil permittivity and soil roughness for each pixel.
The other eight processing modules follow analogous procedures, except that for methods of Sections 3.4, 3.5, 3.8 and 3.9, a file containing the surface DEM is requested as further input (for the considered case the surface is substantially flat and a DEM is not needed). This processing phase is completely unsupervised and no action is required to the user.

Results
The only remaining task is the presentation of the results and calculation of (percentage of) water content for each given pixel. In Figures 10 and 11, dielectric constant, ε, and large-scale roughness σ retrieval maps for Demmin study area for 16 May 2006, with co-pol/cross-pol method, are given. The color scale is showing the strength of the retrieval parameters (i.e., ε and σ), and the black pixels

Comparison with Ground Truth Data
For the 16 May 2006 images, the vegetation average height is low, and, as stressed before, for low vegetation covers, PTSM (and/or PTSTCM) can be safely used [2]. In particular, by using the copol/cross-pol method, a root mean square error (RMSE) of 14 vol% and a percentage of invalid pixels of 31% are obtained; by using the modified copol/correlation method, RMSE and percentage of invalid pixels decrease to 9 vol% and 8%, respectively. These results are in agreement with what is expected considering the results in [2]. In both cases, the entire processing chain requires about 3 min for image size of about 4300 × 1600 pixels per polarimetric channel, on a laptop with an Intel I7 processor (a few minutes more are needed for ALOS and RADARSAT whole images).
A complete comparison of PTSM and PTSTCM results with ground truth data is available in the literature [2,16,19,[23][24][25]27,28]. However, to further verify the correct implementation of these methods in the presented software, we carried out a similar analysis for the results of the software. In particular, we performed the same retrieval procedure described in the previous subsections for the polSAR acquisitions of 19 April 2006, 3 May 2006and 11 May 2006, in order to consider the time period from 19 April to 16 May 2006, during which the vegetation was mostly in an early stage of growth and its average height was low (ranging from nearly zero, for corn fields, to more than 1 m, for rape fields [29]).
For each of these acquisitions, the overall results are similar to the ones obtained for the 16 May one (except for the 11 May acquisition, for which the results are slightly worse). In order to perform a more punctual comparison, we considered separately the soil moisture results for different fields. We underline that both SAR estimates and in situ measurements, relevant to each single field and for each single acquisition, show large standard deviations, often of the same order of mean values. Therefore, for each field and at each acquisition date, we compare average values of retrieved and in situ measured soil moisture. Results are presented in the scatterplots of Figure 13 for the copol/correlation method (a) and the modified co-pol/correlation method (b).

Comparison with Ground Truth Data
For the 16 May 2006 images, the vegetation average height is low, and, as stressed before, for low vegetation covers, PTSM (and/or PTSTCM) can be safely used [2]. In particular, by using the co-pol/cross-pol method, a root mean square error (RMSE) of 14 vol% and a percentage of invalid pixels of 31% are obtained; by using the modified copol/correlation method, RMSE and percentage of invalid pixels decrease to 9 vol% and 8%, respectively. These results are in agreement with what is expected considering the results in [2]. In both cases, the entire processing chain requires about 3 min for image size of about 4300 × 1600 pixels per polarimetric channel, on a laptop with an Intel I7 processor (a few minutes more are needed for ALOS and RADARSAT whole images).
A complete comparison of PTSM and PTSTCM results with ground truth data is available in the literature [2,16,19,[23][24][25]27,28]. However, to further verify the correct implementation of these methods in the presented software, we carried out a similar analysis for the results of the software. In particular, we performed the same retrieval procedure described in the previous subsections for the polSAR acquisitions of 19 April 2006, 3 May 2006 and 11 May 2006, in order to consider the time period from 19 April to 16 May 2006, during which the vegetation was mostly in an early stage of growth and its average height was low (ranging from nearly zero, for corn fields, to more than 1 m, for rape fields [29]).
For each of these acquisitions, the overall results are similar to the ones obtained for the 16 May one (except for the 11 May acquisition, for which the results are slightly worse). In order to perform a more punctual comparison, we considered separately the soil moisture results for different fields. We underline that both SAR estimates and in situ measurements, relevant to each single field and for each single acquisition, show large standard deviations, often of the same order of mean values. Therefore, for each field and at each acquisition date, we compare average values of retrieved and in situ measured soil moisture. Results are presented in the scatterplots of Figure 13 for the co-pol/correlation method (a) and the modified co-pol/correlation method (b). A visual analysis of these scatterplots shows that both methods provide reasonable results, but in the first one for some fields and acquisitions, a strong overestimate is obtained, whereas the second one, accounting for the effect of the (moderate) vegetation, is in better agreement with the ground truth. This is confirmed quantitatively by evaluating the correlation coefficient and RMSE for the two cases. In fact, for the standard co-pol/correlation method, the correlation coefficient is rather low, i.e., 0.34, and RMSE is 13.6 vol%, whereas for the modified co-pol/correlation method the correlation coefficient is 0.64, and RMSE is 9.1 vol%. In addition, these results are practically coincident with those obtained in [2] for the same considered time period, thus confirming the correct implementation of the methods in the software tool.

Conclusions
In this paper, a software tool for the automatic retrieval of soil surface parameters from SAR data and via PTSM and PTSTCM has been presented. The available methodologies for soil surface parameters retrieval have been presented, and the basic theory behind PTSM/PTSTCM has been briefly but thoroughly reviewed.
The presented software tool allows to create maps of dielectric constant , large-scale roughness , and volumetric soil moisture from fully polarimetric and/or dual polarimetric SAR images. In fact, PTSM can be also employed in the absence of cross-polarization (HV), if a bare soil is imaged. In this case, only co-pol ratios and HH-VV correlations are needed. Availability of fully polarimetric data allows obtaining realistic retrieved values also in moderately vegetated areas. The software tool implements nine different algorithms, to be selected according to the typology of the observed scenario, as specified in Section 3. Thanks to its user-friendly interface, the tool is suitable for both expert and non-expert users. With this respect, the presented tool bridges a gap in the available software tools for PolSAR processing: in fact, commercially and freely available software tools for the processing of SAR data that include routines for basic manipulation of polarimetric SAR data (e.g., coherency and covariance matrix calculation, Pauli decomposition, etc.) certainly exist, see Section 1, but they do not implement easy-to-use methods for surface parameter retrieval. This has been the main motivation of our work.
In order to provide a first guide to the use of the tool, its application to a specific example from Germany's Demmin test-site, for which ground truth is available, has been presented. Finally, a comparison of retrieved soil moisture with in-situ measurements has been illustrated, showing a good agreement, at least for the method accounting for the presence of vegetation (modified copol/correlation). Obtained retrieval results are in agreement with what was expected according to the available literature [2]. A visual analysis of these scatterplots shows that both methods provide reasonable results, but in the first one for some fields and acquisitions, a strong overestimate is obtained, whereas the second one, accounting for the effect of the (moderate) vegetation, is in better agreement with the ground truth. This is confirmed quantitatively by evaluating the correlation coefficient and RMSE for the two cases. In fact, for the standard co-pol/correlation method, the correlation coefficient is rather low, i.e., 0.34, and RMSE is 13.6 vol%, whereas for the modified co-pol/correlation method the correlation coefficient is 0.64, and RMSE is 9.1 vol%. In addition, these results are practically coincident with those obtained in [2] for the same considered time period, thus confirming the correct implementation of the methods in the software tool.

Conclusions
In this paper, a software tool for the automatic retrieval of soil surface parameters from SAR data and via PTSM and PTSTCM has been presented. The available methodologies for soil surface parameters retrieval have been presented, and the basic theory behind PTSM/PTSTCM has been briefly but thoroughly reviewed.
The presented software tool allows to create maps of dielectric constant ε, large-scale roughness σ, and volumetric soil moisture m V from fully polarimetric and/or dual polarimetric SAR images. In fact, PTSM can be also employed in the absence of cross-polarization (HV), if a bare soil is imaged. In this case, only co-pol ratios and HH-VV correlations are needed. Availability of fully polarimetric data allows obtaining realistic retrieved values also in moderately vegetated areas. The software tool implements nine different algorithms, to be selected according to the typology of the observed scenario, as specified in Section 3. Thanks to its user-friendly interface, the tool is suitable for both expert and non-expert users. With this respect, the presented tool bridges a gap in the available software tools for PolSAR processing: in fact, commercially and freely available software tools for the processing of SAR data that include routines for basic manipulation of polarimetric SAR data (e.g., coherency and covariance matrix calculation, Pauli decomposition, etc.) certainly exist, see Section 1, but they do not implement easy-to-use methods for surface parameter retrieval. This has been the main motivation of our work.
In order to provide a first guide to the use of the tool, its application to a specific example from Germany's Demmin test-site, for which ground truth is available, has been presented. Finally, a comparison of retrieved soil moisture with in-situ measurements has been illustrated, showing a good agreement, at least for the method accounting for the presence of vegetation (modified co-pol/correlation). Obtained retrieval results are in agreement with what was expected according to the available literature [2].
The presented software is certainly amenable of further developments. For instance, use of compact polarimetric SAR data for soil moisture retrieval has been proposed [31]. At the moment, compact-pol SAR data cannot be employed by the software tool, since the PTSM and PTSTCM should be properly modified to cope with such data: this is left to future work.