Deriving the Main Cultivation Direction from Open Remote Sensing Data to Determine the Support Practice Measure Contouring

: The P-factor for support practice of the Universal Soil Loss Equation (USLE) accounts for soil conservation measures and leads to a signiﬁcant reduction in the modelled soil loss. However, in the practical application, the P-factor is the most neglected factor overall due to high effort for determining or lack of input data. This study provides a new method for automatic derivation of the main cultivation direction from seed rows and tramlines on agricultural land parcels using the Fast Line Detector (FLD) of the Open Computer Vision (OpenCV) package and open remote sensing data from Google Earth™. Comparison of the cultivation direction with the mean aspect for each land parcel allows the determination of a site-speciﬁc P-factor for the soil conservation measure contouring. After calibration of the FLD parameters, the success rate in a ﬁrst application in the low mountain range Fischbach catchment, Germany, was 77.7% for 278 agricultural land parcels. The main reasons for unsuccessful detection were problems with headland detection, existing soil erosion, and widely varying albedo within the plots as well as individual outliers. The use of a corrected mask and enhanced parameterization offers promising improvements for a higher success rate of the FLD.


Introduction
Soil erosion caused by wind or water is a well-known global problem, but a natural process. The term soil erosion includes the processes of detachment, transportation, and sedimentation of soil particles [1]. A massive increase in eroded soil material and land degradation often occurs due to excessive and improper use by humans on agricultural land [2]. The intensity of the removed mass depends on erosive factors such as precipitation and the resulting surface runoff as well as erodible factors that take into account soil properties and land cover as well as practice management [3]. High erosion rates can lead to significant damage, initially affecting the farmland itself, known as on-site damage. The associated loss of the most fertile topsoil causes a decrease in soil quality and productivity and an impairment of important soil functions, which in turn leads to a loss of yield [4]. Furthermore, material deposits and siltation occur downslope within the field. If the detached soil material passes the edge of the field, not only neighboring fields but also infrastructure elements such as canals, traffic areas, or settlement areas can be impaired [1,5]. These effects are known as off-site damage. This includes the quantitative and qualitative impairment of water bodies due to the input of soil material and the water-polluting substances bound to the solid matter, such as heavy metals or nutrients from fertilizers and pesticides [4,6,7].
In order to be able to quantify soil erosion, various models have been developed in recent decades. Among others, the Universal Soil Loss Equation (USLE) is an empirical ods and processing of datasets with spatial reference. As a result, the derivation of the individual factors of the equation becomes significantly better and more accurate depending on the quality and availability of the input data. Methods of image analysis from remote sensing data have been used in the recent past to quantify the P-factor. Karydas et al. [20] used an object-oriented classification to identify terraces based on QuickBird imagery and defined site-specific P-factors. Wang et al. [21] derived P-factor values for contouring by calculating the angle between aspect and ridge direction based on visual satellite image interpretation.
This study adopts the manual approach of Wang et al. [21] to use visual image interpretation to derive the cultivation direction for individual agricultural land parcels from remote sensing data. The main objectives of this study include: • Development of an automated method for determining the main cultivation direction on each agricultural land parcel with seasonal crop; • Site-specific verification of contouring by comparing the main cultivation direction with the mean aspect direction for each agricultural land parcel; • Improvement of the erosion risk assessment by means of USLE equations using an automated derivation of the P-factor based on an airborne image.
To automate this process of visual image interpretation, we apply existing image analysis techniques to high-resolution remote sensing data from Google Earth™ to detect linear structures [22]. A sub-meter image resolution makes it possible to deduce the tramlines of agricultural machinery or even the seed rows of crops or furrows in agricultural fields. In the following, this direction of processing on the field is defined as the cultivation direction. A comparison of this cultivation direction with the mean aspect direction for each land parcel allows the determination of a site-specific P-factor for contouring.

Study Site
To develop the method, a low mountain range area was selected because erosion processes mainly take place in sloped areas, and thus, it represents a suitable landscape area. Here, soil erosion occurs, although erosion control measures are already being implemented. Therefore, we chose the Fischbach catchment located in the southeast of Hesse, Germany, as the research area. For this local application, the German derivation of the USLE with the general soil loss equation ABAG ('Allgemeine Bodenabtragsgleichung') was used [23]. We applied the current technical specifications of the ABAG for Germany, which are normatively regulated in German standard DIN 19708 [14]. The results can be transferred to the USLE and all of its revisions.
Since 2016, the Fischbach catchment has been a part of the field observatory of the Chair of Engineering Hydrology and Water Management (ihwb) at the Technical University of Darmstadt [24][25][26][27][28]. The small catchment has a size of 35.6 km 2 and is delineated by the gauging station Groß-Bieberau 2. At this location, the ihwb monitors the discharge of the Fischbach and qualitative parameters such as water temperature, electric conductivity, pH, dissolved oxygen, and turbidity.
The altitude within the catchment ranges from 162 to 592 m.a.s.l., with a mean of 289 m.a.s.l. (see Figure 1a). The largest elevations are located in the south of the catchment area. From south to north, the elevations decrease on average. Two parallel valleys with streams running in them, Fischbach in the east and Rodauer Bach in the west, extend from the south to the north. Overall, there is a hilly topography typical for the German low mountain ranges, with slopes ranging from 0.01 to a maximum of 214.5%, with a mean slope of 18.0%. Due to this partly steep topography, a large part of the catchment is not suitable for crop cultivation with farm machinery. Therefore, the predominant land use in the Fischbach catchment, shown in Figure 1b, is forest with a share of 50.2%, followed by grassland with 19.7%. The percentage of arable land is almost the same with 19.6%. The low value of 6.7% for settlement illustrates the rural character of the catchment area. Other vegetation, orchards, and tree nurseries account for only 3.7% of the total area, and water bodies for 0.1% [29].
Due to the arable land used for agriculture on areas with steep slopes, there is an increased risk of erosion from precipitation and surface runoff in the Fischbach catchment. The P-factor is most effective between 3 and 8%, although it can be applied to croplands with slopes up to 25%, depending on the slope length [14]. Therefore, the catchment in the low mountain range is very well suited to show the possible potential of applying the P-factor. mountain ranges, with slopes ranging from 0.01 to a maximum of 214.5%, with a mean slope of 18.0%. Due to this partly steep topography, a large part of the catchment is not suitable for crop cultivation with farm machinery. Therefore, the predominant land use in the Fischbach catchment, shown in Figure 1b, is forest with a share of 50.2%, followed by grassland with 19.7%. The percentage of arable land is almost the same with 19.6%. The low value of 6.7% for settlement illustrates the rural character of the catchment area. Other vegetation, orchards, and tree nurseries account for only 3.7% of the total area, and water bodies for 0.1% [29].
Due to the arable land used for agriculture on areas with steep slopes, there is an increased risk of erosion from precipitation and surface runoff in the Fischbach catchment. The P-factor is most effective between 3 and 8%, although it can be applied to croplands with slopes up to 25%, depending on the slope length [14]. Therefore, the catchment in the low mountain range is very well suited to show the possible potential of applying the Pfactor.  [30]; (b) land use based on digital landscape model (DLM) of the official topographic-cartographic information system dataset (ATKIS) [29].
The Hessian Agency for Nature Conservation, Environment and Geology (HLNUG) has shown in a regulatory, statewide survey in Hesse that there is an enhanced risk of erosion for the agricultural land catchment area of the Fischbach river [31,32]. However, this study does not consider the P-factor. The high soil erosion risk results in particular from the brown and parabrown topsoils prevailing in the catchment area, which are susceptible to erosion due to a high loess content [31]. Furthermore, significant parameters for a high or extremely high erosion risk are strongly inclined slopes with long slope lengths and a corresponding cultivation of these areas, for example, with root crops such as corn, potatoes, or sugar beets. It is suspected that the majority of suspended solids that can be measured at the Fischbach stream gauge in Groß-Bieberau 2 are eroded from agricultural land by soil erosion and transported into the stream via connected flow paths.  [30]; (b) land use based on digital landscape model (DLM) of the official topographic-cartographic information system dataset (ATKIS) [29].
The Hessian Agency for Nature Conservation, Environment and Geology (HLNUG) has shown in a regulatory, statewide survey in Hesse that there is an enhanced risk of erosion for the agricultural land catchment area of the Fischbach river [31,32]. However, this study does not consider the P-factor. The high soil erosion risk results in particular from the brown and parabrown topsoils prevailing in the catchment area, which are susceptible to erosion due to a high loess content [31]. Furthermore, significant parameters for a high or extremely high erosion risk are strongly inclined slopes with long slope lengths and a corresponding cultivation of these areas, for example, with root crops such as corn, potatoes, or sugar beets. It is suspected that the majority of suspended solids that can be measured at the Fischbach stream gauge in Groß-Bieberau 2 are eroded from agricultural land by soil erosion and transported into the stream via connected flow paths.

Google Earth™ Remote Sensing Data
Google Earth™ provides an excellent platform for the analysis of high-resolution remote sensing images based on freely available data. In terms of image analysis, the use of the data is limited to non-commercial use and non-mass download or bulk feeds, among others [33]. Available remote sensing data on Google Earth™ consist of satellite and aircraft-based imagery with three bands of red, green, and blue (RGB), whereas the available resolution varies for the different regions of the Earth with the data provider of the imagery. Nevertheless, Google Earth™ has been used in previous studies to map objects such as gullies in the context of soil erosion [34][35][36]. Boardman et al. [37] pointed out the untapped potential of Google Earth™ imagery, which can significantly expand the existing database in the field of soil erosion based on the free availability of data and allows a retrospective view of historical events. However, extracting the georeferenced aerial imagery to process the data is inconvenient. Using a Python script developed for this purpose, we extracted images from Google Earth™ Pro App (7.3.3.7786) and used image data taken on 30 June 2019 with an approximated resolution of about 0.6 m × 0.6 m per pixel (px). Using the script, we set a buffer radius around the catchment area, divided it into a matrix, and saved each image of the matrix after resetting tilt and compass. Additionally, we stored the coordinates of the four corner points for each of the 120 images in a table for later georeferencing.

Digital Elevation Model (DEM)
We used two different digital elevation models to study the influence of raster resolution and accuracy on slope and aspect data for individual agricultural land parcels. The Hessian Administration for Land Management and Geoinformation (HVBG) derived the DEMs from airborne laser scanning data with a resolution of 5 m × 5 m (DEM5) and 1 m × 1 m (DEM1) [30]. The vertical accuracy of DEM1 is quite high at ±0.5 m, while for DEM5 this value is ±1.1 m [38].

Land Use and Land Parcel Information
The basis for determining land use classes in the catchment was the Digital Landscape Model (DLM) from the official topographic-cartographic information system dataset (ATKIS) from 2017 [29]. It contains the classified, topographic features of the surface (see Figure 1b). Since the P-factor is site-specific, the appropriate reference unit for determining the factor is the individual agricultural land parcel with one single crop. The Land Parcel Information System (LPIS) within the European Union was used to determine the eligibility of agricultural land [39]. It accurately maps agricultural land parcels, as subsidies are calculated for each farmer individually. For the definition of the individual agricultural land parcels in the catchment area, we use the clipped and anonymized dataset of the reference annual land uses (RJN) for Hesse [40]. The original dataset contains 1285 agricultural land parcels that lay completely or partly within the catchment area including their associated crops for the year 2019.

