SALBEC – A Python Library and GUI Application to Calculate the Diurnal Variation of the Soil Albedo

The diurnal change in θ_s is an astronomic property and is well described by mathematics. The relationship between α_s and θ_s is computed using the set of empirical equations proposed by Cierniewski et al. (2018a) and broadly discussed ibidem. Here, the basic principles are summarised to illustrate building the computer science implementation. In the first step towards obtaining the diurnal a variation at a specific place and time, we calculate the soil a solely based on the soil spectrum and roughness. Soils spectra can easily be obtained either from European Land Use and Cover Area frame Statistical Survey (LUCAS) (Stevens et al. 2013) or from the Global Soil Spectral Library (Rossel 2009, Rossel et al. 2016). The soil roughness parameters can be calculated based on the direct measurements or estimated based on agricultural treatments according to Table 3 available in the study by Cierniewski et al. (2018b): HSD ∼ 5 mm, 10 mm and 25 mm and T_3D 1.05, 1.1 and 1.25 for smooth harrow, disc harrow and plough, respectively (Cierniewski et al. 2018b). Reflectance spectra in the range of 350–2500 nm was used to describe the a of the studied soils when θ_s = 45°. The first two equations transform the soil parameters into coefficients of the model: (1) α45=0.33−0.1099T3D−5794.4d574−510d1087+7787.2d1355+12161d1656+6932.8d698 \matrix{ {{\alpha _{45}} = 0.33 - 0.1099{T_{3D}} - 5794.4{d_{574}} - 510{d_{1087}} + } \cr {7787.2{d_{1355}} + 12161{d_{1656}} + 6932.8{d_{698}}} \cr }

This equation computes the overall a level of the soils with a given roughness at θ_s = 45° (α₄₅), where T_3D is the roughness index, defined as the ratio of the real surface area within its basic area unit to its flat horizontal area (Taconet, Ciarletti 2007). The index, d, relates to the soil reflectance data transformed to its second derivative for a specified wavelength in nanometres. The indexes are specific for bare soils and were established experimentally on 153 samples of soils from several climatic zones. The second derivative is independent on the absolute values of reflectance and stabilises the results against the total amount of reflected energy. (2) sa=6.26×10−7+0.0043(HSD)−1.418 {s_a} = 6.26 \times {10^{ - 7}} + 0.0043{\left( {HSD} \right)^{ - 1.418}}

This equation calculates s_α as the slope of the a rise relative to θ_s < 75°. The HSD is another roughness index that describes the standard deviation of the relative height of a soil surface within its basic area unit (Taconet, Ciarletti 2007). (3) αθs=α45[ 1+sα(θs−45) ] {\alpha _{\theta s}} = {\alpha _{45}}\left[ {1 + {s_\alpha }\left( {{\theta _s} - 45} \right)} \right]

This equation calculates the a of the studied soil surfaces (α_θs) for θ_s < 75°. (4) αθs=exp(a+cθs1+bθs+d(θs)2) {\alpha _{{\theta _s}}} = {exp}\left( {{{a + c{\theta _s}} \over {1 + b{\theta _s} + d{{\left( {{\theta _s}} \right)}^2}}}} \right)

This equation enables the adjustment of the entire range of the relationship between a of these soils and θ_s, for θ_s ∈ [0°, 90°]. Symbols a, b, c and d denote the fitted parameters of the equation.

The process of converting the input parameters into a soil albedo model is divided into three steps (A–C), as shown in Figure 2. As observed by Cierniewski et al. (2018), the dependency between θ_s and a is almost linear if θ_s∈[0°, 75°] and can be approximated linearly by Eq. (3). The relationship between θ_s∈[75°, 90°] and a is more complex owing to its non-linearity and is approximated by a fitting process on the intermediate data using Eq. (4).

Fig. 2

First, we calculate the values of α₄₅ and s_α from all the input data, i.e. the soil reflectance spectra and soil roughness (T_3D and HSD) using Eqs (1) and (2). Based on this, a list of points (L_{θs, αθs}) is generated, where θ_s < 75° and α_θs are obtained from Eq. (3). The list of points is generated with a density of θ_s = 1°. Finally, the last point, L_{90, 1} = (90, 1), is appended to the list, where L_{90, 1} corresponds to θ_s = 90° and α₉₀ = 1 at the time of sunset (Monteith, Szeicz 1961, Wang et al. 2004, Oguntunde et al. 2006). Subsequently (Fig. 2B), we use the standard least-square fitting algorithm to fit Eq. (4) to L_{θs, αs}. Unfortunately, the relationship between θ_s and α_s in the critical range of θ_s∈[75°, 90°] is not approximated by any auxiliary parameter. Therefore, the result of fitting between 75° and 90° did not fully correspond to the measurements obtained by Cierniewski et al. (2015). Through experiments, we found that the reduction of b in Eq. (4) by 0.01 resolved this problem (Fig. 2C). Finally, we have a tuple of parameters (a, b, c, d) applied to Eq. (4) to create a model: Ã_{a, b, c, d}(θ_s), which is denoted as à henceforth. à facilitates the calculation of a under clear-sky conditions for any θ_s∈[0°, 90°].

