Scale Optimization in Topographic and Hydrographic Feature Mapping Using Fractal Analysis

: A new method for selecting optimal scales when mapping topographic or hydrographic features is introduced. The method employs rank-size partition of heavy-tailed distributions to detect nodes of rescaling invariance in the underlying hierarchy of the dataset. These nodes, known as head / tail breaks, can be used to indicate optimal scales. Then, the Fractal Net Evolution Assessment (FNEA) segmentation algorithm is applied with the topographic or hydrographic surfaces to produce optimally scaled objects. A topological transformation allows linking the two processes (partition and segmentation), while fractal dimension of the rescaling process is employed as an optimality metric. The new method is experimented with the two biggest river basins in Greece, namely Pinios and Acheloos river basins, using a digital elevation model as the only input dataset. The method proved successful in identifying a set of optimal scales for mapping elevation, slope, and ﬂow accumulation. Deviation from the ideal conditions for implementing head / tail breaks are discussed. Implementation of the method requires an object-based analysis program and few common geospatial functions embedded in most GIS programs. The new method will assist in revealing underlying environmental processes in a variety of earth science ﬁelds and, thus, assist in land management decision-making and mapping generalization.


Introduction
Most landscape features or processes, either natural or anthropogenic (e.g., a forest, a rail network, an urban area, a flood, or erosion), are scale-dependent [1]. Scale dependency is due to the fact that every feature or process has its own functional scale, i.e., a specific spatiotemporal interval within which it integrates or smooths to provide a message. The functional scale, in turn, will determine an appropriate mapping scale, i.e., the spatiotemporal unit required to visualize these properties for map readers.
Landscapes are composed of loosely coupled levels which function at distinct spatial and temporal scales and interact as nodes of an inherent hierarchy. In order to interpret possible interactions between these levels, it is necessary to decompose them into their fundamental parts [2].
Scale dependency and necessity of a multi-scale hierarchical approach has been indicated as an important issue in hydrological research for a long time now. Selection of an appropriate scale when mapping topographic or hydrographic features has been one of the trickiest problems, especially in land use and management projects [3].
The major concern about scale is how it may affect geospatial data collection and analysis results with respect to accuracy and reliability [1]. The scale problem has also been expressed by Zhang et al. (2018) [4] through the question: how many scales are necessary (authors' note: to describe a specific intrinsic parameter in fractal geometry is fractal dimension, defined as the ratio of the logarithmic change in detail (n) to the logarithmic change in scale (s) of geographical features [17]: From Equation (1), it is derived that fractal dimension may take decimal values, for example 1.26 or 2.65.
Scale in fractal geometry is primarily defined in a manner in which a series of scales are related to each other in a scaling hierarchy. For example, a coastline is a set of recursively defined bends, forming the scaling hierarchy of far more small bends than large ones. Fractal analysis should not be considered only as a method but mainly as a new paradigm involving the universal scaling pattern across all scales from smallest to largest [1].
In this direction, Jiang and Yin (2014) [18] introduced the 'ht-index' as a metric of complexity in the dataset's hierarchical structure, which fractal dimension alone could not express effectively. The authors experimented successfully with a grey-scale image, a digital elevation model, and road networks of very different cities. More recently, Karydas (2020) [14] employed head/tail breaks to develop a new method for selecting optimal scales for multi-scale segmentation of satellite imagery as a tool for facilitating a bottom-up approach in object-based image classification. The method employs fractal dimension as a metric of scale invariance in segmentation outputs; a topological transformation allows linking the partition with the segmentation process.
The main aim of this research was to develop a new, rapid, easy, objective, and robust methodology for selecting optimal scales when mapping topographic or hydrographic features. According to Meentemeyer (1989) [19], optimal scales are those where the most influencing variables of the mapped features remain invariant with rescaling; Sun et al. (2006) [16] argue that invariance over rescaling can be indicated using fractal dimension. Thus, the basic hypothesis of this research was that selection of optimal scales when mapping topographic or hydrographic features could be guided using fractal dimension as a judging metric, which, if it remains constant over rescaling, is expected to indicate optimal scales. The new method expands the method of Karydas (2020) [14], from object-based image analysis to multi-scale segmentation of digital elevation models (DEM) and their derivatives.

Methodology
The new method detects scale optimality in the structure of topographic or hydrographic datasets and then defines optimally scaled patterns in the spatial context by combining three distinct techniques:

•
Partition of topographic or hydrographic data with the rank-size rule (also known as head/tail breaks), thus determining the hierarchical data structure [20]. • Segmentation of the topographic or hydrographic surfaces using the Fractal Net Evolution Assessment (FNEA), thus converting thematic cells into meaningful objects [21]. • Linkage of the above two processes through a topological transformation, which projects the partitioning results to the segmented topographic or hydrographic data [14].
Partition with the rank-size rule is a process through which a dataset is divided recursively into two subsets: one containing all the values above the mean, called the 'head group', and one containing all the values below the mean, called the 'tail group' [22]. Partitioning allows determination of the inherent hierarchical structure by the data itself. The possibility of partitioning a dataset with the rank-size rule is associated with heavy-tailed distributions, i.e., datasets where far more small values can be found than large ones [6,23]. Types of heavy-tailed distributions include power-law, exponential, or lognormal distributions.
Segmentation is an appropriate rescaling methodology of a raster dataset, as the local heterogeneity of the property of interest is minimized, thus indicating those scales at which meaningful objects (i.e., naturally delineated surfaces) can be detected. One of the most common segmentation algorithms is the Fractal Net Evolution Assessment (FNEA), also known as 'multiresolution segmentation'. With FNEA, objects are created by minimizing the internal heterogeneity of contained cell values through a stepwise region-growing procedure until a user-defined threshold (called the 'scale factor') is reached [21]. The scale factor influences the desired object size and shape, according to a spectral and/or two shape criteria [24]. The spectral criterion applied with image datasets is defined as the standard deviation of the spectral values contained in a candidate object [25]. In this research, however, the so-called 'spectral criterion' is used to measure standard deviation of the relevant thematic values (e.g., elevation, slope, aspect, etc.) instead of spectral values.
FNEA has proven to be one of the most successful image segmentation algorithms in a GEOBIA framework [10]. Bobick and Bolles (1992) [26] argue that the information stored in a layer processed with FNEA is considered as fractal, in terms of holding the same degree of non-regularity at all scales or (inversely) self-similarity across every scale. Previous research findings suggest that when FNEA is applied with satellite images, a power-law equation of scale factor vs. mean object size is derived [14,27,28]: where f : scale factor; s f : the mean object size (at scale f (in m 2 ); a, b: constants (a in m −2 and b dimensionless). Newman (2005) [29] argues that "power-law behavior may be identified empirically using already seen results". Furthermore, Vuduc (1997) [30] states that fractal objects follow, by default, power-law statistical distributions of scaling. Following the paradigm of Karydas (2020) [14] for satellite images, a power-law equation between scale factor and mean object size for each of the topographic or hydrographic features of interest can be extracted for a study area, using a series of data segmentation-for example, the first twelve powers of 2, i.e., {1, 2, 4, 8, . . . , 2048}; the experimentation with imagery shows that such a series is enough to indicate a quite accurate function (with a high coefficient of determination, R 2 > 0.9).
A topological transformation was introduced to transfer optimal scales detected in the frequency distribution with the head/tail breaks to the physical space. This transformation is justified by the fact that the surface corresponding to a head group (produced with the rank-size rule), which in reality is dispersed into smaller pieces throughout the entire studied surface, could form a single object if these small pieces were contiguous (a kind of virtual aggregation); this single object would have the optimal size, provided that it is the result of an optimal scale segmentation. Hofmann et al. (1998) [31] argue that optimal information for further processing and analysis is expected to be provided by the best segmentation results. According to Karydas (2020) [14], the size which an object is supposed to have after the above transformation can be called the 'simulated mean object size' (SMOS), and the number of objects corresponding to SMOS can be called the 'simulated object number' (SON). Practically, SMOS is computed for every partition level by dividing the total extent of the study area over the extent of the head group (corresponding to that partition level), while SON is equal to the ratio of the head group extent over SMOS: where, s n : simulated mean object size (SMOS) at n-th partition (dimensionless); i: extent of the entire study area (m 2 ); h n : extent of the head group at n-th partition (m 2 ); N n : simulated object number (SON). This transformation allows for the projection of the spatial outputs of the partitioning process (after the recursive means) to the extent of the entire study area (where segmentation is applied). In other words, the physical space has been projected to a simulated one and this kind of projection can be considered as a mapping tool, transferring scale optimality from a (global) hierarchical structure to (local) spatial patterns ( Figure 1).
(SON). This transformation allows for the projection of the spatial outputs of the partitioning process (after the recursive means) to the extent of the entire study area (where segmentation is applied). In other words, the physical space has been projected to a simulated one and this kind of projection can be considered as a mapping tool, transferring scale optimality from a (global) hierarchical structure to (local) spatial patterns ( Figure 1). Finally, the computed SMOSs resulting from the topological transformation are inserted in Equation (2) to derive the optimal scale factors for the segmentation. To verify that the computed scales are the optimal ones, the fractal dimensions of all the successive scale factor values are computed and must be found to remain constant.
The whole computational procedure of the new method can be accomplished in three distinct steps (see also Figure 2): 1. Systematic segmentation of the surface of interest and extraction of a trend line (corresponding to Equation (2)) using the Fractal Net Evolution Assessment (FNEA) algorithm (usually embedded in an OBIA program). 2. Implementation of partition with the rank-size rule and computation of the head/tail breaks in a common spreadsheet; simple geospatial functions embedded in a geographic information system (GIS) program can be used to compute the means of the consecutive partitions. 3. Computation of the optimal scale factors using the trend equation on the head/tail breaks that resulted from the partition; the simulated mean object size (SMOS) is computed as an intermediate output; a common spreadsheet is enough to carry out these simple computations. Finally, the computed SMOSs resulting from the topological transformation are inserted in Equation (2) to derive the optimal scale factors for the segmentation. To verify that the computed scales are the optimal ones, the fractal dimensions of all the successive scale factor values are computed and must be found to remain constant.
The whole computational procedure of the new method can be accomplished in three distinct steps (see also Figure 2):

1.
Systematic segmentation of the surface of interest and extraction of a trend line (corresponding to Equation (2)) using the Fractal Net Evolution Assessment (FNEA) algorithm (usually embedded in an OBIA program).

2.
Implementation of partition with the rank-size rule and computation of the head/tail breaks in a common spreadsheet; simple geospatial functions embedded in a geographic information system (GIS) program can be used to compute the means of the consecutive partitions.

3.
Computation of the optimal scale factors using the trend equation on the head/tail breaks that resulted from the partition; the simulated mean object size (SMOS) is computed as an intermediate output; a common spreadsheet is enough to carry out these simple computations.
As it has become clear, partition process by the rank-size rule employs the mean, whereas segmentation process employs standard deviation. The values of the entire dataset under partitioning follow a heavy-tailed distribution, whereas the values contained by every object resulting from segmentation follows a normal distribution. With regard to the spatial context of the two processes, head/tail breaks is by default a non-spatial process applied on a global scale, whereas segmentation is a geospatial process applied on a local scale. Here is where the introduced topological transformation links the two procedures in order to transfer the partition results from the frequency plot to the segmentation outputs, thus converting a hierarchical structure into a spatial one (Table 1).  As it has become clear, partition process by the rank-size rule employs the mean, whereas segmentation process employs standard deviation. The values of the entire dataset under partitioning follow a heavy-tailed distribution, whereas the values contained by every object resulting from segmentation follows a normal distribution. With regard to the spatial context of the two processes, head/tail breaks is by default a non-spatial process applied on a global scale, whereas segmentation is a geospatial process applied on a local scale. Here is where the introduced topological transformation links the two procedures in order to transfer the partition results from the frequency plot to the segmentation outputs, thus converting a hierarchical structure into a spatial one (Table 1).

Study Area and Data Preparation
The two largest river basins in Greece, namely Pinios river basin and Acheloos river basin, located in the central-east and the central-west part of Greece, respectively, were selected to experiment with the new methodology. The Pinios river basin has a projected surface of about 10,701 km 2 with a perimeter of about 888 km and pours out to the Aegean Sea (eastern coast of mainland Greece); the area is mostly plain covered by agricultural crops. The Acheloos river basin has a projected surface of about 5686 km 2 with a perimeter of about 823 km and pours out to the Ionian Sea (western coast of mainland Greece); the area is covered mainly by forested mountains (Figure 3).

Study Area and Data Preparation
The two largest river basins in Greece, namely Pinios river basin and Acheloos river basin, located in the central-east and the central-west part of Greece, respectively, were selected to experiment with the new methodology. The Pinios river basin has a projected surface of about 10,701 km 2 with a perimeter of about 888 km and pours out to the Aegean Sea (eastern coast of mainland Greece); the area is mostly plain covered by agricultural crops. The Acheloos river basin has a projected surface of about 5686 km 2 with a perimeter of about 823 km and pours out to the Ionian Sea (western coast of mainland Greece); the area is covered mainly by forested mountains (Figure 3). An ASTER GDEM V2 (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model, Version 2) product was used as the only input dataset (freely downloaded from: https://earthdata.nasa.gov/). Digital elevation models (DEMs) are used widely for river basins' mapping and analysis [8,32]. An ASTER GDEM V2 (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model, Version 2) product was used as the only input dataset (freely downloaded from: https://earthdata.nasa.gov/). Digital elevation models (DEMs) are used widely for river basins' mapping and analysis [8,32].
The ASTER GDEM product was generated with stereo-pair images acquired with the ASTER instrument onboard of Terra satellite (EOS AM-1) and was released in October 2011 by the Ministry of Economy, Trade, and Industry (METI) of Japan and the National Aeronautics and Space Administration (NASA) of the USA. The ASTER GDEM products are provided in GeoTIFF format in a 30-m pixel size and 1 × 1 degree tiles.
It must be noted, however, that the original values in ASTER GDEM are estimated through optical instruments, i.e., those recording reflectance measurements and thus are less precise compared to LIDAR products (light detection and ranging, or laser imaging, detection, and ranging). LIDAR uses laser remotely sensed data in regular point grids and can be used to create 3D surface representations too [33,34].
The experimentation focused on four common topographic or hydrographic features, namely elevation, slope, aspect, and flow accumulation. The required tiles of the original ASTER GDEM products were mosaicked and clipped to the outlines of the two studied river basins. Elevation values were taken identical to the respective original cell values of the DEMs (in meters of altitude), while slope, aspect, and flow accumulation were computed from the DEM values in 3x3 windows using Spatial Analyst in ArcGIS program; the D8 algorithm was applied for flow accumulation estimation [32,35]). Ariza-Villaverde et al. (2015) [36] argue that the D8 algorithm is a reasonable method when convergent flow conditions predominate and has shown to be consistent among different flow patterns.
Stream and watershed networks in the two river basins were mapped by thresholding the flow accumulation (raster) layers (a process called 'flow accumulation thresholding' or FAT) [37]. The output layers can be as many as the number of the different thresholds selected and usually are converted and stored in vector format (polylines and polygons for streams and watershed, respectively). Here, FAT = 10,000 was applied to extract the river network in both river basins for geovisualization purposes only; selection of any FAT threshold for the display would not affect the analysis (neither the partition nor the segmentation processes), which are applied on the whole of the flow accumulation cells. The Pinios river basin was found to have an average elevation of 421.5 m, a maximum value of 2810 m, and a coefficient of variation of 86.5%; the mean slope of the basin is 16.7%, with a coefficient of variation of 102%; the mean aspect of the basin is 167.5, with a coefficient of variation of 63.2%. The slightly eastward aspect, compared to the expected mean of 180•, means that the mountain ranges lying on the west side are predominant in the basins' delineation. The shape index of the basin was computed to 2.147. Shape index is a metric of geometric complexity of an object and is defined as the ratio of the perimeter over four times the square root of the object's area; higher values correspond to higher complexity; rectangular is the least complex (or most regular) shape, with a shape index value of 1 [25] (Figure 4).
The Acheloos river basin was found to have an average elevation of 811.8 m, a maximum value of 2398 m, and a coefficient of variation of 63.2%%; the mean slope of the basin is 34.8%, with a coefficient of variation of 69.5%; the mean aspect of the basin is 168.9, with a coefficient of variation of 62.5%. The slightly eastward aspect, compared to the expected mean of 180•, also means that the mountains on the west side are predominant in the basins' delineation. The shape index of the basin was computed to 2.729 ( Figure 5). Some descriptive information and metrics of the two studied river basins are summarized in Table 2. analysis (neither the partition nor the segmentation processes), which are applied on the whole of the flow accumulation cells. The Pinios river basin was found to have an average elevation of 421.5 m, a maximum value of 2810 m, and a coefficient of variation of 86.5%; the mean slope of the basin is 16.7%, with a coefficient of variation of 102%; the mean aspect of the basin is 167.5, with a coefficient of variation of 63.2%. The slightly eastward aspect, compared to the expected mean of 180•, means that the mountain ranges lying on the west side are predominant in the basins' delineation. The shape index of the basin was computed to 2.147. Shape index is a metric of geometric complexity of an object and is defined as the ratio of the perimeter over four times the square root of the object's area; higher values correspond to higher complexity; rectangular is the least complex (or most regular) shape, with a shape index value of 1 [25] (Figure 4).  The Acheloos river basin was found to have an average elevation of 811.8 m, a maximum value of 2398 m, and a coefficient of variation of 63.2%%; the mean slope of the basin is 34.8%, with a coefficient of variation of 69.5%; the mean aspect of the basin is 168.9, with a coefficient of variation of 62.5%. The slightly eastward aspect, compared to the expected mean of 180•, also means that the mountains on the west side are predominant in the basins' delineation. The shape index of the basin was computed to 2.729 ( Figure 5). Some descriptive information and metrics of the two studied river basins are summarized in Table 2.

Parameter Pinios River Basin
Acheloos river Basin

Preconditions and Method Implementation
The new methodology was applied only with elevation, slope, and flow accumulation layers, excluding the aspect layers. The reason is that the aspect values were distributed almost constantly (or only with slight differentiations) throughout the entire range of possible values (i.e., 0-360 • ), thus demonstrating an almost straight and horizontal frequency distribution. Therefore, aspect did not meet the fundamental condition for applying the rank-size rule, i.e., that the data should follow a heavy-tailed distribution. This limitation is expected also to exist for a single mountain, and only in case of a small hillslope could the partition be implemented ( Figure 6).
For each of the elevation, slope, and flow accumulation layers, a series of twelve segmentations with the FNEA algorithm was carried out using eCognition Developer program. Plotting scale factor vs. mean object size allowed for the extraction of a predictive equation, in which the computation of the optimal scale factors (dependent variable, y) will be based by inputting SMOS (as the independent variable, x). The plotted segmentation results and the predictive equations (numerically and as trend lines) for the six layers are shown in Figure 7.
excluding the aspect layers. The reason is that the aspect values were distributed almost constantly (or only with slight differentiations) throughout the entire range of possible values (i.e., 0-360 o ), thus demonstrating an almost straight and horizontal frequency distribution. Therefore, aspect did not meet the fundamental condition for applying the rank-size rule, i.e., that the data should follow a heavy-tailed distribution. This limitation is expected also to exist for a single mountain, and only in case of a small hillslope could the partition be implemented (Figure 6). For each of the elevation, slope, and flow accumulation layers, a series of twelve segmentations with the FNEA algorithm was carried out using eCognition Developer program. Plotting scale factor vs. mean object size allowed for the extraction of a predictive equation, in which the computation of the optimal scale factors (dependent variable, y) will be based by inputting SMOS (as the independent variable, x). The plotted segmentation results and the predictive equations (numerically and as trend lines) for the six layers are shown in Figure 7. Then, the successive partition of all the layers according to the rank-size rule was implemented. The extents resulting from the head/tail breaks of the partition were recorded and the SMOS and SON were computed according to Equations (3) and (4), respectively. Three obvious prerequisites were taken into account for the determination of the optimal scale factors: a) SMOS had to be larger than the cell size (here, 900 m 2 ), b) SON had to be equal or larger than 1, and c) the computed scale factors had to be larger than zero.