Methods
The first step of the preprocessing workflow was to define the catchment area for the gauging station Groß-Bieberau 2 with the Hydrology Toolset of ArcGIS. For the definition we only used DEM1 and filled it in advance. For the further analysis, we used the provided DEMs without prior preparation. The year 2019 was selected for analysis due to the current and complete data available in Google Earth™. In order to be able to check soil protective measures of contouring for each agricultural land parcel, the parcel information dataset was imported and analyzed in the first step of data preparation in terms of crop type and geometric shape (see Section 2.3.1). In the second step of data preparation, an analysis of the DEM was carried out and mean slope and mean aspect were determined as in Figure 2a (see Section 2.3.2). LPIS data can optionally be supplemented with land use data. Image data by Google Earth™ had to be georeferenced before the final export of individual masked images for each land parcel (see Figure 2b). After data preparation, a main cultivation direction was derived from the image data using the image analysis techniques in Figure 2c,d (see Section 2.3.3). Subsequently, these directions were compared with each other (see Figure 2e). A site-specific P-factor for the USLE could then be determined from the deviation between these direction values, taking the slope's gradient and length into account. individual masked images for each land parcel (see Figure 2b). After data preparation, a main cultivation direction was derived from the image data using the image analysis techniques in Figure 2c,d (see Section 2.3.3). Subsequently, these directions were compared with each other (see Figure 2e). A site-specific P-factor for the USLE could then be determined from the deviation between these direction values, taking the slope's gradient and length into account.

Land Parcel Geometry Analysis
Agricultural land parcels are based on farmers' property. This results in cultivation with a single crop per parcel with uniform management and growing conditions [1,42]. The topography, adjacent farmland, and the surrounding landscape features such as hedges, shrubs, conservation strips, and fallow land as well as other land uses such as forests, settlements, and water bodies or roads determine the shape of the land parcel. For economic reasons of agricultural machining, many fields resemble a parallelogram in shape, or their polygonal geometries are composed of triangular or quadrangular shapes. The time required for the cultivation decreases with the medium length of the parcel and increases with the irregularity of the shape of the parcel, as the proportion of time required for turns and headland cultivation increases [43]. Thus, from an economic point of view, the geometric shape and the topography of arable fields usually lead to the cultivation direction. However, larger arable fields usually do not have a simple rectangular shape. Moreover, the ratios of the side lengths are often not equal. To quantify the morphological shape form of the agricultural land parcels, we used the isoperimetric inequality. In the two-dimensional plane, the isoperimetric inequality for the perimeter P and the area A of a polygon states that This means that the circle is a shape that encloses the maximum area for a given perimeter [44]. Among others, Osserman [45] derived the isoperimetric quotient (IPQ) from this relationship as a shape compactness value CIPQ with CIPQ = 4πA/P 2 (2) which has become one of the most accepted factors for compactness [46]. The CIPQ is defined as the ratio of a shape area to that of a circle with the same perimeter. Thus, since a circle with a CIPQ of 1.0 is considered the most compact of all shapes, this ratio is called circularity [44]. The CIPQ of shapes ranges in the interval (0, 1] [46]. Figure 3 shows a comparison of the CIPQ values for different geometric shapes. Circle, rectangle, and right triangle are listed as examples. As the ratio of width to height for the rectangle and base to height for the triangle increases, the CIPQ decreases and finally approaches zero. The smaller the value is, the less compact is the geometric shape. The factor is rotationally α Figure 2. Workflow of land-parcel-based application of Fast Line Detector to high-resolution remote sensing data for automated derivation of main cultivation direction and verification of support practice measure contouring: (a) DEM analysis; (b) preparing remote sensing data; (c) applying FLD; (d) deriving main cultivation direction; (e) comparing of directions; data source: [30,40,41].

Land Parcel Geometry Analysis
Agricultural land parcels are based on farmers' property. This results in cultivation with a single crop per parcel with uniform management and growing conditions [1,42]. The topography, adjacent farmland, and the surrounding landscape features such as hedges, shrubs, conservation strips, and fallow land as well as other land uses such as forests, settlements, and water bodies or roads determine the shape of the land parcel. For economic reasons of agricultural machining, many fields resemble a parallelogram in shape, or their polygonal geometries are composed of triangular or quadrangular shapes. The time required for the cultivation decreases with the medium length of the parcel and increases with the irregularity of the shape of the parcel, as the proportion of time required for turns and headland cultivation increases [43]. Thus, from an economic point of view, the geometric shape and the topography of arable fields usually lead to the cultivation direction. However, larger arable fields usually do not have a simple rectangular shape. Moreover, the ratios of the side lengths are often not equal. To quantify the morphological shape form of the agricultural land parcels, we used the isoperimetric inequality. In the two-dimensional plane, the isoperimetric inequality for the perimeter P and the area A of a polygon states that P 2 ≥ 4πA.
This means that the circle is a shape that encloses the maximum area for a given perimeter [44]. Among others, Osserman [45] derived the isoperimetric quotient (IPQ) from this relationship as a shape compactness value C IPQ with C IPQ = 4πA/P 2 (2) which has become one of the most accepted factors for compactness [46]. The C IPQ is defined as the ratio of a shape area to that of a circle with the same perimeter. Thus, since a circle with a C IPQ of 1.0 is considered the most compact of all shapes, this ratio is called circularity [44]. The C IPQ of shapes ranges in the interval (0, 1] [46]. Figure 3 shows a comparison of the C IPQ values for different geometric shapes. Circle, rectangle, and right triangle are listed as examples. As the ratio of width to height for the rectangle and base to height for the triangle increases, the C IPQ decreases and finally approaches zero. The smaller the value is, the less compact is the geometric shape. The factor is rotationally invariant and, by squaring the perimeter, independent of the variable side lengths of the polygons under study [46].

Digital Elevation Model Analysis
Analysis of the digital elevation models was performed using Esri ArcGIS Pro (2.7.3). Calculation of slope and aspect for the whole catchment was conducted by planar method with a 3 × 3 neighborhood for DEM5 and a 9 × 9 neighborhood for DEM1 to consider highresolution elevation data [47]. The slope describes the relationship between horizontal elevation changes on a defined vertical distance. Based on the raster dataset, the slope is defined as medium horizontal inclination of the surrounding raster pixels in the given neighborhood. The aspect specifies the compass direction facing this average downhill slope. In the context of evaluation, we define the zonal mean value of the slope of a land parcel as the main slope for this parcel. The parameter aspect follows the same procedure. Especially in hilly and mountainous areas, the derived mean values of slope and aspect for a plot or parcel do not always correspond to reality since the terrain may incline in several directions and with different gradients. Therefore, we considered the standard deviation (SD) as an index for the complexity of the topography with different slopes or directions. The smaller the value for the SD is, the smaller the dispersion is of the measured features from the mean and the more uniformly the terrain within the plot is sloped in direction considering aspect and in gradient considering slope. Mean and SD for slope raster were figured out with zonal statistics for each land parcel. For the aspect grid, we had to consider the circular statistical functions of the Python package SciPy to determine the mean and SD of the aspect values [48]. For all indicated aspect angles, we used the north azimuth system, measured in clockwise direction from north in units of degrees in the interval [0, 360).

Image Analysis
The usual cultivation of fields by agricultural machinery is carried out in strips, so cultivation of the seed rows and tramlines (tire tracks) runs parallel to the long sides of the field geometries. Tramlines are caused due to regular use for maintenance purposes, e.g., applying pesticides or fertilizer during the growth phase of the crop. In the context of topology, the tramlines form a continuous, mostly parallel network across the entire field and do not end abruptly except at the edge of the field. Usually, the tramlines are not sown. If they are sown, the plant growth within the tramlines is negatively influenced by

Digital Elevation Model Analysis
Analysis of the digital elevation models was performed using Esri ArcGIS Pro (2.7.3). Calculation of slope and aspect for the whole catchment was conducted by planar method with a 3 × 3 neighborhood for DEM5 and a 9 × 9 neighborhood for DEM1 to consider high-resolution elevation data [47]. The slope describes the relationship between horizontal elevation changes on a defined vertical distance. Based on the raster dataset, the slope is defined as medium horizontal inclination of the surrounding raster pixels in the given neighborhood. The aspect specifies the compass direction facing this average downhill slope. In the context of evaluation, we define the zonal mean value of the slope of a land parcel as the main slope for this parcel. The parameter aspect follows the same procedure. Especially in hilly and mountainous areas, the derived mean values of slope and aspect for a plot or parcel do not always correspond to reality since the terrain may incline in several directions and with different gradients. Therefore, we considered the standard deviation (SD) as an index for the complexity of the topography with different slopes or directions. The smaller the value for the SD is, the smaller the dispersion is of the measured features from the mean and the more uniformly the terrain within the plot is sloped in direction considering aspect and in gradient considering slope. Mean and SD for slope raster were figured out with zonal statistics for each land parcel. For the aspect grid, we had to consider the circular statistical functions of the Python package SciPy to determine the mean and SD of the aspect values [48]. For all indicated aspect angles, we used the north azimuth system, measured in clockwise direction from north in units of degrees in the interval [0, 360).

Image Analysis
The usual cultivation of fields by agricultural machinery is carried out in strips, so cultivation of the seed rows and tramlines (tire tracks) runs parallel to the long sides of the field geometries. Tramlines are caused due to regular use for maintenance purposes, e.g., applying pesticides or fertilizer during the growth phase of the crop. In the context of topology, the tramlines form a continuous, mostly parallel network across the entire field and do not end abruptly except at the edge of the field. Usually, the tramlines are not sown. If they are sown, the plant growth within the tramlines is negatively influenced by soil compaction, and the growth progress within the tramlines lags behind that of the other plants [49]. Due to the difference in albedo between plants and soil as well as increasing shading within the tramlines as the plants grow taller, texture gradients are created on the remote sensing images that can be deduced using computer vision [50]. However, for perception of these gradients, a corresponding resolution in the sub-meter range is required for the remote sensing dataset. Therefore, freely available remote sensing data from satellites with a current resolution of 10 m × 10 m (Sentinel-2) or 15 m × 15 m (Landsat 8) per pixel are not suitable [51,52]. Moreover, the intensity of the texture gradient depends on the crop type and the phenological growth stages of the plants [53].
To link the image data to the arable land parcels, the images extracted from Google Earth™ were first loaded into ArcGIS Pro and georeferenced using the four saved corner points for each image. Subsequently, all images were merged to form a mosaic with histogram adaption. Then, we extracted the image grid for each arable land parcel by mask. To obtain a suitable minimum resolution for the images of each land parcel, we only considered areas larger than 1500 m 2 . For further image analysis, we used the Open Source Computer Vision (OpenCV) library with version 4.0.1 for Python [54,55]. The OpenCV library is open source and free to use under license [56]. Among others, it contains the Fast Line Detector algorithm by Lee et al. [22] that works with grayscale images and shows good results for the detection of linear segments [57].
The FLD algorithm uses the Canny edge detector for the effective identification of edges, which operates in the following four steps [50,59,60]:

1.
Image smoothing by Gaussian filter for noise reduction; 2.
Computing local gradient with Sobel operator; 3.
Hysteresis threshold to form contours.
To improve the result of the FLD algorithm, the input parameters for line detection listed in Table 1 can be adjusted. For all subsequent applications of the FLD algorithm, the parameter 'Do Merge' was set to 'False' in order to prevent the merging of segments, which may promote the connection of detected lines of different categories-for example, detected lines from soil erosion and tramlines. In addition, the Canny Aperture Size was set to 5 since the value 3 was usually responsible for the FLD detecting significantly too few lines in preliminary investigations. All other parameters will be further investigated in the adjustment of the FLD in Section 3.3.1. To ensure that the boundaries of the land parcel are not detected as a line, each masked raster image was buffered before line detection. For this operation, a width for the field border should be selected depending on the field area and the image resolution. Especially with a large width-height ratio, the total width of the field should be taken into account so that the buffer width is not too large. The selected buffer width was between 15 and 30 pixels, which corresponds to approximately 9 to 18 m. As a result, landscape features such as individual trees or other vegetation located at the field boundaries or extending the boundaries of the land parcel did not impair line detection. In addition, a part of the headland was not taken into account in this way. The algorithm output the results both visually in the form of an image with detected lines and in the form of a table in which the vertex points of the detected lines are stored in the form of image coordinates. In order Land 2021, 10, 1279 9 of 34 to take into account the continuous topology of the tramlines when evaluating the main cultivation direction, we calculated a weighted average considering the length of the lines. For specifying the main cultivation direction, we used the north azimuth system. Since agricultural machinery can drive along the tramlines in both directions, the interval was limited to angles x with x ∈ [0, 180). All angles x with x ≥ 180 • were converted accordingly by subtracting 180 • . The main cultivation direction for each agricultural land parcel was determined visually in a manual manner using the Google Earth™ dataset in order to be able to check the automatically derived main cultivation direction and calibrate the parameters for FLD algorithm. In addition, we categorized all images regarding whether the cultivation direction was definite or ambiguous by subjective decision. To quantify the result of the FLD algorithm, the mean absolute error (MAE) was applied. It compared the absolute difference between the automatically detected cultivation direction x i and the main cultivation direction visually derived in a manual manner y i for each field. When applied to multiple parameter constellations, the mean was calculated by dividing the sum of all absolute differences by the number n of constellations: A maximum MAE of less than 5 • was defined as a successful detection of the main cultivation direction by the FLD.

