Erosion Map Reliability Using a Geographic Information System (GIS) and Erosion Potential Method (EPM): A Comparison of Mapping Methods, BELGRADE Peri-Urban Area, Serbia

: Soil erosion is a product of natural and anthropogenic factors and, at the same time, an economic and environmental concern. One of the methods applied to calculate the intensity of erosion is the erosion potential method (EPM), with two possible procedures for determining the average erosion coefﬁcient of an area: analytical and graphical. Using GIS and EPM methods, without ﬁeld observations of erosion, based on cartographic materials and satellite images, erosion maps were created for 1970 and 2018, for part of the peri-urban area of Belgrade. Based on the created erosion maps, the values of the mean coefﬁcients of erosion, as well as the arithmetic means for the study area, were determined for the settlements. The aim of the study is to assess the reliability of the mean coefﬁcient of soil erosion, obtained from the erosion map created from the cartographic materials and satellite images, without ﬁeld observations of erosion. Thus, the obtained values of the mean erosion coefﬁcient were compared with the values obtained from the erosion map with ﬁeld observation and the values obtained by the analytical procedure. Statistical analysis ( F test) for 1970 and for 2018 determined a high degree of reliability ( p < 0.05) of the mean erosion coefﬁcients of the area obtained from erosion maps that were created from cartographic materials and satellite images without ﬁeld observation. Regardless of the procedure for determining the mean erosion coefﬁcient, a signiﬁcant decrease in soil losses was observed, from 10.64 to 5.97 t ha − 1 year − 1 (average annual speciﬁc production of sediments, year 1970 and 2018, respectively).


Introduction
Soil erosion is the oldest permanent phenomenon of the collapse of the upper layer of the rock mass of the planet Earth, caused by natural and anthropogenic forces. That is a major problem around the world because of its effects on soil productivity, nutrient loss, siltation in water bodies, and degradation of water quality [1]. Furthermore, if erosion is aided by sudden earthquake movements, the period of landscape change is significantly shortened [2]. On the other hand, the deposition of sediments on the river affects reservoirs and dams, increases their costs of maintenance, and, in the long run, makes them unusable [3]. By increasing urban surfaces, permeable surfaces are decreased while excessive forest exploitation and inappropriate soil cultivation disturb the natural soil structure, making favourable conditions for the intensification of soil erosion, rapid surface runoff on slopes, and maximal discharges in river beds [4]. Understanding how these processes occur and what areas are vulnerable to erosion is paramount to land management. Soil erosion models aid land management by helping us understand sediment transport and Choice of procedure: When should the analytical procedure and when should the graphical one be performed, to determine the average erosion coefficient of an area or basin? 2.
Accuracy of erosion coefficient data: Is it to possible to determine erosion coefficients and create an erosion map from cartographic materials, satellite images, or aerial photographs using a GIS, without surveying? This research is focused on the advantages, disadvantages, and quality assessment mapping of soil erosion applying a GIS, without surveying. The subject of this study is the development of a soil erosion map using GIS technology according to the EPM without a field survey of erosion based solely on cartographic materials and satellite images for a part of the peri-urban area in Belgrade. The aim of this study is to assess the reliability of the mean coefficient of soil erosion in relation to the application of EPM erosion mapping procedures by statistical analysis.
This method would allow for a more convenient and standardized process of creating an erosion map; the majority of the work would be done by analysing cartographic materials and satellite images, whereas a small percentage would be done in the field. Field work would be reduced to verifying the erosion map obtained from the cartographic materials and satellite images. This map would have more accurate boundaries of the erosion areas, as the boundaries of the vegetation areas from land use, measured in ArcGIS, would be used. An erosion map properly prepared by EPM and the calculated erosion coefficients represents the basis for further calculations in which climatic factors have a direct impact on the amount of eroded sediments.