à is valid for the full range of θ_s, but the roughness parameters have limitations (T_3D∈[1.001, 3.5] and HSD∈(0, 100]). These limits do not exceed the measured values of T_3D∈[1.01, 2.99] and HSD∈(0, 60] (Cierniewski et al. 2019).

Although the model à allows mapping the values of θ_s onto α_θs, the values of θ_s are not connected with the actual time of the day. The diurnal change in θ_s is an astronomical expression that can be determined for each location on the Earth on any day of the year. The minimum inputs are the geographic location of the sample expressed as the latitude and longitude (l, ϕ), and the day of a year (DOY). Most astronomical software work primarily in Coordinated Universal Time (UTC). All calculations were made in solar local time (SLT), where the noontime for each location is always at 12:00 (local), as recommended as a time regime for the a calculation (Oguntunde et al. 2006), making the results comparable for each point on the Earth. To recalculate UTC to SLT, a time shift, denoted as D_t, specific for each day and location, is required. This is the difference between the moment of the highest sun position, when θ_s = min (θ_s). For SLT, min (θ_s) is at 12:00 AM, so D_t = t_noonUTC − t_noonSLT.

Modelling the diurnal a variation using the time at a given location as a direct input is impossible. It requires the prior calculation of a hash table, L_{t, θ}, where t denotes the discrete points in SLT between sunrise and sunset, and θ are values mapped onto the specific t keys. The time interval (I_L) of L_{t, θ} defines the density of the key series and affects the balance between the computation time and the accuracy of the results. The result of the application of Eq. (4) on the θ_s component of L_{t, θ} provides a discrete list of α_s assigned to specific keys of L.

Estimation of T_o was the main motivation for creating SALBEC. Cierniewski and Jasiewicz (2020), following the concept of Sellers et al. (1995), had specified T_o for the moment of day when a approximates its mean diurnal value αs¯ \overline {{\alpha _s}} . The diurnal range of θ_s depends on the geographic position and DOY; thus, αs¯ \overline {{\alpha _s}} is calculated by transforming Eq. (4) into the integral. As L_{t, θ, α} is discrete, the following is used to calculate αs¯ \overline {{\alpha _s}} : (5) αs¯=1| L |ΣLαs \overline {{\alpha _s}} = {1 \over {\left| L \right|}}\Sigma {L_{{\alpha _s}}} where:

L_αs are the discrete values of α_s in the given day and location, and

Through experiments, we found that L intervals denoted as I_L = 1 s is a good trade-off between the computation time and prediction accuracy and allows approximation of T_o with an accuracy of 1 s in reasonable computational time. A known αs¯ \overline {{\alpha _s}} allows the estimation of T_o = α_s ± ɛ, where ɛ is the acceptable deviation from αs¯ \overline {{\alpha _s}} . Estimating the T_o requires the discrete set of L_{t, θ, α} for the last L_t, where a<αs¯ a < \overline {{\alpha _s}} . Accordingly, αs¯−ε \overline {{\alpha _s}} - \varepsilon and αs¯+ε \overline {{\alpha _s}} + \varepsilon are used to obtain an acceptable range of T_o. All steps of the calculation are shown in Figure 3.

Fig. 3

SALBEC (Fig. 4) is divided into two distinct parts that are distributed together. The first part is the Python library designed to be imported into the Python environment, and the second part is its graphic user interface (GUI) front-end, which can be started directly from the operating system. The GUI is designed to provide the SALBEC functionality to non-programmers but is separated from the core library. Both environments, the Python library and the GUI, have the same functionality, except some GUI-only features that support processing in the visual environment. The Python library can be easily integrated and extended into more complex data processing scripts, while the GUI has a workflow path limited by its functionality. The core library contains three modules (Fig. 4): Soil Database (soildatabase.py), Soil Albedo (soilalbedo.py) and Diurnal Albedo (diurnalalbedo.py). A working example of running a simple analysis using the Python script is included in the source code directory (example.py).

Fig. 4