Cultivation Statistics
The analysis of different agricultural land uses in 2019 for Germany in Figure 4 shows that arable land accounts for by far the largest share at 70%. For Hesse, the share of arable land is somewhat lower at 61%, and with only 38.5%, this share is significantly lower for the Fischbach catchment. The shares of permanent crops are constant for all regions, being negligible at about 1%. Thus, permanent grassland takes the inverse share of arable land for all regions. As the area of the study region under consideration decreases, the proportion of cropland decreases, while the proportion of permanent grassland increases at the same rate. Thus, the Fischbach catchment eventually had a distinct inverse ratio of arable land to permanent grassland of 60.0%. 15 and 30 pixels, which corresponds to approximately 9 to 18 m. As a result, landscape features such as individual trees or other vegetation located at the field boundaries or extending the boundaries of the land parcel did not impair line detection. In addition, a part of the headland was not taken into account in this way. The algorithm output the results both visually in the form of an image with detected lines and in the form of a table in which the vertex points of the detected lines are stored in the form of image coordinates. In order to take into account the continuous topology of the tramlines when evaluating the main cultivation direction, we calculated a weighted average considering the length of the lines. For specifying the main cultivation direction, we used the north azimuth system. Since agricultural machinery can drive along the tramlines in both directions, the interval was limited to angles x with x ∈ [0, 180). All angles x with x ≥ 180° were converted accordingly by subtracting 180°. The main cultivation direction for each agricultural land parcel was determined visually in a manual manner using the Google Earth™ dataset in order to be able to check the automatically derived main cultivation direction and calibrate the parameters for FLD algorithm. In addition, we categorized all images regarding whether the cultivation direction was definite or ambiguous by subjective decision. To quantify the result of the FLD algorithm, the mean absolute error (MAE) was applied. It compared the absolute difference between the automatically detected cultivation direction xi and the main cultivation direction visually derived in a manual manner yi for each field. When applied to multiple parameter constellations, the mean was calculated by dividing the sum of all absolute differences by the number n of constellations: A maximum MAE of less than 5° was defined as a successful detection of the main cultivation direction by the FLD.

Cultivation Statistics
The analysis of different agricultural land uses in 2019 for Germany in Figure 4 shows that arable land accounts for by far the largest share at 70%. For Hesse, the share of arable land is somewhat lower at 61%, and with only 38.5%, this share is significantly lower for the Fischbach catchment. The shares of permanent crops are constant for all regions, being negligible at about 1%. Thus, permanent grassland takes the inverse share of arable land for all regions. As the area of the study region under consideration decreases, the proportion of cropland decreases, while the proportion of permanent grassland increases at the same rate. Thus, the Fischbach catchment eventually had a distinct inverse ratio of arable land to permanent grassland of 60.0%.   The above distribution shows a different relationship for permanent grassland and arable land than the ATKIS land use data in Section 2.2.3, even though the reference area is of a different size. The share for permanent grassland is significantly higher in the RJN dataset.
Arable land grows crops of various categories, such as cereals, plants for green harvesting (e.g., silage maize), root crops, legumes, oil crops and other cash crops, and horticultural crops [62]. Figure 5 shows the distribution of area shares for the respective categories in comparison for Germany, Hesse, and the Fischbach catchment in 2019. Cereals take the largest share of arable land in all regions at over 50%. In Hesse, this share is 67.3%, and in the Fischbach catchment it is higher at 71.1%. Crops for green harvesting, on the other hand, account for about 15% in Hesse and the Fischbach catchment, which is significantly below the average value of 25.0% for Germany. The shares for root crops are about the same for all regions, at about 5%. The same applies to legumes at about 2%. The area share for oil crops in the Fischbach catchment of 2.0% is significantly lower than those in Hesse (6.3%) and Germany (8.0%). Horticultural crops are significantly less represented in the Fischbach catchment at 0.2%. The share for eligible set-aside land is slightly higher in the study area compared to Hesse and Germany. Even though the proportion of arable land in the Fischbach catchment is lower than in Hesse and Germany, the relative distribution of crop categories is very typical for Hesse and comparable to the average proportions in Germany.
arable land than the ATKIS land use data in Section 2.2.3, even though the reference area is of a different size. The share for permanent grassland is significantly higher in the RJN dataset.
Arable land grows crops of various categories, such as cereals, plants for green harvesting (e.g., silage maize), root crops, legumes, oil crops and other cash crops, and horticultural crops [62]. Figure 5 shows the distribution of area shares for the respective categories in comparison for Germany, Hesse, and the Fischbach catchment in 2019. Cereals take the largest share of arable land in all regions at over 50%. In Hesse, this share is 67.3%, and in the Fischbach catchment it is higher at 71.1%. Crops for green harvesting, on the other hand, account for about 15% in Hesse and the Fischbach catchment, which is significantly below the average value of 25.0% for Germany. The shares for root crops are about the same for all regions, at about 5%. The same applies to legumes at about 2%. The area share for oil crops in the Fischbach catchment of 2.0% is significantly lower than those in Hesse (6.3%) and Germany (8.0%). Horticultural crops are significantly less represented in the Fischbach catchment at 0.2%. The share for eligible set-aside land is slightly higher in the study area compared to Hesse and Germany. Even though the proportion of arable land in the Fischbach catchment is lower than in Hesse and Germany, the relative distribution of crop categories is very typical for Hesse and comparable to the average proportions in Germany. Since the further analysis will examine the cultivation direction in fields with seasonal crops, eligible set-aside land will not be considered anymore. These include buffer strips, ornamental plants, flowering strips, and honey plants as well as fallow and arable land taken from production, as no clear plant type is defined for these areas. In addition, crops from the category of plants for green harvesting are neglected, such arable grass, clover grass, and alfalfa, even if planted seasonally, except maize for silage.
After data processing, the dataset in Figure 6 contains the final geometries for the total of 300 land parcels of arable land bigger than 1500 m 2 that lay completely or partly within the catchment area including their associated crops for the year 2019. The remaining arable land parcels are located especially in the flatter north and, furthermore, in the western part of the catchment.  Since the further analysis will examine the cultivation direction in fields with seasonal crops, eligible set-aside land will not be considered anymore. These include buffer strips, ornamental plants, flowering strips, and honey plants as well as fallow and arable land taken from production, as no clear plant type is defined for these areas. In addition, crops from the category of plants for green harvesting are neglected, such arable grass, clover grass, and alfalfa, even if planted seasonally, except maize for silage.
After data processing, the dataset in Figure 6 contains the final geometries for the total of 300 land parcels of arable land bigger than 1500 m 2 that lay completely or partly within the catchment area including their associated crops for the year 2019. The remaining arable land parcels are located especially in the flatter north and, furthermore, in the western part of the catchment. Table 2 shows an overview of the absolute and relative area shares of the individual crop types in the Fischbach catchment. In addition, the number of land parcels with these crop types is listed there. Winter durum and winter barley together accounted for more than 50% of the total area and, with 157 parcels, more than half of all accounted arable land parcels. Maize plants were subdivided according to the time of harvest, based on processing of the harvested plant in the form of silage or grain corn (category root crops). Combined, they accounted for 15.4% of the total area or 44 land parcels, with each type accounting for one half. Other cereal crops followed, including winter triticale (8.0%), summer oats (4.6%), spring barley (2.6%), and winter mixed grains (1.2%). Additionally, 2.3% contained oil crops such as winter or summer rapeseed. Legumes included broad bean and soybean with a rather low share of 1.5 and 0.8%, respectively. Root crops included sugar beets (4.8%) and potatoes (0.5%). In addition, there were individual fields (1 to 3) with horticultural plants (vegetables, giant pumpkin), other cereals, and other crops, represented by the pseudocereal buckwheat.

Geometry Shape Analysis
Since the minimum area size was limited to 1500 m 2 , the size of the land parcels ranged from 1551 m 2 to a maximum of 73,936 m 2 , with a mean value of 13,745 m 2 . The median of 11.358 m 2 , which is significantly smaller than the mean, resulted in a rightskewed distribution for the shape area of the land parcels. This is evident in the cumulative frequency analysis in Figure 7. The Fischbach catchment contains a very large number of land parcels with a very small area <20,000 m 2 (78.7%). In addition, 19.7% of small-sized parcels are from 20,000 up to 40,000 m 2 in area. Larger parcels bigger than 40,000 m 2 are represented only five times, of which two parcels are larger than 60,000 m 2 . Compared to Hesse with a mean size of 18,423 m 2 per arable land parcel, there are more very small parcels (<20,000 m 2 ) and less small parcels (20,000-40,000 m 2 ) in the Fischbach catchment. Medium (60,000-100,000 m 2 ) and large parcels (>100,000 m 2 ) are clearly underrepresented.   Table 2 shows an overview of the absolute and relative area shares of the individual crop types in the Fischbach catchment. In addition, the number of land parcels with these crop types is listed there. Winter durum and winter barley together accounted for more than 50% of the total area and, with 157 parcels, more than half of all accounted arable land parcels. Maize plants were subdivided according to the time of harvest, based on processing of the harvested plant in the form of silage or grain corn (category root crops).

