**Applications of GIS-Based Software to Improve the Sustainability of a Forwarding Operation in Central Italy**

**Rodolfo Picchio 1, Francesco Latterini 2,\*, Piotr S. Mederski 3, Damiano Tocci 1, Rachele Venanzi 1, Walter Stefanoni <sup>2</sup> and Luigi Pari <sup>2</sup>**


Received: 10 June 2020; Accepted: 14 July 2020; Published: 16 July 2020

**Abstract:** Reducing potential soil damage due to the passing of forest machinery is a key issue in sustainable forest management. Limiting soil compaction has a significant positive impact on forest soil. With this in mind, the aim of this work was the application of precision forestry tools, namely the Global Navigation Satellite System (GNSS) and Geographic Information System (GIS), to improve forwarding operations in hilly areas, thereby reducing the soil surface impacted. Three different forest study areas located on the slopes of Mount Amiata (Tuscany, Italy) were analyzed. Extraction operations were carried out using a John Deere 1410D forwarder. The study was conducted in chestnut (*Castanea sativa* Mill.) coppice, and two coniferous stands: black pine (*Pinus nigra* Arn.) and Monterey pine (*Pinus radiata* D. Don). The first stage of this work consisted of field surveys collecting data concerning new strip roads prepared by the forwarder operator to extract all the wood material from the forest areas. These new strip roads were detected using a GNSS system: specifically, a Trimble Juno Sb handheld data collector. The accumulated field data were recorded in GIS Software Quantum GIS 2.18, allowing the creation of strip road shapefiles followed by a calculation of the soil surface impacted during the extraction operation. In the second phase, various GIS tools were used to define a preliminary strip road network, developed to minimize impact on the surface, and, therefore, environmental disturbance. The results obtained showed the efficiency of precision forestry tools to improve forwarding operations. This electronic component, integrated with the on-board GNSS and GIS systems of the forwarder, could assure that the machine only followed the previously-planned strip roads, leading to a considerable reduction of the soil compaction and topsoil disturbances. The use of such tool can also minimize the risks of accidents in hilly areas operations, thus allowing more sustainable forest operations under all the three pillars of sustainability (economy, environment and society).

**Keywords:** GIS; GNSS; forwarder; precision forestry; sustainable forest operations

## **1. Introduction**

To fulfil the relevant ecological, economic and social functions, sustainable forest management (SFM) [1] should include effective [2,3] and environmentally-acceptable forest operations [4]. Considering the above mentioned functions, SFM should minimize the negative impact of harvesting

on the environment without limiting the productivity while assuring forest workers' safety [4–7]. Modern machines, such as harvesters and forwarders (with wide, rubber tyres), have become common in forest utilization [8], also because they reduce the environmental impact in comparison with others utilization systems characterized by lower mechanization level [9]. The application of precision forest harvesting (PFH) may contribute significantly to the enhancement of efficient cut-to-length technology, i.e., a harvesting system in which trees are delimbed and bucked into assortments prior to subsequent transport to the landing site [10], and optimize SFM. PFH may be implemented by using interdisciplinary concepts, integrating the use of new technologies to create innovative solutions for efficient forest operations [11]. With particular reference to forwarding, the integration of Global Navigation Satellite System (GNSS) technology, Geographic Information System (GIS) and the on-board computing (OBC) hardware and software of modern forwarders, as well as advanced Information and Communications Technology (ICT), could enhance the future development of forest utilization [12–20].

Electronic devices integrated into modern forest machines used for forest operations do not only guarantee higher work productivity, but they could also reduce the environmental impact and enhance the safety of the workers [21–25].

Nowadays, the integration of electronic solutions with forestry practice (which can be seen as part of precision forestry) can contribute significantly to SFM and this creates a new best practice. It is possible for electronic devices to be implemented in all phases of the forest value chain, from the intervention planning to the product traceability. GIS technology could be used to analyze the topographic, ecological and morphological characteristics of the study area. GIS can help to design strip road network for timber harvesting and alternative extraction systems, with particular attention to economic aspects, minimizing negative impact on environment and providing a guarantee safety for operators [26–32]. GIS developed files can be implemented on the modern forwarders' information and communication technology (ICT) system; therefore, the designed strip road pattern can be displayed on the on-board screen and, thanks to the GNSS device, the operator can follow this strip road network, thereby limiting soil compaction [33]. Moreover, geo-data from the GNSS transformed in GIS could be integrated with work productivity and recorded using the standard for forest machine data and communication (StanForD) to carry out an economic evaluation of the entire study area [33,34]. In addition, a radio frequency identification (RFID) system allows for the identification of trees and marking them individually [35]. This technology showed good performances and moreover there are many possibilities of implementation [36].