Research Area
The area of this research extends between the north latitudes 44 • 36.8 and 44 • 44 and the east longitudes 20 • 16 and 20 • 36 . It includes suburbs located 20 km south of the centre of Belgrade. It consists of the settlements Rucka, Rušanj, Pinosava, Beli Potok, and Zuce, which belong to the hilly peri-urban area of Belgrade ( Figure 1). The total area of the research area is 5687.08 ha.
Most of the territory consists of different varieties of sandstone, limestone, and carbonate rocks. These areas are subject to erosion, with visible traces of gully erosion. A smaller part of the territory, about 15%, consists of serpentines.
The average annual precipitation for the area during the period 1960-1970 was 635 mm and the air temperature was 11.7 • C. The average annual precipitation for the area between 1960 and 2018 was 675 mm and the air temperature was 12.3 • C [37].
About 30% of the researched area is under deciduous forest vegetation and orchards, within which residential areas have been developed. The settlement of Beli Potok, located on the slopes of Avala (506 m above sea level), which is the highest point of the exploration area, has a slightly larger area under forest. The lowest point is 85 m above sea level in the village of Rucka. This area is characterised by minor areas of arable land and an increased area of fallow (abandoned) land. The valley sides gravitate towards the basins of two rivers: the Topčiderska and Bolečka. Water flows directly into the Sava River from the area around the settlement of Rucka. The southern part of the research area is furrowed with a 2-5 m deep gully. This area is characterised by a continuing increase in the number of inhabitants, which is accompanied by aggressive urbanization. The number of inhabitants increased in the period between 1971 and 2011 by 27.1%. The population belongs to the category of an old population, and the active agricultural population is constantly declining.  [39], relief map of Belgrade (b) [40], and study area boundary (5 settlements) (c).
About 30% of the researched area is under deciduous forest vegetation and orchards, within which residential areas have been developed. The settlement of Beli Potok, located on the slopes of Avala (506 m above sea level), which is the highest point of the exploration area, has a slightly larger area under forest. The lowest point is 85 m above sea level in the village of Rucka. This area is characterised by minor areas of arable land and an increased area of fallow (abandoned) land. The valley sides gravitate towards the basins of two rivers: the Topčiderska and Bolečka. Water flows directly into the Sava River from the area around the settlement of Rucka. The southern part of the research area is furrowed with a 2-5 m deep gully. This area is characterised by a continuing increase in the number of inhabitants, which is accompanied by aggressive urbanization. The number of inhabitants increased in the period between 1971 and 2011 by 27.1%. The population belongs to the category of an old population, and the active agricultural population is constantly declining.