Geometry Shape Analysis
Since the minimum area size was limited to 1500 m 2 , the size of the land parcels ranged from 1551 m 2 to a maximum of 73,936 m 2 , with a mean value of 13,745 m 2 . The median of 11.358 m 2 , which is significantly smaller than the mean, resulted in a rightskewed distribution for the shape area of the land parcels. This is evident in the cumulative frequency analysis in Figure 7. The Fischbach catchment contains a very large number of land parcels with a very small area <20,000 m 2 (78.7%). In addition, 19.7% of small-sized parcels are from 20,000 up to 40,000 m 2 in area. Larger parcels bigger than 40,000 m 2 are represented only five times, of which two parcels are larger than 60,000 m 2 . Compared to Hesse with a mean size of 18,423 m 2 per arable land parcel, there are more very small parcels (<20,000 m 2 ) and less small parcels (20,000-40,000 m 2 ) in the Fischbach catchment. Medium (60,000-100,000 m 2 ) and large parcels (>100,000 m 2 ) are clearly underrepresented. The plot in Figure 7 shows the compactness of each parcel, calculated according to the isoperimetric quotient CIPQ in Equation (2) as a function of the parcel area size. The results for compactness rank between 0.147 and 0.864. Values smaller than 0.3 occur only for small fields (<6000 m 2 ). In this size range of the fields, moreover, the variance of the possible values for compactness is very large. With increasing area, the minimum compactness increases. If only the compactness and its frequency are considered, a clear leftskewed distribution can be observed with a mean value of 0.625. Only 9.0% of the arable land parcels have a CIPQ greater than 0.785 and, thus, a more compact or circle-like shape than a square (see Figure 3). Some arable land parcels have holes within the geometry, with the area of the holes ranging between 24.98 and 5411.91 m 2 . As the area of the holes increases, the compactness factor decreases due to the increasing circumference and simultaneously decreasing total area (see Figure 7). In these cases, the holes usually represent landscape features such as hedges or shrubs that are not counted as agricultural land. Sometimes there is an unfavorable land parcel definition where one parcel lies within another (see Figure 8a). However, the determination of CIPQ does not work for polygons with holes. These parcels must be evaluated separately, since CIPQ takes into account the perimeter of all outer lines of the polygon, including inner edges. In addition, holes of a certain size have a particular influence on the effort for agricultural cultivation of the parcels, with the result that it is no longer possible to give an indication of the amount of cultivation required on the basis of the geometric shape. There are, altogether, four polygons with holes larger than 200 m 2 , which are separately illustrated in Figure 7. The CIPQ for these parcels is significantly lower depending on the size of the hole. The graphical analysis shows this characteristic especially for the two plots with larger holes. Figure 8 illustrates four examples of geometric shapes of arable land parcels in the Fischbach catchment with increasing values for CIPQ from (a) to (d). In Figure 8a, we first present an example of a polygon with a hole, which gives the very low CIPQ value of 0.203. The parcel in Figure  8b has a trapezoidal shape with a strongly unequal ratio of width to height (1:9). This  The plot in Figure 7 shows the compactness of each parcel, calculated according to the isoperimetric quotient C IPQ in Equation (2) as a function of the parcel area size. The results for compactness rank between 0.147 and 0.864. Values smaller than 0.3 occur only for small fields (<6000 m 2 ). In this size range of the fields, moreover, the variance of the possible values for compactness is very large. With increasing area, the minimum compactness increases. If only the compactness and its frequency are considered, a clear left-skewed distribution can be observed with a mean value of 0.625. Only 9.0% of the arable land parcels have a C IPQ greater than 0.785 and, thus, a more compact or circle-like shape than a square (see Figure 3). Some arable land parcels have holes within the geometry, with the area of the holes ranging between 24.98 and 5411.91 m 2 . As the area of the holes increases, the compactness factor decreases due to the increasing circumference and simultaneously decreasing total area (see Figure 7). In these cases, the holes usually represent landscape features such as hedges or shrubs that are not counted as agricultural land. Sometimes there is an unfavorable land parcel definition where one parcel lies within another (see Figure 8a). However, the determination of C IPQ does not work for polygons with holes. These parcels must be evaluated separately, since C IPQ takes into account the perimeter of all outer lines of the polygon, including inner edges. In addition, holes of a certain size have a particular influence on the effort for agricultural cultivation of the parcels, with the result that it is no longer possible to give an indication of the amount of cultivation required on the basis of the geometric shape. There are, altogether, four polygons with holes larger than 200 m 2 , which are separately illustrated in Figure 7. The C IPQ for these parcels is significantly lower depending on the size of the hole. The graphical analysis shows this characteristic especially for the two plots with larger holes. Figure 8 illustrates four examples of geometric shapes of arable land parcels in the Fischbach catchment with increasing values for C IPQ from (a) to (d). In Figure 8a, we first present an example of a polygon with a hole, which gives the very low C IPQ value of 0.203. The parcel in Figure 8b has a trapezoidal shape with a strongly unequal ratio of width to height (1:9). This results in a low value of 0.252. The width-height ratio of the parcel in Figure 8c is between 1:1 and 1:2. The partially slightly rounded corners of the polygon have a beneficial effect on the C IPQ , which is higher at 0.777. The convex polygon in Figure 8d resembles the circular shape even more, and the C IPQ is 0.857. Convex shapes increase the C IPQ value, while concave shapes decrease it.

DEM Analysis
Regarding the analysis of the digital elevation model, the results include a comparison of DEM1 and DEM5 for slope and aspect, spatially differentiating between the entire Fischbach catchment and the masked area of the arable land parcels. Table 3 compares the results for the slope of the statistical parameters of minimum, maximum, mean, and SD. DEM1 has a significantly higher variance when the entire Fischbach catchment is considered. Furthermore, even when using a larger neighborhood for the slope calculation, the maximum slope for DEM1 is twice as large (439.91%) as that for DEM5 (214.50%). The high maximum slope values result from the edges of an abandoned quarry in the study area. This results in a larger mean for DEM1 and a higher SD in relation to DEM5. However, the differences in slope calculated with DEM1 and DEM5 are only apparent in detail, which is why Figure 9 shows only the slope from DEM5 as an example. Comparison of the land use in Figure 1b with the slope in Figure 9 indicates that areas with steep slopes are predominantly forested or covered with permanent grassland, so significantly lower maximum values of about one-fifth of the total catchment can be observed if only the masked areas of the arable land parcels are considered. This is true for both DEM1 and DEM5. The mean and SD values are significantly lower, too. Furthermore, the difference in statistical values between the two different DEM input datasets is smaller.

DEM Analysis
Regarding the analysis of the digital elevation model, the results include a comparison of DEM1 and DEM5 for slope and aspect, spatially differentiating between the entire Fischbach catchment and the masked area of the arable land parcels. Table 3 compares the results for the slope of the statistical parameters of minimum, maximum, mean, and SD. DEM1 has a significantly higher variance when the entire Fischbach catchment is considered. Furthermore, even when using a larger neighborhood for the slope calculation, the maximum slope for DEM1 is twice as large (439.91%) as that for DEM5 (214.50%). The high maximum slope values result from the edges of an abandoned quarry in the study area. This results in a larger mean for DEM1 and a higher SD in relation to DEM5. However, the differences in slope calculated with DEM1 and DEM5 are only apparent in detail, which is why Figure 9 shows only the slope from DEM5 as an example. Comparison of the land use in Figure 1b with the slope in Figure 9 indicates that areas with steep slopes are predominantly forested or covered with permanent grassland, so significantly lower maximum values of about one-fifth of the total catchment can be observed if only the masked areas of the arable land parcels are considered. This is true for both DEM1 and DEM5. The mean and SD values are significantly lower, too. Furthermore, the difference in statistical values between the two different DEM input datasets is smaller.  Figure 9a shows that extremely steep slopes are not found on the arable land parcels. For the masked arable land parcels in Figure 9b, a mean value was derived from the dataset in Figure 9a and classified according to the P-factor determination in ABAG [14]. The calculated mean values range from 1.52 to 21.59%. Parcels with lower mean values are located in the valleys of the flowing waters and at the edge of the catchment on the upper plateaus. In contrast, parcels with steeper mean slopes are located on the slopes between the outer plateaus and the streams. The spatial distribution of slope values as a function of location within the watershed cannot be determined. plateaus. In contrast, parcels with steeper mean slopes are located on the slopes between the outer plateaus and the streams. The spatial distribution of slope values as a function of location within the watershed cannot be determined.

Slope Analysis
(a) (b) Figure 9. Slope in Fischbach catchment including area of arable land parcels under investigation: (a) Slope DEM5; (b) mean slope of DEM5 for masked arable land parcels following classification of ABAG for determining the P-factor; data source: [14,30,40].
Using the zonal statistics method (see Section 2.3.1), the differences between the parameters for DEM1 and DEM5 in Table 3 were found to be very small. This is true for the population of all elements as well as for each individual arable land parcel in Figure 10.
Here, in ascending order, all mean values for slope per land parcel were calculated based on DEM5 on the x-axis and based on DEM1 on the y-axis. The figure shows that for 17 out of 300 land parcels (5.7%) the absolute difference in mean slope is greater than 0.5%. These points are marked with a red outline in Figure 10. The maximum value for absolute difference is 1.42%. For all other parcels the absolute difference is smaller than 0.5%. Absolute differences greater than 0.5% are distributed over the entire range of mean slope values and are not limited to a low, medium, or high slope. (b) mean slope of DEM5 for masked arable land parcels following classification of ABAG for determining the P-factor; data source: [14,30,40].
Using the zonal statistics method (see Section 2.3.1), the differences between the parameters for DEM1 and DEM5 in Table 3 were found to be very small. This is true for the population of all elements as well as for each individual arable land parcel in Figure 10.
Here, in ascending order, all mean values for slope per land parcel were calculated based on DEM5 on the x-axis and based on DEM1 on the y-axis. The figure shows that for 17 out of 300 land parcels (5.7%) the absolute difference in mean slope is greater than 0.5%. These points are marked with a red outline in Figure 10. The maximum value for absolute difference is 1.42%. For all other parcels the absolute difference is smaller than 0.5%. Absolute differences greater than 0.5% are distributed over the entire range of mean slope values and are not limited to a low, medium, or high slope.

Aspect Analysis
The orientation of the land slope in the catchment is based on the streams that run in the valleys. This characteristic orientation is largely evident for the calculated mean aspect per parcel from DEM5 in Figure 11, where the mean slope direction for each arable land parcel is indicated by an arrow that points downslope. As the distance to the stream increases, the mean slope direction no longer necessarily points to the stream direction. This is especially true for the land parcels located at the border of the catchment on upper plateaus. However, as already indicated in Section 2.1, the topography, especially with respect to the aspect, is often more complex in low mountain ranges and the parcels are often sloped in several different directions. This feature is illustrated in Figure 11 with blue color coding for the SD. The darker the blue shade of the arrows is, the higher is the SD for the respective parcel. Areas with a high SD are located in particular on the upper plateaus at the border of the catchment, although not exclusively there.

Aspect Analysis
The orientation of the land slope in the catchment is based on the streams that run in the valleys. This characteristic orientation is largely evident for the calculated mean aspect per parcel from DEM5 in Figure 11, where the mean slope direction for each arable land parcel is indicated by an arrow that points downslope. As the distance to the stream increases, the mean slope direction no longer necessarily points to the stream direction. This is especially true for the land parcels located at the border of the catchment on upper plateaus. However, as already indicated in Section 2.1, the topography, especially with respect to the aspect, is often more complex in low mountain ranges and the parcels are often sloped in several different directions. This feature is illustrated in Figure 11 with blue color coding for the SD. The darker the blue shade of the arrows is, the higher is the SD for the respective parcel. Areas with a high SD are located in particular on the upper plateaus at the border of the catchment, although not exclusively there.
The SD for all arable land parcels in the catchment is significantly increased, suggesting a complex topography and thus a large proportion of surfaces that are sloped in multiple directions. Comparing the input datasets of DEM1 and DEM5, the correlation for the mean aspect per parcel is even more significant than that for the slope determination when the limit value for an acceptable deviation is set to 5.0 • . This results in only 4 out of 300 land parcels, marked with a red outline in Figure 12, where the absolute difference in mean aspect is greater than 5.0 • . The SD for all arable land parcels in the catchment is significantly increased, suggesting a complex topography and thus a large proportion of surfaces that are sloped in multiple directions. Comparing the input datasets of DEM1 and DEM5, the correlation for the mean aspect per parcel is even more significant than that for the slope determination when the limit value for an acceptable deviation is set to 5.0°. This results in only 4 out of 300 land parcels, marked with a red outline in Figure 12, where the absolute difference in mean aspect is greater than 5.0°. Figure 11. Mean aspect of DEM5 for masked arable land parcels including SD in relation to streams in the Fischbach catchment; data source: [29,30,40]. The final comparison of aspect and slope in Figure 13 indicates a weak linear correlation between mean slope and SD of aspect. At low mean slope <8%, the SD of aspect can reach values above 100°. With increasing values for mean slope, the SD value for aspect decreases and thus indicates the complexity of the topography with respect to the aspect. As a result, the steeper sloped land parcels usually have a more uniform direction of slope.  The final comparison of aspect and slope in Figure 13 indicates a weak linear correlation between mean slope and SD of aspect. At low mean slope <8%, the SD of aspect can reach values above 100 • . With increasing values for mean slope, the SD value for aspect decreases and thus indicates the complexity of the topography with respect to the aspect. As a result, the steeper sloped land parcels usually have a more uniform direction of slope. Figure 12. Comparison of mean aspect per arable parcel for DEM1 and DEM5. Red points indicate an absolute difference greater than 5.0°; data source: [30,40].
The final comparison of aspect and slope in Figure 13 indicates a weak linear correlation between mean slope and SD of aspect. At low mean slope <8%, the SD of aspect can reach values above 100°. With increasing values for mean slope, the SD value for aspect decreases and thus indicates the complexity of the topography with respect to the aspect. As a result, the steeper sloped land parcels usually have a more uniform direction of slope. Figure 13. Comparison of mean slope per arable parcel for DEM5 and SD of aspect for DEM5; data source: [30,40].

