Heuristic and Numerical Geometrical Methods for Estimating the Elevation and Slope at Points Using Level Curves. Application for Embankments

Abstract

Both the calculation of ground slopes at points on the map and the elevation estimation for these points bear significant importance and also have applications in various domains, such as civil engineering, road and railway design. The paper presents two methods that use level curves: one that is fast and approximate and another which is slower, but more precise. The running speed of the two proposed methods and their results are compared by performing 100 million experiments. The paper also presents how these methods can be applied to optimize embankments. An accurate method to calculate the horizontal plane of the excavation/filling when building a new house is also presented.

1. Introduction

The possibility of calculating the slope and elevation using level curves has practical applicability in several fields of activity, such as civil engineering [1], architecture [2,3], railway [4], road and bridge engineering [5,6,7,8], topography [9], etc. Although the results found can be used in various areas, this article is limited to solving two problems related to civil engineering and to railways, roads, and bridges, both ultimately related to the embankments optimization. The problem is well known, but it was not closely examined by design engineers, who preferred an intuitive, approximate solution. The current article presents a rigorous and mathematically accurate solution to find the optimum between excavations and fillings, resulting in a reduction of the embankments related costs. The article presents in detail the solutions of two concrete problems in two different fields of activity.

The first problem investigated is in the field of civil engineering and consists of determining the optimal elevation corresponding to the finished floor of the ground level, when the building has to be positioned on a slope [1,2,3]. The second problem investigated is in the design of a new railway or road when there are several fixed points on the topographic map through which the new road [10,11,12,13] or the new railway [4,11] will pass.

The designer has to know the slope and elevation in each of these points as accurately as possible. Using the methods proposed in this paper, the two parameters (slope and elevation) can be determined fast and with a high level of accuracy using the level curves of the map, without the need for direct measurement on the spot. Furthermore, when designing roads or railways, the risk of failure, the frictional resistance, the safety factor against sliding or failure, the translational fill failure, or the cut volume and the road fill per meter require knowing the slopes (side slope, cut slope, or fill slope).

The paper is organized as follows: an algorithm to calculate the ground slope at a given point, using level curves, is deduced in Section 2. This problem is reduced to estimating the shortest segment that connects two lines through a point. Two methods are proposed, a fast and approximate one and a more exact one, respectively. The two methods are compared, and it is shown that the errors of the first method are low. In Section 3, the optimum excavation/filling plane is calculated when a new house is built. Two methods are proposed and discussed in Section 4, where they are compared. To finalize the paper, conclusions will be drawn based on the arguments presented.

2. Materials and Methods

A method of calculating the soil slope through a given point is being introduced. The input data are: the level curves, in Figure 1, the two level curves located at elevation 610 and, respectively, 620, and the 2D coordinates of a point Q on the topographic map. The direction of the Ox axis is considered from west to east, the direction of the Oy axis is south-north. The output is the slope calculated at point Q. The slope can be determined using two other points collinear with Q for which the elevation is known. The points on the map for which the elevation is known are located on level curves. Therefore, it is enough to find two points collinear with Q located on two different level curves. Of course, for accuracy reasons, it is preferable that Q is located between these 2 points and the distance between the chosen 2 points is minimum. Consequently, the problem can be accurately solved if we determine the closest two level curves to Q (Q is located in between these curves) and two points P₁ and P₂, each located on one of the level curves, and the distance between P₁ and P₂ is minimum (see Figure 1).

Since a level curve is approximated with consecutive coplanar segments of straight lines, our problem can be reduced to finding the shortest length segment |P₁P₂| that connects two coplanar line segments located on two consecutive level curves. When the two 2D points P₁(x₁,y₁) and P₂(x₂,y₂) are determined, they are transformed into 3D points by adding the elevations z₁ and, respectively, z₂ of the level curves on which each of the points are located. Using the 3D coordinates of these 2 points, the slope m = tan(α) at point Q (see Figure 2) can be calculated as follows:

$$\{\begin{array}{c}\mathrm{m}=\frac{\left|{\mathrm{z}}_2-{\mathrm{z}}_1\right|}{\sqrt{{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2}}\\ \mathsf{\alpha }=\mathrm{atan}2\left(\left|{\mathrm{z}}_2-{\mathrm{z}}_1\right|,\sqrt{{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2}\right) \end{array}$$

(1)

The first problem (denoted IPS) intended to be solved is to identify the pair of segments (located on two level curves) to which the finding of the shortest length segment should be applied. Of course, for each of the two level curves, the candidates are the closest segments to point Q. Let us denote by S_i the set of candidates (segments) for curve i (i = 1, 2). For each pair of segments s₁ ∈ S₁ and s₂ ∈ S₂, respectively, the corresponding shortest length segment [P₁P₂] between s₁ and s₂ is calculated. If s₁ contains P₁ and s₂ contains P₂, then the pair (s₁, s₂) is added to the set (denoted C) of solution candidates for IPS (Figure 3).

The set S_i (i = 1, 2) can be built as follows: first, we find the point T belonging to the level curve i that is closest to Q. We consider then k ≥ 1 points (L₁, L₂, …, L_k) to the left of T on the level curve and k points (R₁, R₂, …, R_k) located to right (see Figure 4). Finally, the set S_i consists of 2k segments of straight lines given by the consecutive points on the level curve i: L_k, L_k-1, …, L₁, T, R₁, R₂, …, R_k, i.e.,:

$${\mathrm{S}}_{\mathrm{i}}=\{\left[{\mathrm{L}}_{\mathrm{h}+1}{\mathrm{L}}_{\mathrm{h}}\right]|\mathrm{h}=1,2,\dots ,\mathrm{k}-1\}{{\displaystyle \cup }}^{}\left\{\left[{\mathrm{L}}_1\mathrm{T}\right],\left[{\mathrm{TR}}_1\right]\right\}{{\displaystyle \cup }}^{}\{\left[{\mathrm{R}}_{\mathrm{h}}{\mathrm{R}}_{\mathrm{h}+1}\right]|\mathrm{h}=1,2,\dots ,\mathrm{k}-1\}$$

(2)

Now, let us see how k (giving the number of points in the vicinity of T) is considered. Usually (in real tests), k is enough to be 1. In few cases k is found to be greater. So, the calculations are started with k = 1 and, if necessary, it is increased.

The Algorithm 1 for solving IPS (denoted AIPS) is:

Algorithm 1 Solving IPS (denoted AIPS)

AIPS: Input: the points of the closest two level curves to Q Output: the shortest line segment [P1P2] that connects the given two level curves finished = false; While not finished do   C = ϕ (empty set);   Build the sets of line segments S1 and S2 (see Equation (2));   For each segment s1 from S1 do      For each segment s2 from S2 do         Calculate the shortest length segment P1P2 that connects s1 and s2 and passes through Q;         If P1 ∈ s1 and P2 ∈ s2 then           Add [P1P2] to C;         End if;      End for;   End for;   If C not is empty then      finished = true;   else      k = k+1;   End if; End while;

The shortest line segment that connects both level curves is the shortest segment [P₁P₂] from C.

After applying AIPS, the problem is reduced to calculating the shortest segment that passes through a given point Q and connects two coplanar lines d₁ and d₂ (Figure 4). We shall denote this problem as SSTPC2L (shortest segment through a point that connects 2 lines). In Figure 4 the points A₁₁ and A₁₂ determine the line d₁, and the points A₂₁ and A₂₂ determine the line d₂.

We present two methods of solving SSTPC2L: an approximate method and an exact mathematical method. These two methods (presented in Section 2.1 and Section 2.2) were partially described in [7].

2.1. Heuristic (Approximate) Method for SSTPC2L

From the topographic map, the point Q is known and two segments A₁₁A₁₂ and A₂₁A₂₂ are determined belonging to the level curves from Figure 1 (A₁₁A₁₂ on curve 610 and A₂₁A₂₂ on curve 620). These two segments are found according to Figure 3. The method consists of drawing the perpendiculars from the given point Q on the two given lines d₁ and, respectively, d₂ (Figure 5). We construct a line passing through Q which is parallel to the line determined by the feet of the two perpendiculars. The obtained line is denoted by d₃, and its intersection with the lines d₁ and d₂ is the searched segment [P₁P₂] (approximate solution of SSTPC2L).

Now, let’s present how the points P₁ and P₂ are calculated.

The equation of the line d_i (i = 1, 2) starting from the points A_i1 and A_i2 is:

$$\mathrm{y}={\mathrm{m}}_{\mathrm{i}}\mathrm{x}+{\mathrm{n}}_{\mathrm{i}},{\mathrm{where}\mathrm{m}}_{\mathrm{i}}=\frac{{\mathrm{y}}_{\mathrm{i}2}-{\mathrm{y}}_{\mathrm{i}1}}{{\mathrm{x}}_{\mathrm{i}2}-{\mathrm{x}}_{\mathrm{i}1}}{\mathrm{and}\mathrm{n}}_{\mathrm{i}}={\mathrm{y}}_{\mathrm{i}1}-{\mathrm{m}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}1}$$

(3)

The perpendicular from Q on d_i (i = 1, 2) is:

$$\mathrm{y}=-\frac{1}{{\mathrm{m}}_{\mathrm{i}}}\mathrm{x}+{{\mathrm{n}}^{\prime }}_{\mathrm{i}},{\mathrm{where}{\mathrm{n}}^{\prime }}_{\mathrm{i}}={\mathrm{y}}_{\mathrm{Q}}+\frac{1}{{\mathrm{m}}_{\mathrm{i}}}{\mathrm{x}}_{\mathrm{Q}}$$

(4)

The foot ${\mathrm{Q}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}^{\mathrm{Q}},{\mathrm{y}}_{\mathrm{i}}^{\mathrm{Q}}\right)$ (i = 1, 2) of the perpendicular on d_i is obtained by solving the system of Equations (3) and (4).

The slope of the line Q₁Q₂ is:

$${\mathrm{m}}_{12}=\frac{{\mathrm{y}}_2^{\mathrm{Q}}-{\mathrm{y}}_1^{\mathrm{Q}}}{{\mathrm{x}}_2^{\mathrm{Q}}-{\mathrm{x}}_1^{\mathrm{Q}}}$$

(5)

Since ${\mathrm{m}}_{12}$ is the slope of ${\mathrm{Q}}_1{\mathrm{Q}}_2$, then the parallel through the point $\mathrm{Q}\left({\mathrm{x}}_{\mathrm{Q}},{\mathrm{y}}_{\mathrm{Q}}\right)$ is:

$$\mathrm{y}={\mathrm{y}}_{\mathrm{Q}}+{\mathrm{m}}_{12}\left(\mathrm{x}-{\mathrm{x}}_{\mathrm{Q}}\right)$$

(6)

Finally, the points P_i (i = 1, 2) are obtained by solving the system of Equations (3) and (6).

The following Algorithm 2 (denoted A1SSTPC2L) is obtained to solve the problem SSTPC2L:

Algorithm 2 Solving the problem SSTPC2L

A1SSTPC2L: Input: Q and the points of the closest two level curves to Q Output: The slope at point Q Apply AIPS to find A11, A12, A21, A22; Calculate the point Pi by solving system of Equations (3) and (6) (i = 1, 2). Calculate the slope at point Q using Equation (1).

2.2. The Exact (Mathematical) Method for SSTPC2L

For this method, in order to simplify calculus, we apply some initial transformations (rotations and translations) to the points A₁₁, A₁₂, A₂₁, A₂₂, and Q so that the bisector of the angle between the lines d₁ and d₂ becomes parallel to Oy axis and the point Q is in the origin O, i.e., Q = O (0,0). To do so, we first calculate the intersection I of d₁ and d₂:

$${\mathrm{d}}_1{{\displaystyle \cap }}^{}{\mathrm{d}}_2=\left\{\mathrm{I}\right\}$$

(7)

by solving the system of linear equations:

$$\{\begin{array}{c}{\mathrm{d}}_1:\frac{\mathrm{x}-{\mathrm{x}}_{11}}{{\mathrm{x}}_{12}-{\mathrm{x}}_{11}}=\frac{\mathrm{y}-{\mathrm{y}}_{11}}{{\mathrm{y}}_{12}-{\mathrm{y}}_{11}}\\ {\mathrm{d}}_2:\frac{\mathrm{x}-{\mathrm{x}}_{21}}{{\mathrm{x}}_{22}-{\mathrm{x}}_{21}}=\frac{\mathrm{y}-{\mathrm{y}}_{21}}{{\mathrm{y}}_{22}-{\mathrm{y}}_{21}}\end{array}$$

The following transformations to A₁₁, A₁₂, A₂₁, A₂₂ and Q are applied:

• Translation with (−x_I, −y_I) (I is moved into origin O)
• Rotation with sin(−u) = −sin(u) and cos(−u) = cos(u), where:

  $$\{\begin{array}{c}\mathrm{r}=\sqrt{{\left({\mathrm{x}}_{12}-{\mathrm{x}}_{11}\right)}^2+{\left({\mathrm{y}}_{12}-{\mathrm{y}}_{11}\right)}^2}\\ \cos \left(\mathrm{u}\right)=\left({\mathrm{x}}_{12}-{\mathrm{x}}_{11}\right)/\mathrm{r}\\ \sin \left(\mathrm{u}\right)=\left({\mathrm{y}}_{12}-{\mathrm{y}}_{11}\right)/\mathrm{r}\end{array}$$

  (8)

After this rotation, the line d₁ is over Ox axis of the map’s system of coordinates.

3.

Rotation with sin(v) and cos(v), where:

$$\{\begin{array}{c}\mathrm{w}=\mathrm{atan}2\left({\mathrm{y}}_{22}-{\mathrm{y}}_{21}, {\mathrm{x}}_{22}-{\mathrm{x}}_{21}\right)\\ \cos \left(\mathrm{v}\right)=\cos \left(\mathsf{\pi }/2-\mathrm{w}/2\right)\\ \sin \left(\mathrm{v}\right)=\sin \left(\mathsf{\pi }/2-\mathrm{w}/2\right)\end{array}$$

(9)

After this rotation, the bisector of the angle given by lines d₁ and d₂ is over Oy.

4.

Translation with (−x_Q, −y_Q) (Q is moved into origin O)

The following 2 systems of equations have to be solved in order to find the points P_i(x_i,y_i) (i = 1,2) where the lines d_i (i = 1,2) intersect the line d:

$$\{\begin{array}{c}\mathrm{y}={\mathrm{a}}_{\mathrm{i}}\mathrm{x}+{\mathrm{b}}_{\mathrm{i}}\\ \mathrm{y}=\mathrm{ax}\end{array}$$

(10)

Since A_i1, A_i2 define the line d_i (i = 1,2) in (10) we have:

$${\mathrm{a}}_{\mathrm{i}}=\frac{\left({\mathrm{y}}_{\mathrm{i}2}-{\mathrm{y}}_{\mathrm{i}1}\right)}{\left({\mathrm{x}}_{\mathrm{i}2}-{\mathrm{x}}_{\mathrm{i}1}\right)}$$

(11)

$${\mathrm{b}}_{\mathrm{i}}={\mathrm{y}}_{\mathrm{i}1}-{\mathrm{a}}_{\mathrm{i}1}{\mathrm{x}}_{\mathrm{i}1}$$

(12)

Solving the 2 systems of equations, the points P_i(x_i,y_i) (i = 1,2) are found, where:

$${\mathrm{x}}_{\mathrm{i}}=\frac{{\mathrm{b}}_{\mathrm{i}}}{\mathrm{a}-{\mathrm{a}}_{\mathrm{i}}}$$

(13)

$${\mathrm{y}}_{\mathrm{i}}=\frac{{\mathrm{ab}}_{\mathrm{i}}}{\mathrm{a}-{\mathrm{a}}_{\mathrm{i}}}$$

(14)

The distance between P₁ and P₂ must be minimized, i.e.,:

$$\underset{\mathrm{a}}{\min }\mathrm{f}\left(\mathrm{a}\right)$$

(15)

where:

$$\mathrm{f}\left(\mathrm{a}\right)={\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2+{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2$$

(16)

By replacing (13) and (14) in (16) the following formula for f(a) is obtained:

$$\mathrm{f}\left(\mathrm{a}\right)={\left(\frac{{\mathrm{a}}_2{\mathrm{b}}_2}{\mathrm{a}-{\mathrm{a}}_2}-\frac{{\mathrm{a}}_1{\mathrm{b}}_1}{\mathrm{a}-{\mathrm{a}}_1}+{\mathrm{b}}_2-{\mathrm{b}}_1\right)}^2+{\left(\frac{{\mathrm{b}}_2}{\mathrm{a}-{\mathrm{a}}_2}-\frac{{\mathrm{b}}_1}{\mathrm{a}-{\mathrm{a}}_1}\right)}^2$$

(17)

We make the following notations:

$${\mathrm{a}}_{\mathrm{L}}=\min \left\{\frac{-1}{{\mathrm{a}}_1},\frac{-1}{{\mathrm{a}}_2}\right\}$$

(18)

$${\mathrm{a}}_{\mathrm{R}}=\max \left\{\frac{-1}{{\mathrm{a}}_1},\frac{-1}{{\mathrm{a}}_2}\right\}$$

(19)

It is easy to see that the function f is decreasing for $\mathrm{a}\le {\mathrm{a}}_{\mathrm{L}}$, and is increasing for $\mathrm{a}\ge {\mathrm{a}}_{\mathrm{R}}$, since f is continuous, ${\mathrm{a}}_{\mathrm{L}}$ and ${\mathrm{a}}_{\mathrm{R}}$ are the slopes of the perpendiculars QQ₁ and QQ₂ from Q on the lines d₁ and, respectively, d₂, and “a” is the slope of the line d₃. Consequently, there is a minimum of the function f denoted a_min in the interval $\left[{\mathrm{a}}_{\mathrm{L}},{\mathrm{a}}_{\mathrm{R}}\right]$. So, the value a_min is one of the solutions of the equation:

$${\mathrm{f}}^{\prime }\left(\mathrm{a}\right)=0$$

(20)

on the interval $\left[{\mathrm{a}}_{\mathrm{L}},{\mathrm{a}}_{\mathrm{R}}\right]$, where:

$${\mathrm{f}}^{\prime }\left(\mathrm{a}\right)=2\left(\frac{{\mathrm{a}}_2{\mathrm{b}}_2}{\mathrm{a}-{\mathrm{a}}_2}-\frac{{\mathrm{a}}_1{\mathrm{b}}_1}{\mathrm{a}-{\mathrm{a}}_1}+{\mathrm{b}}_2-{\mathrm{b}}_1\right)\left[\frac{-{\mathrm{a}}_2{\mathrm{b}}_2}{{\left(\mathrm{a}-{\mathrm{a}}_2\right)}^2}+\frac{{\mathrm{a}}_1{\mathrm{b}}_1}{{\left(\mathrm{a}-{\mathrm{a}}_1\right)}^2}\right]+2\left(\frac{{\mathrm{b}}_2}{\mathrm{a}-{\mathrm{a}}_2}-\frac{{\mathrm{b}}_1}{\mathrm{a}-{\mathrm{a}}_1}\right)\left[\frac{-{\mathrm{b}}_2}{{\left(\mathrm{a}-{\mathrm{a}}_2\right)}^2}+\frac{{\mathrm{b}}_1}{{\left(\mathrm{a}-{\mathrm{a}}_1\right)}^2}\right]$$

(21)

Remark 1.

Since the bisector of the angle between d₁ and d₂ is parallel to Oy axis, it is not difficult to see that:

$${\mathrm{a}}_{\mathrm{L}}=-{\mathrm{a}}_{\mathrm{R}}$$

(22)

In order to obtain the solution a_min of the Equation (20), a numerical method [14] can be used such as bisection method [15], or tangent method. Since ${\mathrm{f}}^{\prime }\left({\mathrm{a}}_{\mathrm{L}}\right)<0$ and ${\mathrm{f}}^{\prime }\left({\mathrm{a}}_{\mathrm{R}}\right)>0=>{\mathrm{f}}^{\prime }\left({\mathrm{a}}_{\mathrm{L}}\right)·{\mathrm{f}}^{\prime }\left({\mathrm{a}}_{\mathrm{R}}\right)<0$, the bisection method is suitable to be applied on the interval:

$$\mathrm{a}\in \left[{\mathrm{a}}_{\mathrm{L}},{\mathrm{a}}_{\mathrm{R}}\right]$$

(23)

After a_min is calculated, the points P_i (i = 1, 2) are obtained using (13) and (14) and the distance between these two points is:

$$\sqrt{\mathrm{f}\left({\mathrm{a}}_{\min }\right)}$$

(24)

In order to go back to the initial system of coordinates, the inverse initial transformations must be applied in reverse order to the points P₁ and P₂, i.e.,:

• Translation with (−x_Q, −y_Q)
• Rotation with sin(u) and cos(u) (were calculated in (8))
• Rotation with sin(-v) = sin(v) and cos(-v) = cos(v) (calculated in (9))
• Translation with (x_I, y_I).

The following Algorithm 3 (denoted A2SSTPC2L) is obtained to solve SSTPC2L:

Algorithm 3 Solving SSTPC2L

A2SSTPC2L: Input: Q and the points of the closest two level curves to Q Output: The slope at point Q Apply AIPS to find A11, A12, A21, A22; Transform the points A11, A12, A21, A22 and Q so that the bisector of the angle between the lines d1 and d2 is parallel to Oy axis and the point Q is in the origin O; Find the solution amin of the equation ${\mathrm{f}}^{\prime }\left(\mathrm{a}\right)=0$ using bisection method on the interval $\left({\mathrm{a}}_{\mathrm{L}},{\mathrm{a}}_{\mathrm{R}}\right)$ (see Equations (19) and (20)); Calculate the coordinates of the points Pi (i = 1,2) using (13) and (14), where a = amin; Apply the inverse initial transformations in reverse order to the points P1 and P2; Calculate the slope at point Q using (1).

We implemented the above algorithm in Visual C++ and in the bisection method we set the error err = 0.0001 for finding the value a_min. In Figure 6, a graphical output of our program is presented.

3. Results

Besides estimation of the ground slope for a given point, other two possible problems can be solved since we have now obtained a method to calculate the points P₁, and P₂ on the closest two level curves to a given point Q so that P₁, Q, and P₂ are collinear and the distance between P₁ and P₂ is minimum.

3.1. Elevation Estimation at a Point

The challenge is to estimate as well as possible the elevation z of a point Q. Assuming that z₁ < z₂ (see Figure 3) and since QQ’ is parallel to P₂P′₂, in the triangle P₁P₂P′₂ we have:

$$\frac{{\mathrm{z}}_{\mathrm{Q}}-{\mathrm{z}}_1}{{\mathrm{z}}_2-{\mathrm{z}}_1}=\frac{\mathrm{dist}\left({\mathrm{P}}_1,{\mathrm{Q}}^{\prime }\right)}{\mathrm{dist}\left({\mathrm{P}}_1,{\mathrm{P}}_2^{\prime }\right)}$$

(25)

From (25) we obtain the elevation of Q:

$${\mathrm{z}}_{\mathrm{Q}}={\mathrm{z}}_1+\sqrt{\frac{{\left({\mathrm{x}}_{\mathrm{Q}}-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_{\mathrm{Q}}-{\mathrm{y}}_1\right)}^2}{{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2}}\left({\mathrm{z}}_2-{\mathrm{z}}_1\right)$$

It is immediate that, without the assumption of z₁ < z₂, the following formula can be used to estimate the elevation of Q:

$${\mathrm{z}}_{\mathrm{Q}}=\min \left\{{\mathrm{z}}_1,{\mathrm{z}}_2\right\}+\sqrt{\frac{{\left({\mathrm{x}}_{\mathrm{Q}}-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_{\mathrm{Q}}-{\mathrm{y}}_1\right)}^2}{{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2}}\left|{\mathrm{z}}_2-{\mathrm{z}}_1\right|$$

(26)

3.2. The Slope between Two Points

The problem of estimating the slope between two points is significant, seeing that when designing roads, the slope between any two points on a road cannot exceed a given maximum slope, e.g., for highways this value has to be less than 6% or 7%. Therefore, the chosen points from the future road can be tested for this eligibility. To do that, we consider the elevations of the points Q₁ and Q₂ calculated using formula (26) and we can estimate the slope m(Q₁,Q₂) of the road passing through Q₁(x₁,y₁,z₁) and Q₂(x₂,y₂,z₂) using the following formula (see Figure 7):

$$\mathrm{m}\left({\mathrm{Q}}_1,{\mathrm{Q}}_2\right)=\tan \left(\mathsf{\beta }\right)=\frac{\left|{\mathrm{z}}_2-{\mathrm{z}}_1\right|}{\mathrm{dist}\left({\mathrm{Q}}_1,{\mathrm{Q}}_2^{\prime }\right)}=\frac{\left|{\mathrm{z}}_2-{\mathrm{z}}_1\right|}{\sqrt{{\left({\mathrm{x}}_2-{\mathrm{x}}_1\right)}^2+{\left({\mathrm{y}}_2-{\mathrm{y}}_1\right)}^2}}$$

(27)

3.3. Optimum Excavation/Filling

When a house is built on sloping land and/or with elevations and depressions, the elevation of a horizontal plane must be calculated so that the excavations over the plane would fill the space below the plane resulting in a horizontal platform on which a house can be built. In such a manner, the optimal solution between excavations and fillings is calculated, resulting in a reduction in the costs of embankments. We shall follow up with a presentation of a method to accurately calculate the elevation of this plane. The well-known “flood fill” algorithm [16] from computer graphics is adapted to deal with this problem. Since the algorithm runs in a discrete space, the coordinates will be transformed into their discrete counterparts.

We start with the level curves from the vicinity of the house (Figure 8). Then a discretization is applied (Figure 9), meaning that the plane is divided into equidistant 2D points, e.g., in Figure 9, the distance of 40 cm between points was considered. The points of the contour of the house are transformed into their discrete counterparts (for each such a point, the closest discrete point is determined). Using Bresenham’s algorithm [17], the discrete points of the contour are found (green points in Figure 9). Using a flood fill algorithm [16] in the discrete plane of the house, the discrete points inside the contour are determined (red points from Figure 9). The elevation of each red or green point is calculated using the level curves. Using these elevations, the optimal elevation of the excavation/filling plane is computed. So, the following algorithm denoted AOEF is obtained.

Algorithm 4 Optimum Excavation/Filling (AOEF):

Algorithm 4 Optimum Excavation/Filling (AOEF)

Input: the points Pi (i = 1, 2, ..., n) of the contour; for each point Pi do     Find the closed discrete point P’i; end for; for i = 1 to n − 1 do     Find the set of points Ci using Bresenham’s algorithm from point Pi to point Pi+1; end for; Using Bresenham’s algorithm from point Pn to point P1, find the set of points Cn; Find the points Qi (i = 1, 2, …, m) inside the contour given by C1 ∪ C1 ∪ … ∪ Cn; for each point Qi do     Apply A2SSTPC2L to find the elevation zi of Qi; end for;

Using the elevations z_i, calculate the elevation of the optimal excavation/filling plane (see AHEFP).

Applying the elevations z_i of the discrete points Q_i inside the contour of the house, the elevation of the optimal excavation/filling plane can be determined as follows. The elevations are sorted ascendingly. A horizontal plane π(h) (h is the height/elevation of the plane) is placed consecutively starting with the first (the lowest) elevation, continuing with the second, and ending with the last one (the highest). The following algorithm is obtained:

Algorithm 5 Height of Excavation/Filling Plane (AHEFP):

Algorithm 5 Height of Excavation/Filling Plane (AHEFP)

Input: vector of elevations z = (zi)i = 1, 2, …, m; Sort ascending the vector z; S = 0; for i = 2 to n do     S = S + zi − z1; end for; S1 = 0; S2 = S; for i = 2 to n − 1 do     S1 = S1 + (i − 1) ⋅ (zi − zi-1);     S2 = S2 − (n − i + 1) ⋅ (zi − zi-1);     if S2 ≤ S1 then      opt = i;      break;     end if; end for;

At each iteration of the last “for” loop from AFHEFP, the height of the plane is considered equal to z_i, S₁ is the sum of the distances between the elevation of z_i the plane and the elevations below the plane, and S₂ is the sum of the distances between the elevations over the plane and the elevation z_i of the plane, i.e.,:

$${\mathrm{S}}_1={\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{i}-1}\left({\mathrm{z}}_{\mathrm{i}}-{\mathrm{z}}_{\mathrm{k}}\right)$$

(28)

$${\mathrm{S}}_2={\displaystyle \sum }_{\mathrm{k}=\mathrm{i}+1}^{\mathrm{n}}\left({\mathrm{z}}_{\mathrm{k}}-{\mathrm{z}}_{\mathrm{i}}\right)$$

(29)

when starting a new iteration, S₁ is the sum of i-1 components of z and was calculated in the previous iteration. The new sum S₁ can be obtained from the previous by adding (i − 1) ⋅ (z_i − z_i-1).

At the beginning of each iteration, S₂ is the sum of n-i+1 components of z and was calculated in the previous iteration. Thus, using the value of the previous S₂, the new sum can be obtained from the previous by subtracting (n − i + 1) ⋅ (z_i − z_i-1).

Since S₁ is increasing from 0 to S and S₂ is decreasing from S to 0, the optimum is reached when S₂ becomes less or equal to S₁. More exactly, the optimum height of the plane is inside the interval [z_opt-1, z_opt]. If the values z_opt-1 and z_opt are close enough, i.e., z_opt − z_opt-1 < ε, where ε > 0 is the fixed maximum distance to the optimum, e.g., ε = 5 cm, then any of the values z_opt-1 and z_opt are good approximations of the optimum. If z_opt − z_opt-1 ≥ ε, then a divide and conquer strategy is applied to get closer to optimum. To do that, two initial planes π(z_opt-1) and π(z_opt) are considered. A new plane is placed in the middle. If the distance between S₂ calculated for π((z_opt-1+z_opt)/2) and S₁ for π(z_opt-1) is less than the distance between S₂ for π(z_opt) and S₁ for π((z_opt-1+z_opt)/2), then the optimum is further calculated in the interval [z_opt-1, (z_opt-1+z_opt)/2]. Otherwise, the optimum is calculated in the interval [(z_opt-1+z_opt)/2, z_opt] and so on.

Algorithm 6 divide and conquer for optimum excavation/filling plane (ADCOEFP):

Algorithm 6 Divide and conquer for optimum excavation/filling plane (ADCOEFP)

a = zopt-1; b = zopt; Sa,1 = S1 − (opt − 1) ⋅ (zopt − zopt-1); Sa,2 = S2 + (n – opt + 1) ⋅ (zopt − zopt-1); Sb,1 = S1; Sb,2 = S2; while b − a ≥ ε do     c = (a+b)/2;     Sc,1 = Sa,1 + (opt − 1) ⋅ (c − a);     Sc,2 = Sa,2 − (n − opt + 1) ⋅ (c − a);     if |Sc,2 − Sa,1| < |Sb,2 − Sc,1| then       a = c;       Sa,1 = Sc,1;     else       b = c;       Sb,2 = Sc,2; end if; end while; S1 = Sa,1; S2 = Sa,2;

At the end of the while loop, any of the values, a or b are good approximations of the optimum. We considered a as the optimal solution. The approximate excavation volume is V_e = S₂ ⋅ d₂, and the filling volume is V_f = S₁ ⋅ d₂. Consequently, the difference between the two volumes is:

$$\left|{\mathrm{V}}_{\mathrm{e}}-{\mathrm{V}}_{\mathrm{f}}\right|=\left|{\mathrm{S}}_2-{\mathrm{S}}_1\right|{\mathrm{d}}^2\le {\mathsf{\epsilon }\mathrm{S}}_{\mathrm{house}}\approx {\mathsf{\epsilon }\mathrm{md}}^2$$

(30)

So, the two volumes are very similar.

3.4. Numerical Example

Figure 10 shows a real topographic map for a plot of land on which a house is to be built. The topographic map of the land was drawn up by a topometric engineer, with specific high-precision equipment, called in specialized terms “station.” For small areas of land, the station is successively placed and GPS 3D coordinates are measured relative to a fixed position (in our situation it is a terminal in the Black Sea area located at a distance of approximately 717 km). The first chosen points of the topographic survey are those that describe the contour of the plot of land. In our case, the terrain contour is described by the 15 points (numbered in black). These points are inventoried and noted in a table (coordinate inventory) located on the topographic map. The following points will describe the contour of the level elevations, the positioning of the access road, the position of the neighboring houses (they are marked on the map in red together with the elevation), etc. In the case of lands with a high slope, such as the one exemplified, the measured points must be thickened for a greater accuracy of tracing the level curves (in total we have over 100 such points).

The designer receives in electronic format a topographic map similar to the one in Figure 10. On this topographic support, the designer will place the designed building on the scale, having to respect a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis as in Figure 11. Once the construction is located, the new contour will appear on the topographic map described by a few new points (minimum 4, in the example, for the description of the contour, 6 points are used). A situation plan is drawn up in which, according to the legislation, the 2D coordinates of the corners of the construction must be specified (see

Table 1. The x and y coordinates of the Q points.

Points	E (m)	N (m)
Q1	552,749.248	456,578.380
Q2	552,754.270	456,571.835
Q3	552,766.210	456,580.997
Q4	552,760.244	456,588.772
Q5	552,754.733	456,584.542
Q6	552,753.385	456,581.554

) together with the elevation of the horizontal excavation/filling plane. In Figure 11, this elevation is denoted by CTA and is estimated at 707 m by the designer. A good determination of this elevation is not easy, especially if the terrain has a high slope.

A detail of Figure 11 is presented in Figure 12 illustrating the points ${\mathrm{A}}_{11}^2,{\mathrm{A}}_{12}^2,{\mathrm{A}}_{21}^2,{\mathrm{and}\mathrm{A}}_{22}^2$ found using the algorithm AIPS.

Table 2. 3D coordinates of the Q points.

Points	E (m)	N (m)	Z(m)	Calculated Z(m)
Q1	552,749.2478	456,578.3797		707.8851
${\mathrm{A}}_{11}^1$	552,749.7024	456,578.7371	707.6000
${\mathrm{A}}_{12}^1$	552,749.4725	456,578.9520	707.6000
${\mathrm{A}}_{21}^1$	552,748.8265	456,578.0230	708.0000
${\mathrm{A}}_{22}^1$	552,749.2769	456,577.9832	708.0000
Q2	552,754.2700	456.,571.8350		709.3279
${\mathrm{A}}_{11}^2$	552,754.0980	456,572.9276	709.2000
${\mathrm{A}}_{12}^2$	552,754.4526	456,572.8445	709.2000
${\mathrm{A}}_{21}^2$	552,754.0176	456,571.5626	709.6000
${\mathrm{A}}_{22}^2$	552,754.3366	456,571.5767	709.6000
Q3	552,766.2100	456,580.9970		704.2017
${\mathrm{A}}_{11}^3$	552,768.2179	456,583.3340	704.0000
${\mathrm{A}}_{12}^3$	552,768.7357	456,582.6679	704.0000
${\mathrm{A}}_{21}^3$	552,760.8755	456,578.3463	704.4000
${\mathrm{A}}_{22}^3$	552,760.9848	456,577.8984	704.4000
Q4	552,760.2440	456,588.7720		704.2330
${\mathrm{A}}_{11}^4$	552,762.5489	456,590.7208	704.0000
${\mathrm{A}}_{12}^4$	552,762.9178	456,590.3418	704.0000
${\mathrm{A}}_{21}^4$	552,755.2345	456,585.1132	704.4000
${\mathrm{A}}_{22}^4$	552,754.9360	456,585.3872	704.4000
Q5	552,754.7330	456,584.5420		704.7360
${\mathrm{A}}_{11}^5$	552,755.2345	456,585.1132	704.4000
${\mathrm{A}}_{12}^5$	552,755.1438	456,584.4878	704.4000
${\mathrm{A}}_{21}^5$	552,752.7819	456,582.6377	704.8000
${\mathrm{A}}_{22}^5$	552,752.3362	456,582.9349	704.8000
Q6	552,753.3850	456,581.5540		705.0855
${\mathrm{A}}_{11}^6$	552,753.3671	456,581.7908	704.8000
${\mathrm{A}}_{12}^6$	552,753.7000	456,581.5900	704.8000
${\mathrm{A}}_{21}^6$	552,753.2125	456,581.6078	705.2000
${\mathrm{A}}_{22}^6$	552,753.5536	456,581.3981	705.2000

presents the coordinates of the points ${\mathrm{A}}_{11}^{\mathrm{k}},{\mathrm{A}}_{12}^{\mathrm{k}},{\mathrm{A}}_{21}^{\mathrm{k}},{\mathrm{and}\mathrm{A}}_{22}^{\mathrm{k}}$ calculated for each point Q_k (k = 1, 2, …, 6) by applying AIPS. Using A2SSTPC2L and (26), the elevation for each point Q_k (k = 1, 2, …, 6) is calculated.

Usually, for the dimensions of the project to be in accordance with the real construction, before commencing all the construction works, a plot is made on the field of the future building. For this purpose, the ground tracing is performed. The corners of the house are very precisely marked (the Q points). In our example, the provided measured elevations of the 6 points are presented in

Table 3. The estimated and the measured z-coordinates of the Q points.

Points	Estimated Z (m)	Measured Z (m)
Q1	707.8921	707.8990
Q2	709.3279	709.3260
Q3	704.2017	704.2140
Q4	704.2330	704.2210
Q5	704.7360	704.7500
Q6	705.0855	705.0970

.

As seen in Table 3, the estimated elevation and the measured one for each Q point are very close. The difference between them is less than 2 cm.

Using the coordinates of the points Q_k (k = 1, 2, …, 6) from Table 1, Bresenham, and, then flood fill algorithms were applied resulting in 803 points located inside the perimeter of the house. The elevations of these points were calculated using A2SSTPC2L and (26). The algorithm AHEFP was applied for these 803 elevations and the height of the horizontal excavation/filling plane was accurately estimated as 706.2571 m, instead of the value 707 m proposed by the designer. Out of the 803 points, 412 were located below the horizontal plane, and 391 had the elevation over this plane.

4. Discussion

In order to compare the two proposed methods (from Section 2.1 and Section 2.2) for determining the points P₁ and P₂ 100 million experiments were performed. For each experiment, two random lines d₁ and d₂ were selected and the point Q was also randomly chosen. The tests were performed using an ASUS ROG GL752VW-T4015D laptop with Intel^® Core™ i7-6700HQ 2.60 GHz processor, and 8 GB of RAM.

First, we analyzed the speed of the proposed methods. The total running time of the first method was 5.82 s. The total running time of the second method was 40.14 s. So, the first method is almost 7 times faster than the second one (see Figure 13). The first method is considerably faster because it requires elementary calculations, while the second one runs an iterative algorithm to solve the Equation (18).

We also calculated the accuracy of the first method against the second method since the error for the second method can be set as low as desired. An average error of 0.61% was obtained, the lowest error being 0 and the highest 1.17% (

Table 4. Errors of the first method obtained from 100 million experiments.

	Errors (First Method against Second Method)
Lowest error	0%
Average error	0.61%
Highest error	1.17%

). So, we concluded that the accuracy of the approximate first method is good enough.

The steps of positioning the house will be next discussed, and the utility of the methods proposed in this paper will be shown. The designer receives a topographic map in electronic format. On this topographic support, the designer will place the designed building on the scale, having to comply with a series of rigors imposed by law, among which: the minimum distance from the property limit and/or from the road axis, etc. Once the construction is located, the outline of the house described by several new points will appear on the topographic map (minimum 4). At these new points, it is necessary to specify: the coordinates on the three dimensions according to which the elevation is specified, and the distances of the construction location from the property limits. The determination of the coordinates in the plan does not raise any concerns, but there are difficulties in obtaining the real vertical elevation especially if the terrain has a high slope. The article proposes a method for determining the exact elevation of the land at the corners of the building, and in fact for the entire contour. Depending on the contour dimensions, the designer establishes a position of the horizontal plane, a plane that corresponds to the finished floor of the ground level. The decision to establish this quota should be contingent on the volumes of embankments that will be made. The optimum is obtained when the excavation volume is approximately equal to the filling volume. The article solves this problem in Section 3.3.

For the design of a road, a primary route of the road is made on a topographic map and road recognition is performed. A final design of the road follows both along its entire length and on its cross-sections. Figure 14 shows with a red line the optimal trajectory of a road depending on the elevation of the natural land. Because the considered road is a highway, the rules impose stricter restrictions on the slopes. If the allure of the natural terrain in certain areas is followed, the maximum allowed inclination of the road will inevitably be exceeded. As it can be seen in Figure 14, in these areas, measures must be taken to straighten the slope, such as fillings, excavations or the provision of bridges and viaducts, tunnels respectively. Through detailed awareness of the elevations and the inclination of the land in the points near the road, an optimal vertical tracing of the road can be established so that we have minimum excavations and fillings for an imposed road slope, or we can optimize the lengths of the bridges, respectively tunnels [12].

When designing a new road or railway, a number of initial design factors must be taken into account, such as: road category, design speed, width of a traffic lane, number of lanes, traffic frequency, equipment size, etc. Depending on the route imposed on the road, there are also other specific aspects that must be taken into account:

• Avoid high erosion hazard sites, particularly where mass failure is a possibility.
• Utilize natural terrain features such as stable benches, ridgetops, and low gradient slopes to minimize the area of road disturbance.
• If necessary, include short road segments with steeper gradients to avoid problem areas or to utilize natural terrain features.
• Avoid midslope locations on long, steep, or unstable slopes.
• Locate roads on well-drained soils and rock formations which dip into slopes rather than areas characterized by seeps, highly plastic clays, concave slopes hummocky topography, cracked soil and rock strata dipping parallel to the slope.
• For logging road, utilize natural log landing areas (flatter, benched, well-drained land) to reduce soil disturbance associated with log landings and skid roads.
• Avoid undercutting unstable, moist toe slopes when locating roads in or near a valley bottom.
• Roll or vary road grades where possible to dissipate flow in road drainage ditches and culverts and to reduce surface erosion.
• Select drainage crossings to minimize channel disturbance during construction and to minimize approach cuts and fills.
• Locate roads far enough above streams to provide an adequate buffer, or provide structure or objects to intercept sediment moving down slope below the road.
• If an unstable area such as a headwall must be crossed, consider end hauling excavated material rather than using sidecast methods. Avoid deep fills and compact all fills to accepted engineering standards. Design for close culvert and cross drain spacing to effectively remove water from ditches and provide for adequate energy dissipators below culvert outlets. Horizontal drains or interceptor drains may be necessary to drain excess groundwater [13].

As seen in Figure 15, the angle of inclination of the terrain is an important parameter to know in order to avoid erosion and then landslides [18,19].

From the point of view of the cross-sections, the position of the road in relation to the terrain slope must take into account the above constraints. By the same token, knowing the slope of the land is essential for establishing the position both horizontally and especially vertically of the road, and for achieving a low cost related to the embankment works. This position is established according to the Full Bench Road Prism, which contains the elements in Figure 16.

For very steep road areas such as the one in Figure 17, knowing the slope is imperative because substantial savings in embankment volumes can be made, resulting in decreasing the quantities of gabions or excavations [13].

The problem related to determining the slope of natural land at a point on the topographic map is necessary both in the longitudinal section and in the cross-section of roads or railways and have the same solution described in Section 3.1 and Section 3.2.

The proposed methods can also be applied in environment protection, and hydrotechnical engineering [18,19]. The evaluation of the slope is also necessary for river levees. Slope stability analyses are conducted with rising water levels until certain failure is reached. Discharge occurs at a water height of around 7 cm above the crown, which means the slope would actually fail before reaching the hydraulic heads used for slope stability calculation during the overflow [20].

5. Conclusions

To conclude, the correct estimation of elevation at points has a practical application to calculate the optimum excavation/filling plane when a new house is built. The article presented a method to estimate accurately the elevation of this plane which reduces the costs of embankments.

When designing a new road or railway, designers need to know the ground slope at some points of the future road or rail as exactly as possible, as well as the elevation of each considered point, and, also, the slope of each consecutive two points. Instead of sending workers on the field to perform measurements, these values can be calculated with very good accuracy using the methods described in Section 2.1, and Section 2.2. Two algorithms to solve SSTPC2L are presented, a fast and approximate one and an exact but slower one. The error of the first method is low (see Table 1). It has the advantage of being very fast and accurate enough. Thus, if there are many points for which SSTPC2L is applied, the first method can be preferred. However, if we look for a more accurate solution, the second method is more suitable.