Method
To determine the average erosion coefficient, the graphical and analytical procedure of the EPM were applied in this study. Erosion maps were made by the graphical procedure without a field survey for 1970 and 2018, and based on these maps, the average erosion coefficients were calculated. For the same respective years, the average erosion coefficients were calculated according to the analytical procedure as well as the graphical procedure using the erosion maps and field survey ( 2018. Statistical analysis evaluated the reliability of the obtained erosion coefficients from the erosion maps without the field survey.

Erosion Potential Method-EPM
The erosion potential method (EPM) was developed from the Method for the Quantitative Classification of Erosion (MQCE). [10] The original method was adapted by Prof. Slobodan Gavrilović based on the field research of 17 experimental field stations set up in several regions of Yugoslavia. The regions differed in several parameters (climate, geology, soil, slopes, and relief) and visible traces of erosion. The EPM does not investigate the physics of erosion processes and, as such, is suitable for areas where there is only partial information available and no erosion research has been performed. The method is "the most quantitative of all semi-quantitative methods" [41]. The EPM can provide data on total sediment production, erosion intensity, and suspended as well as transported sediments depending on the following factors: climatological (average annual temperatures and amounts of precipitation), vegetation, land use, types of soils, relief (slope and highest and lowest points of the basin), visible traces of erosion, volume of the basin, and catchment area. The aim of this research is to make an erosion map and determine the mean erosion coefficient, which is a synthetic factor and the most influential for determining the intensity of erosion.
The annual erosion intensity according to the EPM model is estimated according to the following equation: where: W year is the total annual erosion (m 3 /year/km 2 ). T is the temperature coefficient.
H year is the average yearly precipitation (mm). F is the basin area (km 2 ) and Z is the erosion coefficient.
The temperature coefficient (T) is calculated from the following equation: According to the EPM, the erosion coefficient can be determined by analytical and graphical procedures, and its value ranges from 0.01 to 1.5 for natural basins.
The analytical procedure is based on the knowledge of soil resistance, the protection of soil from the effects of the atmosphere, the visible traces of erosion, and the average inclination of the area. The analytical procedure is a numerical model of erosion that determines the average erosion coefficient of the area according to the following equation:

EPM Applied to the Study Area
The average erosion coefficients were determined by the analytical procedure for all five settlements for 1970 and 2018. The coefficients for X, Y, and ϕ are in Tables S1-S3. Tables S1-S4, in the supplementary file, are designed on the basis of the field and laboratory research of the experts on torrential waters and erosion in Serbia.
The value of coefficient Y is the soil erosion resistance coefficient, determined from a pedological map (Table 1, map No. 3).Y is calculated based on the established percentages of the presence of each type of land formation, multiplied by the value for that type, given in Table S1. This is the partial value of Y. Summing those partial values, one obtains the Y coefficient for the area. The value of coefficient X is the numerical value of soil protection from atmospherics and erosion by natural conditions-vegetation and the like-and it is expressed through the way of land use, namely in the following ways:

•
For 1970, X was determined based on a land-use map ( Figure 2, left), which was made using a combination of maps (Table 1, maps No. 4 and 6). They were combined in order to overcome the lack of data in terms of maps, so as to generate a more accurate land-use map. Two maps were combined to avoid the disadvantages of one in relation to the other; • For 2018, the value of coefficient X was determined on the basis of a land-use map, which was made from satellite images ( CORINE maps were not used because they have a high degree of generalization of spatial data [42] through five categories, from which arable areas as very important factors of erosivity are not clearly observed. For 1970, there is no possibility of its use. The mapping scale is 1:100,000 with an accuracy of not less than 100 m. . Example of land-use maps for 1970 (left) and 2018 (right) for the settlement of Zuce. Slope maps ( Figure 3) were made separately for the five locations, based on the map ( Table 1, map No. 7). In Serbia, zoning is most often used in the following ranges: 0° to 3°, 3° to 5°, 5° to 10°, and more than 10°.
According to the map ( The equation labels are the following: I = interval line spacing between contour lines on the map, on the basis of which the slopes of the terrain sides are defined; E = vertical distance between contour lines on a map; ctgα = cotangent of slope angles α = 1°, 2°, 3°,....., 45°. X is calculated based on the established percentages of the presence of each land-use category multiplied by the land-use value for that type, given in Table S2. This is the partial value of X. Summing those partial values, one obtains the X coefficient for the area.
The coefficient ϕ is the numerical equivalent of the visible and clearly expressed erosion processes in the research area according to the following:

•
For the year 1970, coefficient ϕ was determined on the basis of cartographic signs on maps ( ϕ is calculated based on the percent of the representation of visible erosion traces multiplied by the value in Table S3. This is the partial value of ϕ. Summing those partial values, one obtains the ϕ coefficient for the area.
The values of coefficients X, Y, and ϕ are presented in Tables S1-S3, respectively. The erosion intensities, erosion categories, and ranges of Z are presented in Table S4.
The value of coefficient J sr is the mean average slope of the area, determined by settlements based on the lengths of contours from the contents of the map ( where: Land 2022, 11, 1096 7 of 18 F = area of the settlement or territory; h 0 = difference between the lowest point and the first higher contour; L 1 = length of the contour above the lowest point; L 2 -L n−1 = length of contours; h = contour interval; h n = difference between the highest point and the first lower contour; L n = length of the contour below the highest point.
The graphical procedure for determining the average erosion coefficient for an area, via the method of erosion potential, is obtained in existing maps by measuring the areas by erosion type and calculating the percent of their share in the total area. The percentage value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient. value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient.

V2
Groups of buildings or settlements on slopes up to 10° if the buildings are not surrounded by an orchard or a forest, after which it is then mapped as a forest or orchard

Cartographic Method and Realization of an Erosion Map Using GIS Technology
The cartographic method is very common in the two procedures of the EPM model (analytical and graphical) in order to determine the average coefficient of erosion through the following procedures: planimetry, calculating a linear object's length, and mapping. To determine the average erosion coefficient according to the analytical procedure, the following surfaces were measured by planimetry on a georeferenced map in ArcGIS: the Forests, orchards, vineyards, pastures, and meadows on eroded areas with occasionally seen gullies, ravines, furrows, landslides, and scree, regardless of the slope of the terrain value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient. value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient.

V2
Groups of buildings or settlements on slopes up to 10° if the buildings are not surrounded by an orchard or a forest, after which it is then mapped as a forest or orchard

Cartographic Method and Realization of an Erosion Map Using GIS Technology
The cartographic method is very common in the two procedures of the EPM model (analytical and graphical) in order to determine the average coefficient of erosion through the following procedures: planimetry, calculating a linear object's length, and mapping. value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient. value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient. value of each erosion type is multiplied by the average erosion coefficient from Table 2, and thus one partial value of the erosion coefficient is obtained. The average erosion coefficient for the area is obtained by summing the partial values of the erosion coefficient.

V2
Arable land, in an erosion area with an average inclination of up to 3°, and plains with preserved vegetation (meadows and pastures)

V2
Groups of buildings or settlements on slopes up to 10° if the buildings are not surrounded by an orchard or a forest, after which it is then mapped as a forest or orchard

Cartographic Method and Realization of an Erosion Map Using GIS Technology
The cartographic method is very common in the two procedures of the EPM model (analytical and graphical) in order to determine the average coefficient of erosion through the following procedures: planimetry, calculating a linear object's length, and mapping. To determine the average erosion coefficient according to the analytical procedure, the following surfaces were measured by planimetry on a georeferenced map in ArcGIS: the vegetation surfaces, the categories of soil, the visible traces of erosion, and the lengths of contours.
To determine the mean erosion coefficient by the graphical procedure, the existing erosion maps (Table 1, maps No. 1 and 2) were georeferenced in ArcGIS. Considering that previous research for this area was performed in basins, it was necessary to include habitat boundaries and then calculate the areas under different types of erosion. To create a new erosion map, it was necessary to have a slope map, land cover map, and Table 2 (legend symbols for erosion types).
Land-use maps are made for each period for all five settlements as follows: For 1970, the land-use map was made using the maps (Table 1, maps No. 4 and 6). These two maps were combined in order to avoid the flaws of one in relation to the other.
For 2018, the land-use maps were made on the basis of satellite images ( CORINE maps were not used because they have a high degree of generalization of spatial data [42] through five categories, from which arable areas as very important factors of erosivity are not clearly observed. For 1970, there is no possibility of its use. The mapping scale is 1:100,000 with an accuracy of not less than 100 m.

V2
Groups of buildings or settlements on slopes up to 10° if the buildings are not surrounded by an orchard or a forest, after which it is then mapped as a forest or orchard

Cartographic Method and Realization of an Erosion Map Using GIS Technology
The cartographic method is very common in the two procedures of the EPM model (analytical and graphical) in order to determine the average coefficient of erosion through the following procedures: planimetry, calculating a linear object's length, and mapping. To determine the average erosion coefficient according to the analytical procedure, the following surfaces were measured by planimetry on a georeferenced map in ArcGIS: the vegetation surfaces, the categories of soil, the visible traces of erosion, and the lengths of contours.
To determine the mean erosion coefficient by the graphical procedure, the existing erosion maps (Table 1, maps No. 1 and 2) were georeferenced in ArcGIS. Considering that previous research for this area was performed in basins, it was necessary to include habitat boundaries and then calculate the areas under different types of erosion. To create a new erosion map, it was necessary to have a slope map, land cover map, and Table 2 (legend symbols for erosion types).
Land-use maps are made for each period for all five settlements as follows: For 1970, the land-use map was made using the maps (Table 1, maps No. 4 and 6). These two maps were combined in order to avoid the flaws of one in relation to the other.
For 2018, the land-use maps were made on the basis of satellite images ( CORINE maps were not used because they have a high degree of generalization of spatial data [42] through five categories, from which arable areas as very important factors of erosivity are not clearly observed. For 1970, there is no possibility of its use. The mapping scale is 1:100,000 with an accuracy of not less than 100 m.

Cartographic Method and Realization of an Erosion Map Using GIS Technology
The cartographic method is very common in the two procedures of the EPM model (analytical and graphical) in order to determine the average coefficient of erosion through the following procedures: planimetry, calculating a linear object's length, and mapping. To determine the average erosion coefficient according to the analytical procedure, the following surfaces were measured by planimetry on a georeferenced map in ArcGIS: the vegetation surfaces, the categories of soil, the visible traces of erosion, and the lengths of contours.
To determine the mean erosion coefficient by the graphical procedure, the existing erosion maps ( CORINE maps were not used because they have a high degree of generalization of spatial data [42] through five categories, from which arable areas as very important factors of erosivity are not clearly observed. For 1970, there is no possibility of its use. The mapping scale is 1:100,000 with an accuracy of not less than 100 m. Slope maps (Figure 3) were made separately for the five locations, based on the map ( Table 1, map No. 7). In Serbia, zoning is most often used in the following ranges: 0 • to 3 • , 3 • to 5 • , 5 • to 10 • , and more than 10 • .
According to the map (Table 1, map No. 7), with the use of a slope scale, slope areas were separated. If there was no slope scale, then a graph of the function was used: The equation labels are the following: I = interval line spacing between contour lines on the map, on the basis of which the slopes of the terrain sides are defined; In the process of creating erosion maps for 1970 and 2018, the ArcGis software package was used. By overlapping the map of land use and the slope area map with adjusting the transparency of the layer, the methods of land use on different slopes can be noticed. Figure 4 clearly shows the boundaries of arable land areas located on a slope over 10°. The layer "types of erosion" is then included. In this case, it is a "Strong 1" type of erosion according to Table 2, and it is mapped. If a map of transparent material is made, then hatching marks are used for mapping. If an erosion map for any historical period is made, as is the case for 1970, then a coloured sign is used in accordance with the available sources of information on visible traces of erosion ( The map of erosion for 1970 and the draft version of the erosion map on transparent material for 2018 were generated by overlapping the maps of land use and a slope map. The draft version of the erosion map on transparent material was created using hatching symbology and was a basis for the final version of the erosion map for 2018. In the process of creating erosion maps for 1970 and 2018, the ArcGis software package was used. By overlapping the map of land use and the slope area map with adjusting the transparency of the layer, the methods of land use on different slopes can be noticed. Figure 4 clearly shows the boundaries of arable land areas located on a slope over 10 • . The layer "types of erosion" is then included. In this case, it is a "Strong 1" type of erosion according to Table 2, and it is mapped. If a map of transparent material is made, then hatching marks are used for mapping. If an erosion map for any historical period is made, as is the case for 1970, then a coloured sign is used in accordance with the available sources of information on visible traces of erosion ( The map of erosion for 1970 and the draft version of the erosion map on transparent material for 2018 were generated by overlapping the maps of land use and a slope map. The draft version of the erosion map on transparent material was created using hatching symbology and was a basis for the final version of the erosion map for 2018. The draft version of the erosion map should be at the same scale as the topographic map, the land-use map, or the satellite images. Overlaying the land-use map or a satellite image with the erosion map on transparent material in the field, it could be noticed that the mismatch of agricultural or forest area borders and erosion area boundaries needed to be corrected. Based on visible erosion traces in the field, the erosion coefficient and erosion types were corrected on-site. In accordance with Table 2, corrections from the draft version of the erosion map were corrected on the final erosion map ( Figure 6).
The control of the mapped content is carried out on the basis of a plan on the observed discrepancies between the sources: images taken from the Geoserbia or Google Earth platforms and the content on a topographic map at a scale of 1:25,000. The draft version of the erosion map should be at the same scale as the topographic map, the land-use map, or the satellite images. Overlaying the land-use map or a satellite image with the erosion map on transparent material in the field, it could be noticed that the mismatch of agricultural or forest area borders and erosion area boundaries needed to be corrected. Based on visible erosion traces in the field, the erosion coefficient and erosion types were corrected on-site. In accordance with Table 2, corrections from the draft version of the erosion map were corrected on the final erosion map ( Figure 6).
The control of the mapped content is carried out on the basis of a plan on the observed discrepancies between the sources: images taken from the Geoserbia or Google Earth platforms and the content on a topographic map at a scale of 1:25,000. Table 2 was created as a modification of previously used legends of signs on erosion maps in the Republic of Serbia; the descriptions and ranges of the erosion coefficients by the types of erosion are shown in Table S4, based on the research in Serbia by Gavrilović [43], Lazarević [44], and Kostadinov [45].

Statistical Methods
Statistical methods were applied to determine the differences between the average erosion coefficients using arithmetic means in relation to the erosion recording methods and the recording time periods. The F-test was used to test the differences in the arithmetic means of the erosion coefficients of the methods used to determine the average erosion coefficient for 1970. The t-test was used to assess statistical differences between the relevant methods applied in one period and between two periods of erosion recording. Statistical analyses were performed in the IBM SPSS Statistics 26 software package [46].

Research Results
The results of the research are the erosion maps for 1970 ( Figure 5) and the erosion maps for 2018 ( Figure 6). For each year of the research, five erosion maps were made (one for each settlement).
For the 1970 map, there was no possibility of the field control of the mapped contents, so the maps (Table 1, Table 2 was created as a modification of previously used legends of signs on erosion maps in the Republic of Serbia; the descriptions and ranges of the erosion coefficients by the types of erosion are shown in Table S4, based on the research in Serbia by Gavrilović [43], Lazarević [44], and Kostadinov [45].

Statistical Methods
Statistical methods were applied to determine the differences between the average erosion coefficients using arithmetic means in relation to the erosion recording methods and the recording time periods. The F-test was used to test the differences in the arithmetic means of the erosion coefficients of the methods used to determine the average erosion coefficient for 1970. The t-test was used to assess statistical differences between the relevant methods applied in one period and between two periods of erosion recording. Statistical analyses were performed in the IBM SPSS Statistics 26 software package [46].

Research Results
The results of the research are the erosion maps for 1970 ( Figure 5) and the erosion maps for 2018 ( Figure 6). For each year of the research, five erosion maps were made (one for each settlement).
For the 1970 map, there was no possibility of the field control of the mapped contents, so the maps (Table 1, maps No. 3, 5 and 6) were used for correction, respecting the note below Table 2. For 2018, when mapping using GIS technology, two maps were produced, one on the transparent paper, to control the mapped content in the field, and another map after field observations. During the land inventory, no major deviations in the erosion areas or erosion coefficients were observed.  Table S4.
Based on the erosion maps produced by the maps and satellite images, Table 1, using GIS technology, the average erosion coefficients for each settlement were calculated for the periods 1970 and 2018 (Table 3); the arithmetic mean was calculated for the area afterwards.  Table S4.
By comparing the values of the arithmetic means between the two recording periods, a significant difference was noticed. The obtained values were tested by statistical analysis.  Table S4.
A land inventory was performed for all five settlements in the period between April and June 2018, after which the mapping was performed. The following is an example of the results of a land inventory of the settlement of Zuce: On the map (Table 1, map No. 7) for the settlement of Zuce, a vineyard is shown when it is not in the satellite image ( Table 1, No. 8) or in nature. In the southern part of the settlement of Zuce, there is a continuous forest on the Geoserbia portal, and hangars and a road network can be seen on the Google Earth website. The existence of special-purpose facilities was confirmed in the field. A meadow is shown in the northern part of the Zuce settlement, and a large complex of an "Ikea" shopping centre with a large parking facility was noticed at the spot. After the field inventory was completed, corrections were drawn and the final appearance of the erosion map was formed.
Based on the erosion maps produced by the maps and satellite images, Table 1, using GIS technology, the average erosion coefficients for each settlement were calculated for the periods 1970 and 2018 (Table 3); the arithmetic mean was calculated for the area afterwards. By comparing the values of the arithmetic means between the two recording periods, a significant difference was noticed. The obtained values were tested by statistical analysis.
For 1970, a total of four values of the average erosion coefficient of the area was determined: one by the analytical and three by the graphical procedure. The values of the average erosion coefficient were graphically determined on the basis of erosion maps ( Table 1, maps No. 1 and 2) as well as the erosion map without field observations. The analytical procedure determined the average erosion coefficient based on maps (Table 1,  For 2018, three values of the average erosion coefficient were determined: one by the analytical and two by the graphical procedure. The following were used for the analytical procedure: land-use maps (Figure 3, on the right) and maps (Table 1, maps No. 3 and  7). The value of the erosion coefficient on the draft version of the erosion map and the final erosion map was determined by the graphical procedure. Between these two maps, minimal deviations by settlement have been determined, and no significant difference between the arithmetic means was found. Table 4 shows the values of the arithmetic means of the average erosion coefficient range obtained by the analytical procedure and the graphical procedure for both periods. To control the results of mapping based on cartographic materials and satellite images using the GIS, the areas of erosion types were measured from existing erosion maps (Table 1, maps No. 1 and 2) and maps that are necessary for determining the average erosion coefficient via the analytical procedure. Table 4 shows the existence of differences between the arithmetic means. The observed differences were statistically significant at p < 0.05.

Statistical Evaluation of Results
The statistical analysis of the arithmetic means of the mean erosion coefficients determined by the analytical procedure for 1970 and 2018 (columns one and five of Table 4) shows that (M = 18.20, SD = 3.96); t (4) = 10.27, p = 0.001. By testing the arithmetic means of the average erosion coefficients determined from the maps obtained from the cartographic materials and satellite images for 1970 and 2018 (columns four and seven of Table 4), we discerned that (M = 17.40, SD = 4.56); t (4) = 8.53, p = 0.001. The results obtained from the average erosion coefficient between periods by any method are statistically significant and differ from each other with a probability of more than 95%.
An analysis of four arithmetic means (columns one, two, three, and four, Table 4) using the F-test, for the period 1970, revealed that the differences in the arithmetic means of the average coefficients, compared to the method of recording erosion during 1970, are not statistically significant because (M tot = 53.85, SD tot = 7.00); F (3, 16) = 2.39, p = 0.107. It was determined that each of the results obtained, regardless of the method, can be used in other research as well.
Using the t-test for the pairs of the arithmetic means of the average erosion coefficients, for the year 1970, the results were compared by column four with columns one and two ( Table 4). The test determined that maps created on the basis of cartographic materials and satellite images have acceptable results regarding the average erosion coefficient, since the differences in the arithmetic means are not statistically significant.
By comparative analysis, for the year 2018, it was determined that the arithmetic means of the average erosion coefficients for columns six and seven are "equal", and therefore the t-test between the two arithmetic means between columns five and seven was applied. It was found that there was no statistically significant difference between the arithmetic means of the average erosion coefficients because (M = −0.60, SD = 4.83); t (4) = −0.28, p = 0.795.
The results obtained within either 1970 or 2018, regardless of the method of determining the average erosion coefficient, have no statistically significant difference. This was determined with a probability of more than 95%.

Discussion
Prior to the development of a GIS, erosion mapping required extensive field work and highly trained personnel, who would assess the terrain slope, observe erosion patterns and the types of vegetation, and, based on knowledge of the geology of the land, determine the erosion coefficient according to Tables S1-S3. This way of working contained a lot of subjective assessments and required personnel who had the same criterion for evaluating erosive surfaces in the field.
The process of mapping erosion using cartographic materials, satellite images, and GIS technology is a specific graphical procedure. In this research, that procedure was used to create an erosion map for the year 1970, as well as for 2018. The precondition for making an erosion map for those years is to already have a land-use map and relief slope map for the area. If a relief slope map is not obtainable, it is necessary to have a map showing the relief by contours in order to produce a relief slope map. If a land-use map is not available for a certain period, it is necessary to have adequate topographic maps, such as maps and satellite images (  The erosion map for the current period should be constructed in two phases: First, the erosion map on the transparent paper is made, where the erosion area is presented by hatching symbology. Field corrections and erosion coefficients are entered into the map from field observations, and then in the second phase, using GIS, these field corrections are entered into the final erosion map. This procedure is called research center-fieldresearch center.
This research shows that the erosion coefficient obtained from the erosion map, made without visually verifying the site, is statistically acceptable. The accuracy of the presented contents is higher because the process is performed on the basis of the proportional and instrumentally determined boundaries of land use and clearly defined slope area.
The process of mapping, based on cartographic materials and satellite images, has a number of significant positive outcomes, including the following: the time required for field work is reduced, causing a subsequent reduction in erosion map production costs; preparations for field trips can be more comprehensive; projects are easier to implement by understanding erosion hazards and investment risks for projects; and proposals for anti-erosion measures are easier to define. Field checking is focused on critical and unclear erosion areas that are assessed using cartographic materials, so the situation is verified on the terrain and necessary anti-erosion measures are proposed.
Over 80% of the mapping of erosion areas can be completed based on cartographic materials and satellite images. The production of a draft version of a map can be performed by professionals who might not be experts in calculating erosion (technical personnel) as a preparatory document for the inventory of a study area. Field observations can be performed by students of forestry, geography, agriculture, construction engineering, etc., i.e., persons trained in identifying the visible traces of erosion. In observing the terrain, land use should be identified, such as arable land, forests, pastures, vineyards, etc., as should the visible traces of erosion; in accordance with the observations, erosion boundaries and coefficients should be corrected.
The draft version of an erosion map can be used in public discussions when proposing projects of anti-erosion measures and making a final erosion map of an area.
The weakness of the method lies in the dependency on reliable cartographic sources and aerial and satellite images. Making historical erosion maps using this method would depend on available sources for that period of time.

Conclusions
The research clearly indicated that, with the help of high-quality large-scale topographic maps and adequate aerial imaging, mapping can be successfully performed as preparation for field work or the production of broad-scale erosion maps for any historical period.
Creating a map based on cartographic materials and satellite images without field observation, it is possible to assess soil losses in previous periods and create a draft version of the map for the field observations of erosion. Based on obtained mean erosion coefficients and climatological factors, soil losses (Wsp) were calculated, amounting to 10.64 t ha −1 year −1 in 1970 and 5.79 t ha −1 year −1 in 2018.
A detailed analysis of anthropogenic factors on the intensity of erosion determined that in the period from 1970 to 2018 there was a decrease in arable land (by 59.7%) and livestock (by 42.4%). In the same period, there was an increase in the areas of forest (by 37.6%), meadows (by 219.3%), and infrastructure facilities (by 56.8%).
The results obtained using statistical analysis (arithmetic means) indicate that the value of the average erosion coefficient obtained by the cartographic materials and satellite images does not statistically differ from the results obtained by graphical and analytical procedures.
Statistical analysis (F-test) proved that the obtained value of the coefficient of erosion from cartographic materials and satellite images without going to the field is reliable, compared to the values obtained by other methods (Table 4)  The assumptions in this study have been tested experimentally by processing measurable cartographic and aerial photo documentation and should be understood as conceptual for the application of this research in other geographical areas, especially in areas with excessive erosion. We feel that it would be a good idea to compare the results with detailed field measurements, but there are no data for our study area at this time.