Hence, aims of the present study were: (1) to apply GNSS and GIS technologies for the design of strip roads for forwarding operations in central Italy, and (2) to compare the net of the electronically-designed strip roads with those established in the forest by a forwarder operator according to his experience. Thus, it aimed to evaluate the effectiveness of precision forestry technology in the improvement of forwarding.

#### **2. Materials and Methods**

#### *2.1. Study Areas*

Three study areas were located on the slopes of Mount Amiata (Tuscany, Italy) in the Piancastagnaio district, in the Province of Siena.

Average annual rainfall in the study areas is approximately 1400 mm/yr−<sup>1</sup> with an average annual temperature of 11.5 ◦C. The stands were dominated by *Pinus radiata* D. Don, *Pinus nigra* Arn. and *Castanea sativa* Mill. in which three different forest treatments were conducted: clear-cutting (CC), thinning (TH) and coppicing with standards (CS), respectively (Table 1). In all three study areas, felling and processing were performed with a John Deere 1070G harvester, whilst extraction was carried out using a John Deere 1410D forwarder.


**Table 1.** Main characteristics of three analyzed study areas.

## 2.1.1. Clear Cutting (CC)

The total study area amounted to 0.31 ha dominated by *Pinus radiata* D.Don. (Figure 1), with a standing volume of 400 Mg·ha−1; all the timber was harvested as clear-cut (Table 1). The stand was located in a hilly area, the prevalent elevation was 618 m a.s.l. with a maximum elevation of 622 m a.s.l. and a minimum of 611 m a.s.l., with a predominately north-westerly aspect. The degree of the slope was between 5% (I Class) and 22% (II Class) (Figure 1).

**Figure 1.** Land cover and slope map of clear cutting (CC) and thinning (TH) study areas. CRS: WGS84-UTM32T. EPSG 32632.

## 2.1.2. Thinning (TH)

The total area of the forest sub-compartment equaled 2.30 ha, from which 1.10 ha was dominated by *Pinus nigra* Arn. and 1.20 ha of high stands of *Castanea sativa* Mill., *Quercus cerris* L. and *Fraxinus ornus* L. derived from a natural regeneration after an artificial pine stand (Figure 1). The total standing mass was 365 Mg·ha<sup>−</sup>1, from which 148 Mg·ha−<sup>1</sup> were harvested as thinning (Table 1). Only pine trees were cut, while all the broadleaved individuals were left upstanding. The area was also hilly, with a prevalent elevation of 626 m a.s.l. (maximum 653 m a.s.l. and minimum 613 m a.s.l.), mainly with a northwest aspect. The slope degree was between 5% (I Class) and 33% (II Class) (Figure 1).

## 2.1.3. Coppicing with Standards (CS)

The total study area was 1.93 ha, fully covered with *Castanea sativa* Mill coppice with standards (Figure 2) with a standing mass of 220 Mg·ha−1, from which coppicing of 200 Mg·ha−<sup>1</sup> had been conducted with the release of 50 standard trees per hectare. This area had the steepest slopes, with a prevalent elevation of 1,019 m a.s.l. (from 999 m a.s.l. to 1077 m a.s.l.), and a predominant northeasterly aspect. The degree of the slope ranged between 5% (I Class) and 45% (III Class) (Figure 2).

**Figure 2.** Land cover and slope map of coppice with standards (CS) study area. CRS: WGS84-UTM32T. EPSG 32632.

### *2.2. Field Reliefs*