Image Analysis
In the initial application, the FLD algorithm was designed to determine the main cultivation direction of the arable land parcels in a simplified way. Among the 300 images of arable land parcels extracted from Google Earth™, there were 18 images with an ambiguous main cultivation direction. Therefore, these were not considered further in the analysis. Similarly, the four arable land parcels with holes in their geometry (see Section 3.1.2) were not included in the initial application. The remaining 278 masked arable land parcels with a definite main cultivation direction and without holes were categorized by crop type. Initially, the physical reasonable range of the parameters Length Threshold t, Distance Threshold dt, Canny Threshold 1 cth1, and Canny Threshold 2 cth2 was assessed.

Image Analysis
In the initial application, the FLD algorithm was designed to determine the main cultivation direction of the arable land parcels in a simplified way. Among the 300 images of arable land parcels extracted from Google Earth™, there were 18 images with an ambiguous main cultivation direction. Therefore, these were not considered further in the analysis. Similarly, the four arable land parcels with holes in their geometry (see Section 3.1.2) were not included in the initial application. The remaining 278 masked arable land parcels with a definite main cultivation direction and without holes were categorized by crop type. Initially, the physical reasonable range of the parameters Length Threshold t, Distance Threshold dt, Canny Threshold 1 cth1, and Canny Threshold 2 cth2 was assessed. Then, the parameters were assigned to achieve the best line detection results. For both procedures, the image dataset of the crop type spring barley with its eight images was used because it contained images with well-visible and less well-visible tramlines. Furthermore, it included images with partly extending trees and their shadows at the parcel border. Additionally, we used the smallest arable land parcel of the investigation with the crop type winter oats. Finally, the results for the application of the FLD with the parameters applied to all 278 images and parameter constellations for each crop type are shown.

Adjustment of FLD Parameters
First, we investigated the influence of the individual parameters on the detection of single lines as well as their average orientation. For this purpose, we varied the value for the parameter Length Threshold t in Figure 14 for two arable land parcels with different area sizes while keeping the other parameters constant. An increasing value for t filtered out shorter lines, so the average length of the detected lines increased for an increasing t. Filtering out smaller lines ensured that fewer lines would be detected by the FLD in total. This applied to all images regardless of shape and size.
The variation of the Length Threshold t in Figure 14 shows that lines less than or equal to 10 px (6 m) have a significantly larger variance in the orientation. Comparing the area of the parcels, this variance is more pronounced for larger fields (see bottom row of Figure 14) than for smaller fields (see upper row of Figure 14). This behavior is evident in the analysis of the SD for the detected lines in Figure 15. With increasing t, the SD of all detected lines decreases. The SD for small fields is generally lower than the SD for large fields. Therefore, the minimum value for t was set to 10 px (6 m).
single lines as well as their average orientation. For this purpose, we varied the value for the parameter Length Threshold t in Figure 14 for two arable land parcels with different area sizes while keeping the other parameters constant. An increasing value for t filtered out shorter lines, so the average length of the detected lines increased for an increasing t. Filtering out smaller lines ensured that fewer lines would be detected by the FLD in total. This applied to all images regardless of shape and size. The variation of the Length Threshold t in Figure 14 shows that lines less than or equal to 10 px (6 m) have a significantly larger variance in the orientation. Comparing the area of the parcels, this variance is more pronounced for larger fields (see bottom row of Figure 14) than for smaller fields (see upper row of Figure 14). This behavior is evident in the analysis of the SD for the detected lines in Figure 15. With increasing t, the SD of all detected lines decreases. The SD for small fields is generally lower than the SD for large fields. Therefore, the minimum value for t was set to 10 px (6 m). However, for a value of t greater than 50 px, lines tend not to be detected anymore because the gradients in image values are too large (see Figure 14e). This value is correspondingly smaller for parcels with a smaller area due to the reduced number of pixels. Due to this relationship, it is suggested that t (px) is best determined as a function of field area according to the following Equation (5): Afterwards, the calculated value was converted to an integer as input for the FLD algorithm. Here, 'fieldarea' is specified as a number of pixels in order to be able to define t independently of the image resolution. Figure 16 shows the graph of the function with converted integers. Results that are smaller or larger than the two threshold values take the corresponding lower value of t = 10 px or upper value of t = 50 px. ngth Threshold t (px) Figure 15. SD for orientation of detected lines depending on Length Threshold t for one small (green, 1421) and one medium-large parcel (blue, 1324); data source: [40,41].
However, for a value of t greater than 50 px, lines tend not to be detected anymore because the gradients in image values are too large (see Figure 14e). This value is correspondingly smaller for parcels with a smaller area due to the reduced number of pixels. Due to this relationship, it is suggested that t (px) is best determined as a function of field area according to the following Equation (5): Afterwards, the calculated value was converted to an integer as input for the FLD algorithm. Here, 'fieldarea' is specified as a number of pixels in order to be able to define t independently of the image resolution. Figure 16 shows the graph of the function with converted integers. Results that are smaller or larger than the two threshold values take the corresponding lower value of t = 10 px or upper value of t = 50 px.
However, this behavior of not detecting lines for large t values depends on the Distance Threshold dt, among others. The smaller the threshold value for dt is, the more likely the detection of a line is interrupted. Figure 17 visualizes the bivariate relationship between the combination of parameters t and dt exemplified on an image with spring barley (Parcel 1323). On the horizontal axis, dt increases from left to right. The value for t increases on the vertical axis from top to bottom. Afterwards, the calculated value was converted to an integer as input for the FLD algorithm. Here, 'fieldarea' is specified as a number of pixels in order to be able to define t independently of the image resolution. Figure 16 shows the graph of the function with converted integers. Results that are smaller or larger than the two threshold values take the corresponding lower value of t = 10 px or upper value of t = 50 px. However, this behavior of not detecting lines for large t values depends on the Distance Threshold dt, among others. The smaller the threshold value for dt is, the more likely the detection of a line is interrupted. Figure 17 visualizes the bivariate relationship between the combination of parameters t and dt exemplified on an image with spring barley (Parcel 1323). On the horizontal axis, dt increases from left to right. The value for t increases on the vertical axis from top to bottom.   As already shown in Figure 14, the number of detected lines decreases, the average length of lines increases, and the SD for the orientation of the lines decreases as the value for t increases in Figure 17. These trends can be observed for different values of dt. All three trends are strongest for an increasing t value for dt = 2 1/2 px. The FLD algorithm already detected an SD of 1.6° at t = 20 px, for dt = 2 1/2 px. For all larger values of dt, the SD is clearly higher at 14.7° or even 16.0°. For dt > 2 1/2 px, the FLD detected significantly  As already shown in Figure 14, the number of detected lines decreases, the average length of lines increases, and the SD for the orientation of the lines decreases as the value for t increases in Figure 17. These trends can be observed for different values of dt. All three trends are strongest for an increasing t value for dt = 2 1/2 px. The FLD algorithm already detected an SD of 1.6 • at t = 20 px, for dt = 2 1/2 px. For all larger values of dt, the SD is clearly higher at 14.7 • or even 16.0 • . For dt > 2 1/2 px, the FLD detected significantly more lines that, however, show a greater variance in the orientation. In this plot, only minimal differences can be detected visually and quantitatively for dt > 2 1/2 px. Results for dt > 10 px show no further changes in line detection both visually and quantitatively. Therefore, a value of dt = 2 1/2 px provides low SDs even at low values for t. However, a strongly decreasing number of lines has to be considered in this case. Thus, it occurs with other images from the dataset with application of the FLD that with a very small number of lines already, small deviations in the orientation of the single lines are responsible for an increased SD. This results in an unclear or only impaired determination of the main cultivation direction. In addition, the FLD no longer detects any lines if the values for t are too high in some cases, depending on the parcel area. Therefore, a relatively small value for dt should be chosen (dt ≤ 5.0 px) in order to connect continuous objects. The minimum was set as the distance from a neighboring pixel (dt = 2 1/2 px), which is the standard value for dt in the FLD.
The parameters cth1 and cth2 describe the upper and lower threshold values, respectively, in the final step of Canny edge detection (see Section 2.3.3). The Canny Threshold parameters have a significantly less sensitive effect on the results than the parameters t and dt. Therefore, we used statistical parameters instead of image results in Table 4 to present the results from investigating these parameters. Upper and lower values were automatically assigned so that cth1 and cth2 are interchangeable without changing the results of the line detection algorithm. Canny [60] recommends a threshold ratio for cth2 and cth1 between 1:2 and 1:3. In the following, cth1 describes the upper threshold, and for cth2, we multiply cth1 by the corresponding ratio. For this representation, we used Equation (5) for calculating the parameter t in combination with dt = 2 1/2 px. For one image of the category spring barley, no tramlines were visually recognizable. The quality of the results for this image was rather random regardless of the parameters and therefore not used for further adjustment. Table 4 shows the influence of cth1 and cth2 on the results of FLD in terms of number of lines as well as SD of detected line orientation at the different ratios of 1:4, 1:3, 2:5, 1:2, and 2:3. Values of 25, 50, and 100 for cth1 did not result in any changes for all ratios regarding number of lines and SD for images of spring barley compared with the results for cth1 = 150. Therefore, Table 4 only contains the results for cth1 ≥ 150 with a stepwise increase up to a maximum of cth1 = 300. With the increasing value for the upper threshold cth1, the mean number of recognized lines per image, and therefore the sum of lines, decreased, while the mean of the SD for line orientation decreased slightly. The mean line length did not show any changes. The influence of different ratios on the results for the crop type spring barley only became apparent for cth1 ≥ 200 in combination with high ratios of cth2:cth1 > 1:3. Considering all images simultaneously, better results became apparent for cth1 = 300, since the SD for the lines was lowest here. At the same time, however, the number of detected lines became smaller. Considering the single images separately, different results were obtained. In most cases, the number of detected lines changed, mostly in the form of a decrease. At the same time, the SD changed, yet with different effects. In some cases, the changes were very small (<2 • ), while in other cases there were larger positive or even negative effects on the SD. Table 5 shows the best results of the FLD for different Canny Threshold parameter constellations considering all seven remaining images for spring barley. It contains the absolute number of images with successful detection of the main cultivation direction by FLD (MAE ≤ 5 • ) compared with the visual detection of Google Earth™ images. Furthermore, the numbers of images with unsuccessful detection and with no line detection are listed. Table 5 shows similar results to Table 4. The FLD successfully determined the main cultivation direction for six of the seven images at cth1 ≤ 200 regardless of the ratio. The FLD showed the best result for cth1 = 300 and a ratio of 2:5. In none of the seven fields did the FLD detect no lines at all. However, no optimal parameter constellation for the Canny Threshold parameters can be derived from the results.