Results
The outputs of the new method implementation for each of the studied properties (elevation, slope, and flow accumulation) for the studied river basins (Pinios and Acheloos) are presented in Tables 3, 4, 5, 6, 7, and 8.  Then, the successive partition of all the layers according to the rank-size rule was implemented. The extents resulting from the head/tail breaks of the partition were recorded and the SMOS and SON were computed according to Equations (3) and (4), respectively. Three obvious prerequisites were taken into account for the determination of the optimal scale factors: (a) SMOS had to be larger than the cell size (here, 900 m 2 ), (b) SON had to be equal or larger than 1, and (c) the computed scale factors had to be larger than zero.

Results
The outputs of the new method implementation for each of the studied properties (elevation, slope, and flow accumulation) for the studied river basins (Pinios and Acheloos) are presented in Tables 3-8. * optimal scale factors in bold. The fractal dimension of the consecutive partition pairs for all the properties of interest were derived using the Mandelbrot equation (Mandelbrot 1982) [17], adapted to this method by Karydas 2020 [14], as follows:     * optimal scale factors in bold.

Discussion
The computed optimal scale factors for elevation or slope segmentations as well as the corresponding mean object size were found to be quite small for landform units. Indicatively, even for the highest computed scale factor for elevation in the Pinios river basin (i.e., 16.9), the corresponding mean object size was found to be only about 16 ha (Figure 8). Similarly, for the Pinios slope, the highest computed scale factor (i.e., 8.4) resulted in a 7-ha mean object size. In contrast to that, the computed optimal scale factors for flow accumulation are shown to be much larger than those computed for elevation and slope for both river basins. This can be attributed to the fact that flow accumulation layers are 'empty' of high values or inversely full of 0 values to a large extent. In other words, high values follow linear patterns, leaving big homogeneous objects of low values to fill the surface everywhere else. As a result of this particular spatial pattern, the mean object size for flow accumulation was found to be much larger than that of the other two feature types, and consequently, the scale factor was found to be similarly large. The total number of rows in the computational tables (Tables 3-8) corresponds to another fractal metric known as the ht-index, introduced by Jiang and Yin (2014) [18]. Ht-index is an integer indicating the number of hierarchical levels of a dataset and, as such, a metric of its hierarchical structure complexity. In their work, the authors estimated the ht-index of the elevation in the United States of America to 16 using a DEM of 1800-m resolution. Similar ht-index values (14 or 15) were computed in the current study for elevation and slope too. The difference between the two resolutions (30 vs. 1800 meters) should not be considered as an obstacle for comparison, as the geographic scale between the two case studies also differs greatly. Significantly lower values were found for flow accumulation in the current study (four in both cases). However, no other similar studies have been made available to do a comparison. The overall results from the two river basins indicate high complexity for the elevation and slope surfaces and low complexity for the flow accumulation surfaces (Table 9).  The total number of rows in the computational tables (Tables 3-8) corresponds to another fractal metric known as the ht-index, introduced by Jiang and Yin (2014) [18]. Ht-index is an integer indicating the number of hierarchical levels of a dataset and, as such, a metric of its hierarchical structure complexity. In their work, the authors estimated the ht-index of the elevation in the United States of America to 16 using a DEM of 1800-m resolution. Similar ht-index values (14 or 15) were computed in the current study for elevation and slope too. The difference between the two resolutions (30 vs. 1800 m) should not be considered as an obstacle for comparison, as the geographic scale between the two case studies also differs greatly. Significantly lower values were found for flow accumulation in the current study (four in both cases). However, no other similar studies have been made available to do a comparison. The overall results from the two river basins indicate high complexity for the elevation and slope surfaces and low complexity for the flow accumulation surfaces (Table 9).
In some cases (for example, in Pinios slope mapping), noticeable divergences were detected between the fractal dimension derived from the Mandelbrot (1982) [17] equation (refer to Equation (3)) and that inferred from the 'slope' of the trend line, as defined by the plots of segmentations series (refer to Figure 4). This deviation was expected and could be attributed to the lower (compared to that of the other features) coefficient of determination (R 2 ) of the equations derived from the plots. In all cases, the fractal dimension values computed with Equation (3) are the most precise ones as they are derived from the optimal data, whereas the trend lines correspond only to some systematic but not necessarily optimal segmentations. Similar deviations were noticed by Karydas (2020) [14] in the original implementation of the method with a variety of satellite images. However, the implementation of the new method with the DEM and its derivatives revealed more problems than that with the imagery. As it can be understood from the distribution plots of Acheloos elevation and slope, their frequency distribution does not follow an ideal heavy-tailed shape. More specifically, the Acheloos elevation distribution demonstrates two additional peaks and valleys, while the Acheloos slope distribution imitates a normal distribution skewed to the left. In contrast, Pinios elevation demonstrates only a small extra peak, while both the slope distributions are closer to heavy-tailed shapes than any other else ( Figure 9). Thus, only the flow accumulation distributions of both river basins were of an ideal heavy-tailed shape.
Another problem was that the rank-size rule in both river basins resulted in many ups and downs towards reaching the end-up threshold, while some of the head groups were detected close to 50% early in the partitioning process ( Figure 10). In such cases, it is debatable if partitioning could be applied or better be avoided. Note that for digital elevation models, the end-up threshold is set usually to less than or equal to 40% (Jiang 2015) [19].
Nonetheless, regardless of the above abnormalities, the current results indicate clearly that the new method was successful in all cases. Moreover, plotting segmentation vs. partitioning data of any of the studied features shows that optimal scale factors could be computed directly from the mean values derived from the partitions. Indicatively, the two studied river basins resulted in the following equations for elevation (see also Figure 11): For Pinios : f = 4·10 −13 e 0.0111M n R 2 = 0.9023 (6) For Acheloos : f = 2·10 −18 e 0.0184M n R 2 = 0.8046 (7) where f : scale factor; M n : mean at partition n. If Equations (4) and (5) were generalized into a single one for any river basin, optimal scale factors for each feature could be computed directly from the successive mean values of the partition process alone and, thus, the segmentation process could be skipped. Considering that scale factor is equivalent to standard deviation within candidate objects, Equations (4) and (5), in reality, link the two main statistics in every dataset. Furthermore, high correlation between other series of data used in this study has also been detected (e.g., head extents with computed scale factors). This could be attributed to the fact that head groups correspond to the successive means of the dataset, while segmented objects at different scales correspond to the successive standard deviations of the created objects.
shape. More specifically, the Acheloos elevation distribution demonstrates two additional peaks and valleys, while the Acheloos slope distribution imitates a normal distribution skewed to the left. In contrast, Pinios elevation demonstrates only a small extra peak, while both the slope distributions are closer to heavy-tailed shapes than any other else ( Figure 9). Thus, only the flow accumulation distributions of both river basins were of an ideal heavy-tailed shape. Another problem was that the rank-size rule in both river basins resulted in many ups and downs towards reaching the end-up threshold, while some of the head groups were detected close to 50% early in the partitioning process ( Figure 10). In such cases, it is debatable if partitioning could be applied or better be avoided. Note that for digital elevation models, the end-up threshold is set usually to less than or equal to 40% (Jiang 2015) [19].   Flow accumulation means Equations (4) and (5), in reality, link the two main statistics in every dataset. Furthermore, high correlation between other series of data used in this study has also been detected (e.g., head extents with computed scale factors). This could be attributed to the fact that head groups correspond to the successive means of the dataset, while segmented objects at different scales correspond to the successive standard deviations of the created objects.
(a) (b) Figure 11. Indicative log-log plots of optimal scale factor vs. mean elevation for the two studied river basins; a) Pinios and b) Acheloos.