A preliminary survey was conducted by recording the existing road network with a GNSS device. Then, in each of the three forest areas, six sample plots were randomly identified. Each sample plot had a circular shape and a radius of 10 m (surface of ca. 314 m2). Within the plots, all the strip roads open for harvesting (and forwarding) were detected using a GNSS device. Additionally, the length and average width of the strip roads were measured by using a measuring tape. Between the two systems used for distance measurement, there were no statistical differences. Finally, the coordinates of the center point of each sample area were recorded using the GNSS device to make it possible to transfer the sample area surfaces and locations on GIS. For data collection, a Trimble Juno Sb handheld was used. The Trimble Juno Sb was powered by the Windows Mobile 6.1 operating system and a 533 MHz Samsung S3C3443 processor. According to device's specifications, a real-time accuracy of 2 to 5 m was possible thanks to the integrated SBAS receiver. Subsequent post processing of the data using Trimble Delta Phase technology made it possible to reach a positional accuracy of 1–3 m [37]. However, for the aims of this study, post processing the data was unnecessary; therefore, the field data showed a positional accuracy of 2–5 m.

The data characterizing the relief were collected and recorded by GIS software: in particular, the open-source software Quantum GIS 2.18 Las Palmas, that allowed the creation of a line shape file of the strip road pattern, within the six sample areas for the three study areas. All the GIS files were geo-referenced in WGS84-UTM32T CRS (EPSG 32632). The data collected from the field surveys and elaborated using GIS technology were used to calculate three crucial parameters for the designed experiment: the surface impacted by forwarder passes, the length of the strip roads and strip road density.

## *2.3. GIS Implementation*

## 2.3.1. Preliminary GIS Steps

The GIS procedure developed for and applied in the design of a new improved network of strip roads needed two basic elements, i.e., a line shape file of the existing forest road network and a digital terrain model (DTM) of the area. The line shape file of the existing forest road network was derived from the GNSS survey. The DTM was built based on a topographic vector map of Tuscany, with a scale of 1:5,000. More precisely, a 2 m resolution DTM was derived with a QGIS plugin, using triangulated irregular network (TIN) interpolation. It should be noted here that the best DTM resolution freely available for the whole of Italy is currently 10 m [38,39]. Considering the size of the three study areas, it was decided in this case that it was inappropriate to use a 10 m DTM resolution; therefore, a DTM with a2mresolution was built using local geo-data.

## 2.3.2. New Strip Road Pattern Development and Determination of the Forwarder Passes Needed for Extraction

For the creation of the GIS-planned strip road network, the QGIS tool, Forest Road Designer (FRD), [40] was used. This is a GIS plugin that relies on a DTM and on points or a line shapefile reporting the zones. It generates another polyline meeting a series of design requirements established by the user (longitudinal slope and curvature radius among others) [40].

One of the most important parameters to be set using the FRD was the maximum slope gradient characterizing the new strip roads. This parameter had to be defined taking into account the characteristics of the forest machine that is supposed to be used in the strip roads. Considering the 1410D forwarder, a maximum slope of 45% was defined for stretches perpendicular to the contour lines and 25% for stretches parallel to them. In this way, machine-tipping risk was minimized. Moreover, with the FRD plugin, it was also possible to indicate some areas over which the newly designed strip roads should not pass (for example, high-value conservation areas). This was done by simply indicating such areas with a polygon shapefile. It was necessary in this case to indicate certain areas over which driving was forbidden as the study areas were surrounded by the properties of other owners. Another GIS procedure was developed and implemented to define the number of forwarder's passes needed for the extraction of all the timber from the three study areas, according to the GIS-planned strip road pattern. The first step was a calculation of the forwarding areas (i.e., forwarder accessible areas) where timber was within reach of the forwarder's boom. This was 12 m from the middle of each strip road, taking into account the fact that the working distance of the forwarder boom was 12 m. Once the forwarding areas were identified using the QGIS plugin fixed-distance buffer, it was possible to divide them using another QGIS tool: the polygon divider, which differentiated an input polygon layer (forwarding area) into a number of squarish polygons of a defined size. Knowing the forwarder loading capacity, which was 13 tons for the 1410D, and the harvested mass for each study area, it was possible to define the dimension of the sub-polygons or forwarding pixels (FPx), into which each forwarding area was divided. Each FPx surface corresponded to one forwarder load. The number of forwarder passes (NPs) corresponded to:

$$\text{NPs: } 2\text{-FPx - 2} \tag{1}$$

Finally, using the dedicated QGIS tool, the GIS-based strip road network was converted from an ESRI shape file format to .kml and .gpx ones, in order for them to be compatible with the forwarder's computer and to be visible on the screen.