Application of the Adjusted Parameters to the Image Data Set
For the spring barley image dataset, the parameters t and dt were identified as very sensitive. Cth1 and cth2 showed less sensitive effects on the FLD results. Therefore, Table 6 presents the best results for different values for cth1 and ratios of the Canny Threshold parameters to account for their impact on the entire dataset. In this case, we used all 278 images even if no tramlines were visually detectable. The table contains the absolute number of images with successful detection of the main cultivation direction by FLD (MAE ≤ 5 • ). Furthermore, we list the number of images with unsuccessful detection and with no line detection. Even though the best results were obtained for different upper thresholds in combination with different ratios, the rate of successfully detected main cultivation direction was the same for all applied thresholds, with a success rate of 77.7%. In addition, it is obvious that the FLD did not detect any lines at all for significantly more images at an upper threshold of cth1 = 300, which were only unsuccessfully detected at lower thresholds. Changing the Distance Threshold to dt = 3 px, the results became  Selecting the optimal Canny Threshold parameters cth1 and cth2 for the respective crop type resulted in a success rate of 79.9%. Table 7 contains the individual results for this application. However, when analyzing these results, the very small number of image acquisitions for some crop types must be taken into account, which means that the results are not representative for all crop types. The FLD did not detect lines for one image each of winter spelt, winter barley, potatoes, and sugar beet. Furthermore, detection of the main cultivation direction was less successful for maize (non-silage) and winter rapeseed. Considering the crop categories (see Table 2), different trends can be observed for the upper threshold value cth1. For root crops, summer cereals, and other crops, the FLD provided better results for smaller values of cth1 from 25 to 150, while for oil crops and winter cereals, it achieved better results for a large value of cth1 = 300. In this study, the value for cth1 showed little to no influence on horticultural crops and plants for green harvesting. Selecting the optimal Canny Threshold parameters cth1 and cth2 for the respective crop type resulted in a success rate of 79.9%. Table 7 contains the individual results for this application. However, when analyzing these results, the very small number of image acquisitions for some crop types must be taken into account, which means that the results are not representative for all crop types. The FLD did not detect lines for one image each of winter spelt, winter barley, potatoes, and sugar beet. Furthermore, detection of the main cultivation direction was less successful for maize (non-silage) and winter rapeseed. Considering the crop categories (see Table 2), different trends can be observed for the upper threshold value cth1. For root crops, summer cereals, and other crops, the FLD provided better results for smaller values of cth1 from 25 to 150, while for oil crops and winter cereals, it achieved better results for a large value of cth1 = 300. In this study, the value for cth1 showed little to no influence on horticultural crops and plants for green harvesting. Selecting the optimal Canny Threshold parameters cth1 and cth2 for the respective crop type resulted in a success rate of 79.9%. Table 7 contains the individual results for this application. However, when analyzing these results, the very small number of image acquisitions for some crop types must be taken into account, which means that the results are not representative for all crop types. The FLD did not detect lines for one image each of winter spelt, winter barley, potatoes, and sugar beet. Furthermore, detection of the main cultivation direction was less successful for maize (non-silage) and winter rapeseed. Considering the crop categories (see Table 2), different trends can be observed for the upper threshold value cth1. For root crops, summer cereals, and other crops, the FLD provided better results for smaller values of cth1 from 25 to 150, while for oil crops and winter cereals, it achieved better results for a large value of cth1 = 300. In this study, the value for cth1 showed little to no influence on horticultural crops and plants for green harvesting. Selecting the optimal Canny Threshold parameters cth1 and cth2 for the respective crop type resulted in a success rate of 79.9%. Table 7 contains the individual results for this application. However, when analyzing these results, the very small number of image acquisitions for some crop types must be taken into account, which means that the results are not representative for all crop types. The FLD did not detect lines for one image each of winter spelt, winter barley, potatoes, and sugar beet. Furthermore, detection of the main cultivation direction was less successful for maize (non-silage) and winter rapeseed. Considering the crop categories (see Table 2), different trends can be observed for the upper threshold value cth1. For root crops, summer cereals, and other crops, the FLD provided better results for smaller values of cth1 from 25 to 150, while for oil crops and winter cereals, it achieved better results for a large value of cth1 = 300. In this study, the value for cth1 showed little to no influence on horticultural crops and plants for green harvesting. Selecting the optimal Canny Threshold parameters cth1 and cth2 for the respective crop type resulted in a success rate of 79.9%. Table 7 contains the individual results for this application. However, when analyzing these results, the very small number of image acquisitions for some crop types must be taken into account, which means that the results are not representative for all crop types. The FLD did not detect lines for one image each of winter spelt, winter barley, potatoes, and sugar beet. Furthermore, detection of the main cultivation direction was less successful for maize (non-silage) and winter rapeseed. Considering the crop categories (see Table 2), different trends can be observed for the upper threshold value cth1. For root crops, summer cereals, and other crops, the FLD provided better results for smaller values of cth1 from 25 to 150, while for oil crops and winter cereals, it achieved better results for a large value of cth1 = 300. In this study, the value for cth1 showed little to no influence on horticultural crops and plants for green harvesting.  Figure 18 shows the evaluation of the results taking into account the parcel geometry with the parameters parcel area (a) and the compactness factor CIPQ (b). For a comparative presentation, both parameters were classified and the absolute values for successful and unsuccessful line detection are listed as well as the number of images for which the FLD did not detect any lines at all. With increasing area, the relative amount of successfully 0% 50% 100%  Figure 18 shows the evaluation of the results taking into account the parcel geometry with the parameters parcel area (a) and the compactness factor C IPQ (b). For a comparative presentation, both parameters were classified and the absolute values for successful and unsuccessful line detection are listed as well as the number of images for which the FLD did not detect any lines at all. With increasing area, the relative amount of successfully detected main cultivation directions increased from 75 to 100%. The three images for which the FLD did not detect any lines were within the classes with a smaller parcel area ≤20,000 m 2 . Considering the C IPQ , the success rate ranged from 68.8% for very large C IPQ values to 100% for values smaller than 0.3. However, the success rate was apparently independent from the compactness factor, as no functional relationship between successful detection of the main cultivation direction and the C IPQ could be established. detected main cultivation directions increased from 75 to 100%. The three images for which the FLD did not detect any lines were within the classes with a smaller parcel area ≤20,000 m 2 . Considering the CIPQ, the success rate ranged from 68.8% for very large CIPQ values to 100% for values smaller than 0.3. However, the success rate was apparently independent from the compactness factor, as no functional relationship between successful detection of the main cultivation direction and the CIPQ could be established.
FLD Successful FLD Not Successful FLD No Line Detected (a) (b) Figure 18. Results of FLD with parameters of (a) the shape area (different class widths) and (b) the compactness factor CIPQ; percentages describe the success rate of the FLD per interval.
Finally, we present the results of the comparison of mean aspect and main cultivation direction in Table 8 for all 278 arable land parcels. From the smallest difference of these two directions, we categorized the support practice measure contouring into three categories. Optimal contouring occurs when the two lines are orthogonal and the difference is in the interval [67.5°, 90°]. The category of no contouring exists for a cultivation along the direction of the aspect if the directions are almost parallel and the difference is in the interval of [0°, 22.5°). For values of the difference that lie in the interval [22.5°, 67.5°), we define an oblique cultivation. From this comparison, we deduced for the individual arable land parcels that the support practice measure contouring is used for 149 parcels (53.6%) within the Fischbach catchment. For 33 (11.9%) of the analyzed arable land parcels, the mean aspect corresponds to the main cultivation direction (no contouring). The conclusion for the FLD is comparable, only rather worse. At 50.0%, slightly fewer arable land parcels are categorized as having contouring, while 13.7% belong to the category no contouring. The value for oblique cultivation is almost the same at 35.3%. The FLD was not able to determine a cultivation direction at all for three images. Table 8. Final results for the application of the support practice measure contouring comparing the mean aspect vs. the main cultivation direction. Example illustration on the left: Representation of the classification for contouring with mean aspect (black arrow) and cultivation direction manually derived (blue) and automatically derived (green).  Finally, we present the results of the comparison of mean aspect and main cultivation direction in Table 8 for all 278 arable land parcels. From the smallest difference of these two directions, we categorized the support practice measure contouring into three categories. Optimal contouring occurs when the two lines are orthogonal and the difference is in the interval [67.5 • , 90 • ]. The category of no contouring exists for a cultivation along the direction of the aspect if the directions are almost parallel and the difference is in the interval of [0 • , 22.5 • ). For values of the difference that lie in the interval [22.5 • , 67.5 • ), we define an oblique cultivation. From this comparison, we deduced for the individual arable land parcels that the support practice measure contouring is used for 149 parcels (53.6%) within the Fischbach catchment. For 33 (11.9%) of the analyzed arable land parcels, the mean aspect corresponds to the main cultivation direction (no contouring). The conclusion for the FLD is comparable, only rather worse. At 50.0%, slightly fewer arable land parcels are categorized as having contouring, while 13.7% belong to the category no contouring. The value for oblique cultivation is almost the same at 35.3%. The FLD was not able to determine a cultivation direction at all for three images. Table 8. Final results for the application of the support practice measure contouring comparing the mean aspect vs. the main cultivation direction. Example illustration on the left: Representation of the classification for contouring with mean aspect (black arrow) and cultivation direction manually derived (blue) and automatically derived (green). detected main cultivation directions increased from 75 to 100%. The three images for which the FLD did not detect any lines were within the classes with a smaller parcel area ≤20,000 m 2 . Considering the CIPQ, the success rate ranged from 68.8% for very large CIPQ values to 100% for values smaller than 0.3. However, the success rate was apparently independent from the compactness factor, as no functional relationship between successful detection of the main cultivation direction and the CIPQ could be established.

Absolute Difference between Mean Aspect and Main Cultivation Direction
FLD Successful FLD Not Successful FLD No Line Detected (a) (b) Figure 18. Results of FLD with parameters of (a) the shape area (different class widths) and (b) the compactness factor CIPQ; percentages describe the success rate of the FLD per interval.
Finally, we present the results of the comparison of mean aspect and main cultivation direction in Table 8 for all 278 arable land parcels. From the smallest difference of these two directions, we categorized the support practice measure contouring into three categories. Optimal contouring occurs when the two lines are orthogonal and the difference is in the interval [67.5°, 90°]. The category of no contouring exists for a cultivation along the direction of the aspect if the directions are almost parallel and the difference is in the interval of [0°, 22.5°). For values of the difference that lie in the interval [22.5°, 67.5°), we define an oblique cultivation. From this comparison, we deduced for the individual arable land parcels that the support practice measure contouring is used for 149 parcels (53.6%) within the Fischbach catchment. For 33 (11.9%) of the analyzed arable land parcels, the mean aspect corresponds to the main cultivation direction (no contouring). The conclusion for the FLD is comparable, only rather worse. At 50.0%, slightly fewer arable land parcels are categorized as having contouring, while 13.7% belong to the category no contouring. The value for oblique cultivation is almost the same at 35.3%. The FLD was not able to determine a cultivation direction at all for three images. Table 8. Final results for the application of the support practice measure contouring comparing the mean aspect vs. the main cultivation direction. Example illustration on the left: Representation of the classification for contouring with mean aspect (black arrow) and cultivation direction manually derived (blue) and automatically derived (green).