The soil database module is a supporting module designed to store and manage soil spectra. Soil spectra can be imported as text or spreadsheet files and can be imported separately during the run-time. However, considering that individual imports for each calculation could be an error-prone task, a separate module was created to import and manage the data. The database is designed as a list of soils, where each entity is stored as a single serialised binary object. Each object contains the spectrum as a wavelength–reflectance pair, geo-location of the soil sample (if available) and pre-computed spectrum derivatives necessary for Eq. (1). Soil names are a convenient way to access spectral data with intelligible names. The module is designed as a single class that facilitates basic management tasks – adding, removing, modifying and renaming the individual entities – over the database. The soil database is independent of the other modules and can easily be expanded into a more universal, domain-agnostic data management tool. The software is distributed with the reflectance spectra for the 33 most extensive agricultural regions globally, averaged based on the georeferenced topsoil samples acquired from national, continental and global soil spectral libraries. These include the European LUCAS (Stevens et al. 2013) and the Global Soil Spectral Library (Rossel 2009, Rossel et al. 2016). The data set and the estimation of optimal observation time for the above-mentioned regions are extensively discussed in the study by Cierniewski and Jasiewicz (2020).

The soil surface albedo module is a step towards calculating the diurnal a at specific locations and times. The module soilalbdo facilitates the calculation of the a variation for the full range, [0°, 90°] of θ_s. The module contains two classes (soil and soilCurve classes) and a few auxiliary functions, mostly to conduct the import and export of soil spectra from and to the soil database.

Class soil is designed to work directly with the spectra. It facilitates reading the spectra directly from the text file or reads pre-processed data from the soil database. An individual spectral curve can be imported from the wavelength–reflectance tabular data or read from the soil database. For newly imported data, the class extracts the local second-derivative values of the wavelength/reflectance curve. The second derivative is calculated from the empirical values, initially interpolated to the standardised resolution of 1 nm and subsequently fitted to the curve at around 10 nm. The calculations are based on the SciPy built-in functions (interp1d and derivative). The class also provides an auxiliary method for drawing a spectrum in a publishable format and methods to communicate (read and write data) with the soil database.

Class soilCurve is designed to calculate the clear-sky a of the soil. It fits the free parameters of Eq. (4) to the full range of θ_s. It accepts input values provided by the soil class and the two roughness coefficients (T_3D and HSD) provided by the user. Other class capabilities include the auxiliary functions to visualise the curve (Fig. 5A), along with the export parameters and the albedo model as a spreadsheet file.

The third module, diurnal albedo, contains one main class albedo and a batch of auxiliary functions. The class is the central part of the computational routines of the software and accepts l, ϕ in decimal degrees along with the DOY. These parameters are used to calculate the hash table L(t, θ, a) for a given location and day. Internally, class albedo uses the astral package to calculate specific points during the day and heights of the sun in terms of the solar height angle (θ_h), and it internally recalculates θ_h = 90°−θ_s. The class also calculates αs¯ \overline {{\alpha _s}} and T_o.

The complete results of the albedo class include a pandas data record. The record contains all critical moments of the day, including T_o ± ɛ. The class also provides methods to display the results in a graphical format (Fig. 5B) and to export them as spreadsheets. External batch functions facilitate data processing in a batch mode, for a specific period, or even an entire year. A detailed description of the class functionality is described in the software documentation and source code.

The GUI is designed in the ‘ribbon’ style to provide quick and intuitive access to all the parameters accepted by the Python modules described above. The application is just one window (Fig. 6) divided into the results (Fig. 6A) and the input ‘ribbon’ bars (Fig. 6B1–B3). The input bar is visually divided into three separate parts.

Fig. 5

Soil surface (B1) is where the user can choose the soil spectrum and set the soil roughness parameters (‘T3D’ and ‘HSD’ fields). The soil spectrum as an input to the model can only be selected from the list of soils represented by their names. The roughness parameters, T_3D and HSD, must be set manually, but the GUI prevents entering values outside of the acceptable range. This part of the software facilitates the investigation of the relationship of θ_s ∼ α_s as a line plot (Fig. 5A). The model is always re-fitted whenever a user decides to change the input parameters.

Geo-location (B2) is where the user can set the coordinates (latitude and longitude expressed as decimal degrees) of the expected location. If the soil spectrum has a geo-location, it is automatically set when a soil is selected from the list, but users can adjust the location manually, either by setting values in the input boxes or by clicking on the world map. Only valid geographical coordinates are accepted.

Days of the analysis (B3) is intended to set the date range. The day of the analysis can be set as a single day (current day by default), a range of dates or the entire year. If a range of dates or the entire year is selected, the number of days to be analysed can be reduced by increasing the daily interval (in days). The error field allows acceptable deviations from αs¯ \overline {{\alpha _s}} to be included in the results. The error values must be set as a list of percentages separated by a comma.