#### 2.3.3. GIS Data Elaboration and Statistical Analysis

Having defined the GIS-based new strip road network, it was possible to calculate the soil disturbance parameters: the impacted surface, strip road length and strip road density for the GIS-planned study areas. These were calculated using the geo-fences of six sample plots from each study area. Following this, the data obtained were analyzed with Statistica 7.0 software: after checking for normality and homoscedasticity, both a one-way ANOVA and HSD Tukey test were used to find out if there were statistically significant differences among data collected manually in the forest and data obtained from the GIS.

## **3. Results**

## *3.1. Clear Cutting (CC)*

In the clear-cutting study area, the mean strip road length in reality was 31.93 m, corresponding to 1017 m· ha−<sup>1</sup> of strip road density and 35% of impacted surface. With GIS planning, the mean strip road length obtained from the sample plots was only 15.50 m, with a strip road density of 494 m·ha−<sup>1</sup> and 17% of surface impacted.

With regard to the harvested volume in the CC study area, the surface of each FPx was 325 m2, GIS analysis returned 9 FPx corresponding to 16 NPs. The forest road (fuchsia line, 3D model developed using QGIS tool QGIS2threejs, Figure 3b) is located along the southeast side of the study area and connected to a secondary truck road (light blue line) located to the northeast (Figure 3a).

**Figure 3.** (**a**) Actual forest road network in clear-cutting (CC) study area; (**b**) CC on 3D-model built using QGIS plugin QGIS2threejsbased on digital terrain model (DTM) and orthophoto map from 2013; (**c**) CC GIS (Geographic Information System)-planned strip road network; (**d**) CC study area divided into small areas contributing to one load, forwarder pixel (FPx). CRS: WGS84-UTM32T. EPSG 32632.

The GIS-planned strip road network (red lines) was created in a fir-shape (Figure 3c) reaching basically all the FPx (orange rectangles) of the CC study area (Figure 3d). The GIS-planned strip road network started from the existing forest road with a central axis, from which various branches departed to extract all the timber from the whole sub-compartment.

## *3.2. Thinning (TH)*

In the thinning study area, the manually measured mean strip road length was 25.00 m, which corresponded to a strip road density of 796 m·ha−<sup>1</sup> and 28% of surface impacted. With GIS planning (Figure 4), lower values were obtained likewise: mean strip road length was 12.50 m, strip road density equaled 398 m·ha−<sup>1</sup> and the surface impacted was 14%.

**Figure 4.** (**a**) Actual forest road network in thinning (TH) study area; (**b**) TH on 3D-model built using QGIS plugin QGIS2threejs based on DTM and orthophoto map from 2013; (**c**) TH GIS-planned strip road network; (**d**) TH study area divided into small areas contributing to one load, forwarder pixel (FPx). CRS: WGS84-UTM32T. EPSG 32632.

Considering the harvested volume in the TH study area, the surface of each FP was 866 m2, the GIS analysis returned 14 FPx corresponding to 26 NPs. The forest road was attached to the northwest side of the study area (fuchsia line, Figure 4a) and it was linked to a secondary truck road (light blue line), located to the northeast (3D model developed using QGIS tool QGIS2threejs, Figure 4b). The GIS-planned strip road network (blue lines) also had a fir-shape (Figure 4c), which gave access to all the FPx (blue rectangles, Figure 4d). The GIS-planned strip road network started from the existing

forest road, but in contrast to the CC, it was not possible to only develop a single central axis because of the presence of broadleaf groups (green shading), which had to remain upstanding. As a consequence, the GIS-planned strip road network had a more developed dendritic pattern.

## *3.3. Coppice with Standards (CS)*