Land Use and Land Parcel Analysis
The analysis of land use and arable land parcels has shown that areas with an aboveaverage slope dominate the Fischbach catchment and are predominantly forested. The share of arable land is therefore lower than that in Hesse and Germany. In the category of agricultural land, less steeply sloping terrain offers more opportunity for arable land, while permanent grassland is more likely to be found in regions with steeper slopes. On the one hand, this is due to the risk of soil erosion on steeper slopes, but on the other hand, it is more difficult to machine steeper sloped plots [1]. This phenomenon is applicable for Hesse and especially for the investigated catchment area. The study area shows, in accordance with the statistics for Hesse and Germany, typical cultivation structures for arable land with a large share of cereals containing typical grain crops such as winter wheat and winter barley. These accounted for a share of over 55% of the cultivated area in 2019. The arable land parcels in the Fischbach catchment are, on average, smaller in relation to Hesse, which is due to the topography. The steeper and more complex terrain is responsible for the frequently alternating land uses in the form of forest and permanent grassland. Fault scarps in the terrain as well as additional landscape features such as hedges, hems, copses, or small water bodies limit the extent of the arable land parcels. This leads to smaller parcels on average. Smaller parcels in turn have a positive effect on the slope length, since this is automatically accompanied by a shortening of the maximum possible slope length. This in turn results in a shorter flow path and thus in a lower transport force of the surface runoff. However, there is still an increased risk of erosion in the catchment area due to an above-average slope and a high share (>70%) of erosion-prone loess soils such as brown earth and parabrown earth in the uppermost soil horizon with a high share of silt [63][64][65].
Depending on the geometric shape and the orientation of a parcel in relation to the terrain topography, the direction of cultivation must be considered from an economic point of view. With an unfavorable ratio of width to length of a field, machining is only possible in the direction of the longer side, as frequent turning maneuvers would make machining of the field uneconomical [43]. If narrow fields or fields with an unfavorable width-to-length ratio are mismanaged, then the likelihood of soil erosion from downslope tramlines is significantly increased. Other land uses, such as permanent grassland, or land consolidation may need to be considered. This would transpose the field boundaries to achieve a more favorable arrangement of individual croplands, taking into account the topography. Furthermore, a joint, cross-border field management within a field block is conceivable in cooperation with the field neighbor [1].
The compactness factor C IPQ was introduced with the aim of taking the geometric shape and the area of the individual land parcels into account. In addition, the factor should be used to examine the influence of the area and geometric shape during line detection. For this purpose, the C IPQ was chosen because it is easy to calculate using perimeter and area data, thus allowing quantification of the polygon shape. The factor could be applied to all studied geometric shapes of arable land parcels since there was no overlap of the outer edges of the polygons. The application worked well, with a few exceptions, yet was not always successful. It was found that the results do not allow an explicit conclusion on the simplified, geometric basic shape in every case. This is due to the topology of the fields, which, in practice, have significantly more nodes and edges than regular triangles or quadrilaterals. More nodes lead to a concave shape of land parcels. That is, they approximate a circular shape. Convex structures, on the other hand, reduce the C IPQ . Larger objects with irregular boundaries tend to have a higher compactness [66]. Although the C IPQ can clearly distinguish between areas with unfavorable aspect ratios and areas with favorable aspect ratios, there are some exceptions. Holes within polygons significantly reduce the C IPQ by increasing the perimeter. For a corrected calculation, only the outer perimeter has to be considered. Nevertheless, the automatic identification of these polygons with holes must be the first step. Therefore, the results of the C IPQ can only serve as a first approximation to estimate the geometric shape. In the future, it would be better to include other shape descriptors, such as the elongation shape factor or the aspect ratio [67,68]. In particular, this is to identify narrow fields with unfavorable aspect ratios in order to be able to check their orientation with respect to the topography. If necessary, a combination of several factors can be used to make an appropriate classification of the geometric shapes and to check for a correlation between the detection of the main cultivation direction by the FLD and the geometric shape [69]. The results should be applied to a study area with larger parcels to cover the full range of possible values for the area of land parcels and not only rather small parcels as used in this study. This would allow us to test the hypothesis of whether C IPQ actually tends to be higher for larger plots.

DEM Analysis
The comparison of the statistical parameters for slope and aspect has shown that, with respect to the entire catchment, there are significant differences in the calculated values for the different datasets DEM1 and DEM5. These result from the different accuracies of the derived terrain models. If the spatial extent is restricted to the arable land parcel, the differences between the results of DEM1 and DEM5 are reduced. This applies to both parameters, slope and aspect. Among other things, this is due to the fact that areas with very steep slopes are unsuitable for agricultural cultivation and therefore tend to be forested or covered with permanent grassland in the Fischbach catchment (see Section 4.1). However, with the choice of an arable land parcel as the reference area or spatial unit for the calculation of the P-factor, the determination of an average value for this area unit is sufficient. In this case, the deviations between the different datasets are much smaller and therefore negligible. The correlation for mean slope and mean aspect of the two DEMs is very strong, although it depends on the defined threshold for a deviation. The parameter SD shows the same behavior. It can be stated that a DEM with a resolution of 5 m × 5 m is completely sufficient for terrain analysis according to this method in the Fischbach catchment. However, this decision is individual for each study area and is based on the topography of the catchment. At least the SD should be used as an indicator for the dispersion of these parameters. Especially for calculation of aspect, the applied method has to be questioned critically. In areas with medium slopes and complex topography, there may be different aspects within one land parcel (see Figure 13). In this case, the SD only serves as an indicator for the uniqueness of the aspect. Figure 13 shows that the complexity of the topography in relation to the SD of aspect decreases with an increasing slope. Nevertheless, this complexity has to be considered especially for steep slopes during an evaluation. The alternative is to use the individual contour lines within the land parcel to compare the main cultivation direction with the local aspect [70].
For future applications, a split of the land parcels depending on different aspects would be conceivable. In terms of mean slope, the Fischbach catchment is ideally suited as a test area for the application of the presented method since contouring is recommended as effective for slopes between 3 and 15% [1,14]. Overall, 84.7% of the arable land parcels have slopes between 3 and 15%, with a mean slope of 10.5% for all arable land parcels. Nevertheless, there are some arable land parcels (14.7%) with increased slope values greater than 15%.

Image Analysis
The method presented in this study is a simplified application of line detection on arable land parcels to automatically derive a main cultivation direction from the image data. The main cultivation direction was determined by an average orientation of the lines, weighted by their length. The aim of line detection is to detect a certain number of lines with a certain minimum length in order to be able to make a valid statement about the orientation of the lines. Only if the FLD algorithm detects a minimum number of lines does the standard deviation become meaningful. With a smaller number of lines, even a few outliers ensure that the determined orientation of the lines has a high variability and is therefore not reliable. In addition, the lines should not be too short since tramlines usually run across the entire field without interruption. However, since the FLD can only detect linear elements, it was necessary to choose a threshold value for the line length depending on the field size. The adjustment in Section 3.3.1 illustrates that the parameters Length Threshold t and Distance Threshold dt have a very sensitive effect on the result of the line detection. When varying these parameters, the FLD showed significant changes in the results with respect to the number of lines, average line length, and the SD of the alignment of the lines. For the parameter Length Threshold, a function for t depending on the field size could be successfully defined with Equation (5). It was taken into account to define an upper and a lower limit for the line length. This ensured a minimum line length of 6 m for small fields. The required line length must not exceed the field length, but only represent a certain percentage of it. On larger fields, the required minimum length is increased depending on the area. However, the required length should not be too long. It was limited by an upper threshold value of 30 m so that tramlines with a slight curve could be detected (see Section 3.3.1). In addition, variations in the gradient of the gray levels can occur, resulting in an interruption or splitting of the line elements. If the threshold value is too high, both linear elements would then be evaluated as too short and not detected. Furthermore, we found a suitable default value for the parameter Distance Threshold dt for the application on arable land parcel images with the FLD of dt = 2 1/2 px. The Distance Threshold shows effects on the result in interaction with t, yet is a little less sensitive than t. Furthermore, there is a dependency between the parameters dt and the Canny Aperture Size cas, yet this was not investigated further in this study.
The upper threshold for Canny edge detection cth1 and, depending on the ratio, the lower threshold cth2 are less sensitive than t and dt. Different ratios were chosen based on Canny's recommendation [60]. However, no clear values could be derived from the adjustment and from the application on the entire image dataset. As a conclusion, it can be stated that the combination of a high, upper threshold cth1 and large ratios (1:2 and 2:3) leading to a high cth2 value often reacted with a higher sensitivity than other parameter constellations. However, the result for each image was often individual, so the FLD results were both better and worse in terms of the number of lines and SD for alignment. The final recommendation for the Canny Threshold parameter is to choose a ratio of 1:2 or smaller (2:5, 1:3, or 1:4) to avoid negative effects on line detection. For the application of the adjusted parameters to the entire image dataset in Section 3.2.2, we set the maximum absolute deviation MAE to 5 • so that the automatic derivation of the main cultivation direction is evaluated as successful. This included both a positive and a negative deviation of 5 • with respect to the manually determined main cultivation direction from the Google Earth™ image dataset. We set this limit of 5 • because the measured deviation in the manual detection was approximately ±2 • . In addition, the deviation should not be too large as there is a risk of lateral runoff and thus linear erosion if the field is worked oblique to the contour lines [1]. Figure 19 presents the influence of the selection of the limit value MAE for a successful line detection by the FLD. Figure 19a shows the results for different upper thresholds cth1 with a ratio of 1:2 and Figure 19b plots the results for an upper threshold of cth1 = 300 with different ratios. The x-axis shows the absolute limits for MAE, while the y-axis contains the shares of a successful detection. Both graphs clearly show an increase in the FLD success rate with increasing MAE limit value. This increase is stronger for smaller limit values. However, it decreases with increasing MAE values and corresponds in its form to a logarithmic relationship. In addition, Figure 19a shows that no clear statement can be made about an optimal value for the parameter cth1 since the success rate for different absolute MAE thresholds provides an optimal result for different upper thresholds in each case. Figure 19b shows a similar result for the lower cth1 threshold. This plot for cth1 = 300 contains the largest variation in the success rate for the different ratios among all the values examined for cth1. Once again, no optimal ratio was evaluated as a function of MAE. However, the detection success rate was already above 50% for a MAE of 2 • , which corresponds to the estimated measurement error. ent ratios among all the values examined for cth1. Once again, no optimal ratio was evaluated as a function of MAE. However, the detection success rate was already above 50% for a MAE of 2°, which corresponds to the estimated measurement error. Crops of different categories differ significantly in their morphology, their seasonally recurring development, and, thus, their phenological properties. This results in some differences in the optimal values of the Canny Threshold parameters that, related to the small size of the dataset, indicate a first tendency for the crop categories. Cereals show the clearest formation of seed rows, including small row spacing and very little individual grain spacing within the seed row. Only a few leaves grow directly on the woody stem. Herbaceous plants such as rapeseed grow without lignification of the stem [71]. Leaf and stem parts are rather unsorted, with multiple junctions on which small leaves grow. From a plan view, these plants tend to show patchy, gapless growth between the seed rows. Tramlines are visible depending on the stage of development. Due to the single-grain sowing method for root crops and the large grain spacing within the seed row, the seed rows are not yet linear elements and the tramlines are thus impaired [72]. In addition, the headland for maize and other root crops varies in size depending on the polygon shape. However, it is in some cases significantly larger than for other crop types. These morphological influences are evident both in terms of crop categories and the time of sowing. For summer cereals, root crops, and others, which are all sown in spring, growth is not yet sufficiently advanced [53,72]. These crops show the best result for low cth1 values between 25 and 150 (see Section 3.3.2). Winter cereals and oil crops are more advanced in their growth towards the end of June. They show an optimal result for a large threshold of cth1 = 300. However, seed row detection is limited for rapeseed and all root crops because of the non-linear shape in the top view [53,72]. The recording time of the images seems to be well suited on average, although possibly a little later, around mid-July, is optimal for this catchment. However, calibration of the individual crop types does not seem to be possible with the available parameters. The result is only minimally better for an individual adjustment for each crop type. However, the number of images for the individual crop types is too small for a general statement. Only winter wheat and winter barley allow a meaningful result based on the large number of images. The success rate of the FLD for these crop types was over 75%.  Crops of different categories differ significantly in their morphology, their seasonally recurring development, and, thus, their phenological properties. This results in some differences in the optimal values of the Canny Threshold parameters that, related to the small size of the dataset, indicate a first tendency for the crop categories. Cereals show the clearest formation of seed rows, including small row spacing and very little individual grain spacing within the seed row. Only a few leaves grow directly on the woody stem. Herbaceous plants such as rapeseed grow without lignification of the stem [71]. Leaf and stem parts are rather unsorted, with multiple junctions on which small leaves grow. From a plan view, these plants tend to show patchy, gapless growth between the seed rows. Tramlines are visible depending on the stage of development. Due to the single-grain sowing method for root crops and the large grain spacing within the seed row, the seed rows are not yet linear elements and the tramlines are thus impaired [72]. In addition, the headland for maize and other root crops varies in size depending on the polygon shape. However, it is in some cases significantly larger than for other crop types. These morphological influences are evident both in terms of crop categories and the time of sowing. For summer cereals, root crops, and others, which are all sown in spring, growth is not yet sufficiently advanced [53,72]. These crops show the best result for low cth1 values between 25 and 150 (see Section 3.3.2). Winter cereals and oil crops are more advanced in their growth towards the end of June. They show an optimal result for a large threshold of cth1 = 300. However, seed row detection is limited for rapeseed and all root crops because of the non-linear shape in the top view [53,72]. The recording time of the images seems to be well suited on average, although possibly a little later, around mid-July, is optimal for this catchment. However, calibration of the individual crop types does not seem to be possible with the available parameters. The result is only minimally better for an individual adjustment for each crop type. However, the number of images for the individual crop types is too small for a general statement. Only winter wheat and winter barley allow a meaningful result based on the large number of images. The success rate of the FLD for these crop types was over 75%. The results of the FLD for maize (non-silage), potatoes, winter spelt, and winter rapeseed were below average poor. These findings need further investigation based on a larger population.
No clear influence on the result was found for the compactness factor C IPQ . Perhaps the application of an alternative geometric shape factor can produce a correlation between FLD success rate and the polygon shape. In terms of land parcel size, a higher success rate was obtained especially for parcels ≥20,000 m 2 . However, the distribution of land parcels for the defined classes in Figure 18a with a positive skewness is very unequal. Therefore, it is necessary to conduct an independent study with a much larger population for an investigation of the influence of area size on the FLD success rate.
Sources of error in the context of image analysis are partly due to the crop type. However, a major part is caused by the individual characteristics of the images. We identified the following causes for unsuccessful line detection during the analysis:  Table 9 shows the causes for unsuccessful line detection differentiated for the categories 'Unsuccessful'-for example, line detection with a deviating main cultivation direction-and 'No lines detected'-for example, no detection of a line at all. Considering that Table 9 sometimes lists multiple reasons for one image, the main problems for determining a deviating main cultivation direction are identifying headland, existing soil erosion, and widely varying albedo as well as outliers. Figure 20 illustrates examples of the abovementioned causes of line detection errors for the different categories. We discuss the individual categories separately below. potatoes, winter spelt, and winter rapeseed were below average poor. These findings need further investigation based on a larger population. No clear influence on the result was found for the compactness factor CIPQ. Perhaps the application of an alternative geometric shape factor can produce a correlation between FLD success rate and the polygon shape. In terms of land parcel size, a higher success rate was obtained especially for parcels ≥20,000 m 2 . However, the distribution of land parcels for the defined classes in Figure 18a with a positive skewness is very unequal. Therefore, it is necessary to conduct an independent study with a much larger population for an investigation of the influence of area size on the FLD success rate.
Sources of error in the context of image analysis are partly due to the crop type. However, a major part is caused by the individual characteristics of the images. We identified the following causes for unsuccessful line detection during the analysis:  Table 9 shows the causes for unsuccessful line detection differentiated for the categories 'Unsuccessful'-for example, line detection with a deviating main cultivation direction-and 'No lines detected'-for example, no detection of a line at all. Considering that Table 9 sometimes lists multiple reasons for one image, the main problems for determining a deviating main cultivation direction are identifying headland, existing soil erosion, and widely varying albedo as well as outliers. Figure 20 illustrates examples of the abovementioned causes of line detection errors for the different categories. We discuss the individual categories separately below.