Conclusion
In this research, a new method for scale optimization when mapping elevation, slope, or flow accumulation was developed and experimented with the two biggest river basins in Greece. Aspect does not meet fundamental conditions for applying the rank-size rule and thus the new method could not be applied with aspect layers.
The new method uses head/tail breaks to indicate the optimal scales in the hierarchical data structure and a topological transformation to project the results to the spatial context. Fractal dimension is employed to confirm optimality in the spatial data. The method expands the work of Karydas (2020) [14], for multi-scale segmentation of image data, to topographic and hydrographic raster datasets.
The experimentation revealed some deviations from the ideal conditions regarding heavy-tailed distributions and partitioning. Regardless of the possibility of not meeting the ideal conditions, however, the method proved successful to identify a set of optimal scales in all cases. Therefore, the basic hypothesis of the research, that selection of optimal scales in topographic and hydrographic mapping could be guided by fractal dimension as a scale optimality metric, was verified in all terms. Scale factor (f) Mean elevation at partition n (m) Figure 11. Indicative log-log plots of optimal scale factor vs. mean elevation for the two studied river basins; (a) Pinios and (b) Acheloos.

Conclusions
In this research, a new method for scale optimization when mapping elevation, slope, or flow accumulation was developed and experimented with the two biggest river basins in Greece. Aspect does not meet fundamental conditions for applying the rank-size rule and thus the new method could not be applied with aspect layers.
The new method uses head/tail breaks to indicate the optimal scales in the hierarchical data structure and a topological transformation to project the results to the spatial context. Fractal dimension is employed to confirm optimality in the spatial data. The method expands the work of Karydas (2020) [14], for multi-scale segmentation of image data, to topographic and hydrographic raster datasets.
The experimentation revealed some deviations from the ideal conditions regarding heavy-tailed distributions and partitioning. Regardless of the possibility of not meeting the ideal conditions, however, the method proved successful to identify a set of optimal scales in all cases. Therefore, the basic hypothesis of the research, that selection of optimal scales in topographic and hydrographic mapping could be guided by fractal dimension as a scale optimality metric, was verified in all terms.
Optimal scales known in advance of environmental problem solving and land management tasks would reveal underlying (or even hidden) processes, thus supporting decision making. Furthermore, the new method may define appropriate mapping scales for object generalization within the framework of topographic and terrain analysis, where an appropriate proportion of variability among feature densities and textures has to be retained. For example, the object-scale selection process could be equated to a filtering process for near real-time web-mapping.
Implementation of the new method prerequisites an object-based analysis program (to conduct segmentation tasks), a GIS-program (to prepare the required layers and compute spatial statistics), and a spreadsheet (to carry out the computations of the partial and final outputs).
Regarding the studied river basins, it was shown that the Pinios river basin was more heterogeneous for elevation and less for slope from that of Acheloos, while they were similarly heterogeneous for the aspect. Both river basins were found to be similarly complicated, as they had similar values of ht-index. The fractal dimension differed for slope between the two basins, whereas it was similar for elevation and flow accumulation.
In summary, the new method proved to be rapid (short time for accomplishment), easy (simple processes), objective (fully numerical), and robust (results were derived even with non-ideal preconditions). The method could be further simplified by directly using the mean values of the partition to compute the optimal scale factors, provided that more tests with diverse river basins would take place towards generalization. In such a case, partition would be the only required process to carry out, so that data plotting would result in the necessary non-linear equations. In its general form, the method could help identify the best datasets for extracting features of a particular size.
Finally, new light must be shed on the properties of heavy-tailed distributions with a view to clarify and further specify the required preconditions. Future research also should expand the method to the vector data type of hydrographic features, for example, to stream networks and watersheds. Acknowledgments: Segmentation tasks were conducted with the eCognition Developer program in demonstration mode; geospatial analysis was conducted with ArcGIS Pro; statistics and graphs were produced with Excel spreadsheets. Many thanks to the anonymous reviewers, whose substantial critique helped us improve the paper a lot.