In the coppice study area (six sample plots), the manually measured mean strip road length was 36.50 m, the strip road density was 1,162 m·ha−<sup>1</sup> and 41% of the area was impacted by the forwarder's driving. It was possible to reach the CS study area from driving from two sides: from the top of the slope and from the valley. Therefore, two GIS-planned strip road networks were designed: herringbone (CS\_Herr) and high-low (CS\_HL). The CS\_Herr strip roads basically ran parallel to the contour lines, while the CS\_HL ran in a perpendicular direction towards the contour lines. The GIS-planned strip road mean length (from the sample plots) came to only 13.00 m and 10.83 m, for CS\_Herr and CS\_HL, respectively. The strip road density was 414 and 345 m·ha−<sup>1</sup> for each design, while the impacted surface amounted to only 14.49% and 12% of the area for Herr and HL, respectively. The harvested timber mass of 200 Mg·ha−<sup>1</sup> required 650 m<sup>2</sup> of FPx, which in the GIS analysis amounted to 30 FPx, corresponding to 58 NPs for CS\_Herr, and 27 FPx corresponding to 52 NPs for CS\_HL. The study area visible on the topographic map was surrounded by the existing forest road network: from the west and east sides (Figure 5a). To the west, there was the main truck road (green line), from which a forest road (fuchsia line) departed and ran near the western side of the sub-compartment. To the east, there was the presence of a secondary truck road (light blue line), from which another forest road (red dotted line) started running along the southern side of the study area. A 3D model of road locations was developed using the QGIS tool QGIS2threejs (Figure 5b). The CS\_Herr GIS-planned strip road network (blue lines) had a perpendicular layout in comparison with the CS\_HL one (red lines, Figure 5c,d). The FPx of CS\_Herr (purple rectangles) were slightly different to the CS\_HL FPx (orange rectangles, Figure 5e,f). The GIS-planned strip road network departed from the existing forest road, developing in a parallel manner to the contour lines in CS\_Herr and perpendicularly to those in CS\_HL.

## *3.4. Statistical Analysis*

There were statistically significant differences between the real and the GIS-planned strip road networks in the CC and CS study areas, for both the CS\_Herr and CS\_HL models (Table 2). At the same time, there were no statistically significant differences in the TH area (the ANOVA p-value was 0.07, close to the level of significance).

However, for all three study areas, GIS-planning brought a considerable reduction in the area of impacted soil, with a percentage reduction between 50% and 71%. In the CS study area, there were no significant differences between the two GIS-based strip road patterns (CS\_Herr and CS\_HL). Therefore, GIS-planning considerably decreased the area of impacted soil.

**Figure 5.** (**a**) Actual forest road network in coppice with standards (CS)study area; (**b**) CS on 3D-model built using QGIS plugin QGIS2threejs based on DTM and Orthophoto from 2013; (**c**) herringbone (CS\_Herr) GIS-planned strip road network. (**d**) high–low (CS\_HL) GIS-planned strip road network; (**e**) CS study area divided into small areas contributing to one load, forwarder pixel (FPx) according to CS\_Herr strip roads network; (**f**) CS study area divided into small areas contributing to one load, forwarder pixel (FPx) according to CS\_HL strip roads network CRS: WGS84-UTM32T. EPSG 32632.


 separately for each forest study area. "\*\*\*" in the first column indicates p-value

**Table 2.** Overall results of one-way ANOVA and HSD Tukey test conducted

## **4. Discussion**

As found and demonstrated in several other studies, reducing the area of impacted soil during forest utilization is a good indication of SFM standards [41–43]. In this study, the effectiveness of the advanced electronic systems in reducing soil impact has been demonstrated. Thanks to the application of GNSS and GIS precision forestry tools for the planning of strip road networks, there was a reduction of 50%–70% in the area impacted in comparison with the plots on which the strip roads were created during the harvesting operation.

A GIS planned strip road pattern can also be beneficial from a social point of view. For instance, it was helpful to plot strip roads on slopes with a limited gradient, improving safety and maneuverability. Thanks to technological progress, which in the last years has led to an efficient integration of electronic devices in modern forest machines, such as harvesters and forwarders, it is possible to take one step further and transfer GIS files onto these machines.

Integrated GNSS technologies and modern ICT systems can visualize an optimal strip road pattern on the machine's display and help the operator drive in a comfortable, safe and efficient way in the forest (Figure 6).

**Figure 6.** Integrated GNSS, GIS, OBC and ICT—precision forestry optimizing workflow and contributing to good SFM standards.