Soil Erosion/Albedo
Existing soil erosion as well as varying albedo of soils, often caused by a spatially differentiated growth stage within a land parcel, pose a major problem in line detection (see Figure 20a). In this case, irregularities in the grayscale of the image are responsible for the lines being interrupted at different points. In addition, the FLD detects both linear erosion structures and the outer edges of sheet erosion as well as the deposits. The orientation of the detected lines of the erosion structures significantly influences the result of the main cultivation direction. If the growth stage of the plants is sufficiently advanced, they will quickly cover the erosion structures present at the soil surface. However, soil erosion often has a negative impact on the structure of plant growth and is still clearly visible at later growth stages. These problems caused the FLD to be unsuccessful in detection for about one-third of the images (see Table 9). Erosion is present regardless of geometry or plant type.

No Tramlines/Polygon Shape
The image in Figure 20b shows no tramlines. Seed rows or tramlines must be present and identifiable for successful detection of the main cultivation direction. For all three images for which the FLD could not detect any lines (100%), no tramlines were visually detectable. All in all, the detection was unsuccessful for 11.9% of all images due to missing tramlines (see Table 9). The reason for the lack of visually detectable tramlines or seed rows may be the early stage of the plant growth for certain crop types [53]. This means that the plants have just been sown and no fertilization or other mechanical treatment has taken place yet, so no tramlines are present. In this case, plant growth is usually not advanced and the FLD is not able to detect lines on the remote sensing images. Furthermore, identification of lines is not possible with herbaceous, ground-covering vegetation such as grass. No-till farming or undersown crops can have an influence on the appearance of tramlines and seed rows. Here, the timing of the image recording in relation to the developmental progress of the plants is decisive, too. Furthermore, the poor line detection can be attributed to a too-narrow field. The reason for unsuccessful detection was due to a too-narrow polygon shape in 3.4% of the images. In some cases, a very complex shape with convex formations was the source of error. Furthermore, the FLD did not detect any lines at all for a total of two images because of the polygon shape.

Landscape Features/Headland
Arable land parcels usually do not contain landscape features [39]. However, due to inaccuracies in digitizing and defining these areas, elements such as copses, shrubs, or forests may extend into the field (see Figure 20c). In addition, depending on the position in the sky, the sun creates shadows on the ground due to these objects when recording these images. In this study, the application of an additional buffer prior to application of the FLD proved to be effective (see Section 2.3.3). The buffer size depends on the field area such that a large part of the headland as well as landscape features and their cast shadows at the field edges are not considered in the line detection. However, for some parcels, the buffer was not large enough, so partial headlands were considered during line detection (see Figure 20d). This was the main reason for unsuccessful detection. In only one case landscape features still had a disturbing effect on the detection. However, in 44.1% of the unsuccessful cases, headland was responsible for incorrect line detection (see Table 9).

Curves
The FLD algorithm can only detect linear elements. Therefore, we removed all parcels that had multiple cultivation directions from the analysis before applying the FLD. On a few remaining parcels, slightly curved tracks were visible because the farmer adapted the cultivation to the field boundaries. Since the FLD does not detect the existing lines to the same extent on the whole field in Figure 20e, the final weighted mean of the direction does not correspond to the main cultivation direction of the field. Tramlines with curved elements were responsible for 15.3% of all unsuccessful detections (see Table 9).

Outlier
In combination with existing soil erosion, landscape features, or albedo gradients, the FLD detected individual lines in some images that deviated significantly from the main cultivation direction (see Figure 20f). With a sufficiently large number of detected lines, individual outliers do not have such a strong effect on the weighted result. However, with fewer lines and particularly strong deviation, the influence of one single outlier line is so strong that a successful determination of the main cultivation direction fails. Outliers caused unsuccessful detection for 30.5% of all images (see Table 9).
Most of the problems encountered in applying the FLD to arable land parcels can be solved using improvements within the algorithm or previously in data preparation. By applying a mask with an adjusted buffer radius based on a new shape factor and field area, the landscape features and headland at the edge of the field can be excluded from the line detection. With a sufficient number of detected lines per parcel, individual outliers can be identified as significant and disregarded in the subsequent calculation of the mean line orientation. For this application, however, the number of detected lines must be higher for some parcels for which the FLD has only detected a very small number of lines so far. The consideration of curve elements for seed rows and tramlines is more difficult since they require either automatic splitting of the parcel into different sections based on parameters to be defined or another detection algorithm for curve elements [73]. Improvement within the FLD algorithm or calculating the mean line orientation is not possible for problems of soil erosion, albedo, and non-detectable tramlines. We think that an improvement of the detection for these categories can be achieved by applying filter methods, such as sharpening [50]. If available, higher-resolution images can be used in the analysis. In this study, we did not choose the best possible resolution for the data export from Google Earth™ for the image analysis. Although a better resolution is available for the study area and the time of the analysis, this results in a larger dataset and thus more time-consuming export and data preparation. In addition, the availability of data resolution varies by region with an increasing trend towards better resolution [37]. However, due to data processing, the images are available with a time delay. In recent years, the availability of images for this study area has increased significantly: for the period 2015 to 2020, one image is available for every year except 2017, although the month of recording varies. The latest image of 2020 has a maximum resolution of about 0.2 m × 0.2 m per pixel [41]. However, as the size of the study area increases, the probability increases that different recording times and thus different light conditions or even resolutions are present in Google Earth™. When considering the available resolution, it is necessary to take a possible loss of quality during data processing into account. During data export and import as well as associated georeferencing and image mosaicking, there may be a noticeable loss of data quality. Image distortion during rectification or even a reduction in the original resolution due to compression should be avoided at all costs [50]. In further studies, the influence of the image resolution should be investigated stepwise in order to define a minimum resolution for the application with an adequate result. Furthermore, the time of recording is an interesting point of investigation for follow-up studies with regard to the plant growth stages.
With an improved method of calculating the buffer and successful detection of the outliers, it is estimated that 35 to 50% of previously unsuccessful detections could be improved in the adapted application of the method. This could increase the detection success rate up to 88%.

Conclusions and Outlook
The objective of this study was to develop an automated method for determining the main cultivation direction of agricultural land parcels with seasonal crops. To detect the main cultivation direction on each land parcel, the Fast Line Detector (FLD) for detecting linear elements was applied to open remote sensing data retrieved from Google Earth™. The method was tested for the low mountain range catchment of the Fischbach, Germany.
The mean slope and mean aspect from the digital elevation models DEM1 and DEM5 were derived for each parcel. Subsequently, the main cultivation direction was compared to the mean aspect with the objective of obtaining information about the soil conservation practice of contouring for each arable land parcel. The results were able to provide a helpful tool for soil loss estimation and for the determination of the P-factor within the USLE.
The research area, namely the Fischbach catchment, located in the low mountain range with typical medium slopes and erosion-prone loess soils was proven well suitable for the application of contouring and the potential of applying the P-factor. The crop types in the study area are typical for Hesse. They encompass a high proportion of cereals, although the parcel sizes are rather small compared to other patches in Hesse and Germany. The differences in the results for DEM1 and DEM5 for the mean slope and mean aspect were negligible given that this study utilized the land parcel as a spatial unit. Thus, DEM5 was shown to be suitable for this application.
The initial application of the FLD for the automated derivation of the main cultivation direction achieved appropriate results with a success rate of 77.7% for a uniform specification of parameter settings for all images of 278 arable land parcels. The FLD parameters Length Threshold t and Distance Threshold dt were particularly sensitive. In addition, the FLD parameters for the Canny Threshold indicated tendencies with respect to the appearance of unique seed rows and tramlines in the remote sensing images depending on the time of sowing or the plant growth stage in relation to the time of image recording. However, the application of the FLD with crop-type-specific parameters yielded only a slightly higher success rate of 79.9% and, therefore, did not result in any significant improvements.
The main reasons for the determination of a deviating main cultivation direction within the FLD are due to problems with detected headland, existing soil erosion, and strongly varying albedo within the parcels as well as single outliers in line detection. The application of a corrected buffer to the remote sensing images with enhanced parameterization offers promising improvements for a higher success rate of the FLD.
For further development of this new method, it is necessary to verify the results on a more extensive dataset with further crop types and larger land parcels. The requirements for a possible application of the FLD algorithm include available information on the topography (DEM5) and the geometry of the arable land parcels, including their associated crops, if available. In addition, images for line detection should be available as remote sensing data with a sub-meter resolution, preferably taken in the same year as the shape data of the plots. Furthermore, the influence of the image resolution and the time of image recording on the result of the line detection should be examined more closely. Finally, calculation of the annual soil loss estimation considering the P-factor would show its importance for the application of the USLE. Using this method, we identified a large number of fields that apply contouring. Therefore, this method contributes to a realistic determination of the factor for support measures and to better estimates of soil conservation by using real data. In addition, governmental institutions could use the presented method to identify land parcels that are not managed optimally in terms of contouring. This may allow the institutions to suggest appropriate measures to farmers, such as land consolidation or joint, cross-border field management for the affected land parcel.
This study takes a first step towards deriving the measure of contouring in an automated way based on open data, making it easier to anticipate the P-factor in soil erosion studies.