Pop-up dialogue windows are reduced to the minimum, except for the Soil Manager, which is a separate, GUI-only sub-module. This separate sub-module is divided into two tabs – the collection manager and the spectrum manager (Fig. 6C1 and C2). These do not cover some functions of the underlying core library layer and are designated to manage the soil spectra visually. The additional benefit of the soil manager is the soil spectra browser and visual exporter. Each spectrum must receive a human-readable name during the import to the soil database. The collection manager offers an easy tool to limit the number of named spectra displayed in the soil surface.

Fig. 6

The software presented in this study allows the modelling of the a variation of bare soil surfaces with differing roughness. The modelling is simplified to a clear-sky condition and refers to the optimal setup when soil surfaces can be viewed from a satellite. Eqs (1)–(4) used in this study are models fitted to the measurements acquired during several field campaigns in three countries (Poland, France and Israel; see Cierniewski et al. 2015, 2018b, 2019 for details). To demonstrate how the software refers to the real situation, we used the data acquired in Israel during the field campaign in 2015 when near clear-sky conditions were observed.

These data were collected at two locations: N30°59′16″, E34°42′15″ (Haplic Xerosol) and N30°59′16″, E34°42′15″ (Calcic Xerosol) on July 5th and July 7th, 2015, respectively. Each location was tested for the roughness formed by a smoothing harrow (smooth), disc harrow (rough) and plough (very rough), and its surface reached the air-dried state. Data to characterise the roughness of the tested soil surfaces by the HSD and T_3D indices were obtained using stereo photographs from a digital camera. The diurnal albedo variation of these surfaces with a size of at least 30 × 30 m was measured in the range of 35–2800 nm by an albedometer, LP PYRA 06 (by Delta OHM), placed ∼1.5 m above them. The albedo measurements started soon before 12:00 SLT and were continued with a 1-minute interval until sunset. In the laboratory, a reflectance spectrum was additionally obtained for soil material collected from each surface. A Vis/NIR FieldSpec 4 spectrometer with the Hig-Brite Muglight receptor working in the range of 335–2500 nm at 1 nm intervals was used for this purpose.

Figure 7 shows that there are differences between the real measurements and the results obtained by their adequate models. The shape of the measured and modelled curves is similar, implying that the differences are systematic. Models concerning the measured data are either under- or overestimated. Such results suggest that the error does not stem from the uncertainty of measurements. As the impact of atmospheric conditions was negligible, we presume that the mentioned differences result from the imprecision of the two roughness measures – T_3D, which impacts α₄₅ and the height of the curve over the x-axis, and HSD, which affects the curve slope for θ_s∈[0°, 70°]. The field of view of the albedometer is larger than the area for which the roughness was estimated. Moreover, the original data collected by Cierniewski et al. (2018) show the relationship between soil T_3D and measured α₄₅. Similarly, the relationship between the HSD and the slope of the models fitted to the measured data is not ideal. Differences shown in Figure 8 do not exceed the variations of T_3D and HSD shown by Cierniewski et al. (2018) in Figures 6 and 8.

Fig. 7

Fig. 8

At that moment, we could not compare our results with those of other similar tools, since such tools did not exist to our knowledge. The only known software pertains to urban areas (Chimklai et al. 2004) and not cultivated soils, but the results reported by its authors show that differences between the modelled and measured values are even higher than those reported in this study. Eq. (4), which is the core of the software described in our study, works with almost the same error with all types of soils and accepted ranges of surface roughness. Future observations and more data will facilitate the formation of new, potentially more detailed models tailored to specific conditions. The general approach proposed by Ziar et al. (2019) is very promising; however, at the present moment a comparison with our proposal would be disputable. Ziar's model is based on different assumptions and requires several, not always, available data, but our model is intended to work using solely data obtained from the laboratory or public databases.

Two factors potentially affect the computational performance of the SALBEC library. First is the procedure of model fitting accomplished by the soil albedo module that is performed once for each set of surface parameters to avoid redundant calculations. The second challenge is the calculation of the daily variation of a. Here, the computation time depends on the length of daylight for a given DOY and, thus, indirectly on the length of the L_{t, θ}. On a desktop computer using Intel Core i5-9400F CPU 2.90 GHz, the time of the fitting procedure varied between 0.14 s and 0.171 s. Diurnal variation of a for a daylight of 12 h varied between 0.0101 s and 0.0291 s for 30 repetitions. Calculations of the entire year take <15 s inside the GUI interface, including the time taken to display a widget filled with calculation results.