Such a possibility is still available in modern harvester and forwarder models, which have OBC with dedicated software, such as John Deere TimberMatic Maps or TimberOffice [44,45], Ponsse Opti2 [46] or Komatsu MaxiXT [47]. These OBC systems can record the data from the harvested or processed timber through the StanForD standard, thus also providing the operator information about the work productivity and quality [34]. Furthermore, integrating the positioning data from GNSS and GIS, with productivity data from the StanForD data, acquired by the harvester or forwarder OBC system, can be helpful for the forest inventory [48] or for building decision support systems [49]. Moreover, vibration and ultrasonic sensors applied to a forwarder's OBC can record data on vehicle stability [50] and rut depth [51]. Although the suitability of modern technologies to improve the sustainability of forest operations has already been highlighted by scientific research, very little has been conducted in the Mediterranean region. This study therefore aimed to be a starting point in central Italy, demonstrating the effectiveness of a GIS-GNSS approach in decreasing the negative impact of forwarding.

Considering the above, another important aspect to be underlined is the possibility of using these technologies in small-scale forestry, though with rather lower level of accuracy. A feasible example of this could be smartphone use for improving forest utilization [52]. Smartphones are able to act as low-cost GNSS receivers, also under forest canopy cover, with sufficient precision, i.e., about 9 m of accuracy, which should be sufficient for small-scale forestry use [52,53]. Many smartphone applications, developed both for Android and for iOS systems, are able to display geo-data, geopoints, geolines and geo-fences files in .kml or .gpx format, and locate the operator's position. However, even if ca. 9 m accuracy is not sufficient for a forestry-fitted farm tractor driving (following a GIS-based strip road pattern displayed on the smartphone's screen), there are other useful functions which may be available. It can be very helpful for forest workers, for example, to display the geo-fence of the treatment area on the smartphone screen, allowing them to remain within the land boundaries or to avoid restricted areas, such as biodiversity hotspots.

A further step ahead in the integration of navigation technologies on forest machines could be represented by the development of tele-operated or unmanned forest vehicles. To reach this goal, which has been achieved in agriculture [54], there is the need to integrate in forest machines differential GNSS (DGNSS) technology, such as radio-beacon differential GNSS (RBDGNSS, or real time kinematic (RTK) [55,56], inertial measurement unit (IMU) sensors [18,57] and simultaneous localization and mapping (SLAM) algorithms [19].

## **5. Conclusions**

In recent years, several improvements have been observed in forestry, mainly the growing interest in sustainability, due to the importance of forests as environmental and social value [58,59].

Consequently, one of the most important purposes of the scientific research on forest utilization is to minimize the negative impact of forest operations, specifically on soil disturbances. Cutting edge technology and electronic devices could be powerful instruments used to reach this goal. In fact, using technological innovations, which are often mistakenly considered as something negative for the environment, may turn very helpful in forest operations.

In the presented research, the study confirmed that GNSS and GIS were useful technologies for forest operations and could improve SFM. GNSS and GIS resulted in being very helpful, both for real strip roads (established in the forest) detection and for electronically-designed strip road network.

Additionally, the use of the GIS-planned strip roads showed that soil impact due to forwarding may be decreased by 50%–70%. The herein presented precision forestry approach may be considered an efficient strategy for improving forest operations and SFM. Obviously, the practical implementation of such an approach in real forest yards in central Italy requires further steps in the training of operators, but could be very helpful in improving the sustainability of forest operations. Therefore, the presented findings can be used to improve forest utilization, also with the application of advanced technology, such as GIS and GNSS, in order to reach the effective sustainability of the whole value chain.

**Author Contributions:** Conceptualization, R.P. and F.L.; methodology, R.P. and F.L.; software, F.L.; validation, R.P., P.S.M., R.V.; formal analysis, F.L., D.T., W.S.; data curation, F.L., D.T.; writing—original draft preparation, F.L. and R.V.; writing—review and editing, R.P., F.L., P.S.M., R.V., D.T., W.S. and L.P.; supervision, R.P. and L.P.; funding acquisition, M.P. and L.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study was performed within the framework of the Polish Ministry of Science and Higher Education's programme: "Regional Initiative Excellence" in the years 2019–2022, Project No. 005/RID/2018/19. The APC was funded by the project AGROENER ("Energia dall'agricoltura: innovazioni sostenibili per la bioeconomia"), founded by Italian Ministry of Agriculture and Forestry DD n. 26329/2016.

**Acknowledgments:** The authors gratefully acknowledge MIUR (the Italian Ministry for Education, University and Research) for support (Law 232/2016, Italian University Departments of Excellence).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
