Terrain Analysis in Google Earth Engine: A Method Adapted for High-Performance Global-Scale Analysis

Abstract

Terrain analysis is an important tool for modeling environmental systems. Aiming to use the cloud-based computing capabilities of Google Earth Engine (GEE), we customized an algorithm for calculating terrain attributes, such as slope, aspect, and curvatures, for different resolution and geographical extents. The calculation method is based on geometry and elevation values estimated within a 3 × 3 spheroidal window, and it does not rely on projected elevation data. Thus, partial derivatives of terrain are calculated considering the great circle distances of reference nodes of the topographic surface. The algorithm was developed using the JavaScript programming interface of the online code editor of GEE and can be loaded as a custom package. The algorithm also provides an additional feature for making the visualization of terrain maps with a dynamic legend scale, which is useful for mapping different extents: from local to global. We compared the consistency of the proposed method with an available but limited terrain analysis tool of GEE, which resulted in a correlation of 0.89 and 0.96 for aspect and slope over a near-global scale, respectively. In addition to this, we compared the slope, aspect, horizontal, and vertical curvature of a reference site (Mount Ararat) to their equivalent attributes estimated on the System for Automated Geospatial Analysis (SAGA), which achieved a correlation between 0.96 and 0.98. The visual correspondence of TAGEE and SAGA confirms its potential for terrain analysis. The proposed algorithm can be useful for making terrain analysis scalable and adapted to customized needs, benefiting from the high-performance interface of GEE.

1. Introduction

Terrain analysis is essential for modeling environmental systems [1,2,3]. The variability of landforms is frequently used to understand, map or model geomorphological, hydrological, and biological processes [4,5,6,7]. Elevation has a strong relationship with terrestrial temperature, vegetation type, and with the potential energy accumulated on a slope. The aspect and derived products, such as Northernness and Easternness attributes, can be linked to the potential solar irradiation on terrain. The Slope gradient, for example, controls the overland and subsurface flow velocity and runoff rate. Similarly, curvatures are associated with acceleration and dispersion of water and sediment flows, which impacts the erosion and soil water content [8].

The public availability of elevation data with global coverage, such as the digital elevation model (DEM) derived from NASA’s Shuttle Radar Topography Mission (SRTM DEM, [9]) and the digital surface model from the Advanced Land Observing Satellite (AW3D30 DSM, [10]), has promoted the exploration of topographic features in different contexts using processing tools available in several geographic information systems (GIS) [4,11,12]. However, despite the popularization of many global elevation datasets, it is important to pay attention to their quality when used for modelling purposes, as the acquisition mean and other production aspects can significantly impact the outputs [13,14]. In addition, analyzing big geospatial datasets can pose some limitations to traditional GIS. This becomes more critical with the availability of new digital datasets, which are providing better temporal and spatial resolutions due to advances in sensor technologies [15].

The Global Multi-resolution Terrain Elevation Data 2010 [16] and the global suit of terrain attributes [2] are examples of datasets that were produced using large computational tasks for mapping the global extent and in different spatial resolutions, which demanded optimized processing architectures. In general, high performance architectures are based on splitting the data in smaller subsets (tiles) to take the advantage of distributed computing operations. Recently, with the advent and popularization of cloud-based interfaces for processing big geospatial data, e.g., Google Earth Engine [17], the Pangeo software packages [18], and Actinia REST service [19], computational tasks applied to terrain analysis could be scaled and customized directly by the user.

Earth Engine (GEE) is a cloud-based platform developed by Google that supports the global-scale analysis of big catalogs of Earth Observation data [17]. It has been used to map global forest change in the 21st century [20], Earth’s surface water change [21], global urban areas [11], wildfire progression [22], global bare surface change [23], and others. In this sense, GEE becomes compelling not because the distributed processing tasks are executed on the server-side of Google, but also due to the increasing availability of many global geospatial datasets that could be explored in topographic mapping. There exist several available topographical data within GEE, such as the global SRTM DEM, AW3D30 DSM, Global 30 Arc-Second Elevation data (GTOPO30 DEM, [24]), and others. Thus, GEE characteristics could permit the customization of high-performance terrain analysis with minimal user input and any computational processing on the user side. In fact, GEE provides three algorithms for calculating slope, illumination, and aspect of terrain, but lacks in providing calculation methods of other terrain information, such as the curvatures and landscape characterization.

In addition, a common obstacle of global terrain analysis in common GIS is the need for projecting DEMs onto projected coordinate systems, which ensures the elevation data is equally spaced on a plane square grid [25]. This step is complicated because it is difficult to define a projected system that minimizes terrain distortions over a global extent [26]. Moreover, as many available global DEMs are referenced by geographical coordinate systems and some researchers continue to apply square-grid algorithms to them, the algorithms should consider the geometry and specificity of global spheroidal DEMs [25]. This aspect is important because the application of square-grid methods to spheroidal equal angular DEMs leads to substantial computational errors in models of morphometric variables [25].

In this paper, we aimed at describing and making available a user-friendly processing algorithm for performing terrain analysis in GEE. This algorithm takes advantage of GEE’s high-performance architecture for making the computational analysis scalable, adapted to customized needs, and requiring minimal user input. For this, the proposed package takes advantage of a calculation method adapted for spheroidal elevation grids, which favors the global-scale analysis of different DEM resolutions without projecting elevation data.

2. Material and Methods

2.1. Algorithm Description

The Terrain Analysis in GEE (TAGEE) package use calculation methods adapted to spheroidal angular grids, i.e., the DEM can be referenced in a geographical coordinate system, e.g., the World Geodetic System (WGS84). The following paragraphs briefly describe the calculation methods performed by TAGEE package. The readers are referred to [8] for the mathematical concepts of geomorphometry, a historical overview of the progress of digital terrain modelling, and the notion of the topographic surface and its limitations.

2.1.1. Topographic Surface

The land topography can be approximated by a topographic surface defined by a continuous, single-valued bivariate function (Equation (1)) [8]:

$$z=f\left(x,y\right)$$

(1)

where $z$ is elevation (meters), and $x$ and $y$ are the coordinates in geographical coordinates (degrees).

The local morphometric variables are functions of the partial derivatives of elevation. Using the Evans–Young method, the function $z=f\left(x,y\right)$ is expressed as the second-order bivariate Taylor polynomial (Equation (2)):

$$z=\frac{r{x}^2}{2}+\frac{t{y}^2}{2}+sxy+px+qy+u$$

(2)

where $r$, $t$, $s$, $p$, and $q$ are the partial derivatives, and $u$ is the residual term.

Differently from a digital elevation model projected on a plane square grid, where the partial derivatives of terrain are estimated by finite differences, the processing and analysis of a spheroidal equal angular DEM must consider the spheroidal geometry. In such case, a grid spacing with approximately equal linear units along meridians and parallels exists only at the equator. To estimate the parameters of a spheroidal grid, a 3 × 3 moving window must retrieve both the geometry elements and the elevation values of the window nodes (Figure 1).

2.1.2. Terrain Parameters: Neighbor Elevations and Geometries

The elevation values of a 3 × 3 moving window are estimated by convolution kernels. For geometries, the Haversine formula is used to determine the great-circle distances between two neighbor nodes within the spheroidal window, given their latitude and longitude geographical positions (Equations (3)–(5)):

$$j={\sin }^2\left(\frac{\mathsf{\Delta }\varphi }{2}\right)+\cos {\varphi }_1\cdot \cos {\varphi }_2\cdot {\sin }^2\left(\frac{\mathsf{\Delta }\lambda }{2}\right)$$

(3)

$$k=2\cdot atan2\left(\sqrt{j},\sqrt{\left(1-j\right)}\right)$$

(4)

$$l=R\cdot k$$

(5)

where ${\varphi }_1$ is latitude for the first given node in radians, ${\varphi }_2$ is the latitude for the second given node in radians, ${\lambda }_1$ is the longitude for the first given node in radians, ${\lambda }_2$ is the longitude for the second given node in radians, $\mathsf{\Delta }\varphi$ and $\mathsf{\Delta }\lambda$ are the respective differences of latitude and longitude between the given nodes, and $R$ is the mean radius of Earth equals to 6,371,000 meters. The linear distance $l$ is given in meters.

Knowing the latitude and longitude of the window nodes (Figure 1), the Haversine formula allows the calculation of linear distances of a, b, c, d, and e, which are used with the neighbor elevation values (from $z_1$ to $z_9$) to calculate the partial derivatives of terrain.

2.1.3. Terrain Derivatives

To estimate the first and second-order partial derivatives r, t, s, p and q, the polynomial model is fitted by least squares and results in the following estimations (Equations (6)–(10)) [8]:

$$p=\frac{a^{2}cd\left(d+e\right)\left(z_3-z_1\right)+b\left(a^2{d}^2+c^2{e}^2\right)\left(z_6-z_4\right)+a{c}^{2}e\left(d+e\right)\left(z_9-z_7\right)}{2[a^2{c}^2\left(d+e{)}^2+b^2\left(a^2{d}^2+c^2{e}^2\right)\right]}$$

(6)

$$\begin{array}{ll}\hfill q=& \frac{1}{3de\left(d+e\right)\left(a^4+b^4+c^4\right)}\hfill \\ \hfill & \cdot \{\left[d^2\left(a^4+b^4+b^2{c}^2\right)+c^2{e}^2\left(a^2-b^2\right)\right]\left(z_1+z_3\right)\hfill \\ \hfill & -\left[d^2\left(a^4+c^4+b^2{c}^2\right)-e^2\left(a^4+c^4+a^2{b}^2\right)\right]\left(z_4+z_6\right)\hfill \\ \hfill & -\left[e^2\left(b^4+c^4+a^2{b}^2\right)-a^2{d}^2\left(b^2-c^2\right)\right]\left(z_7+z_9\right)\hfill \\ \hfill & +d^2\left[b^4\left(z_2-3z_5\right)+c^4\left(3z_2-z_5\right)+\left(a^4-2b^2{c}^2\right)\left(z_2-z_5\right)\right]\hfill \\ \hfill & +e^2\left[a^4\left(z_5-3z_8\right)+b^4\left(3z_5-z_8\right)+\left(c^4-2a^2{b}^2\right)\left(z_5-z_8\right)\right]\hfill \\ \hfill & -2\left[a^2{d}^2\left(b^2-c^2\right)z_8+c^2{e}^2\left(a^2-b^2\right)z_2\right]\}.\hfill \end{array}$$

(7)

$$r=\frac{c^2\left(z_1+z_3-2z_2\right)+b^2\left(z_4+z_6-2z_5\right)+a^2\left(z_7+z_9-2z_8\right)}{a^4+b^4+c^4}$$

(8)

$$\begin{array}{ll}\hfill s=& \{c\left(a{[}^2\left(d+e\right)+b^{2}e\right]\left(z_3-z_1\right)\hfill \\ \hfill & -b\left(a^{2}d-c^{2}e\right)\left(z_4-z_6\right)\hfill \\ \hfill & +a\left[c^2\left(d+e\right)+b^{2}d\right]\left(z_7-z_9\right)\}\hfill \\ \hfill & \cdot \frac{1}{2[a^2{c}^2\left(d+e{)}^2+b^2\left(a^2{d}^2+c^2{e}^2\right)\right]}\hfill \end{array}$$

(9)

$$\begin{array}{ll}\hfill t=& \frac{2}{3de\left(d+e\right)\left(a^4+b^4+c^4\right)}\hfill \\ \hfill & \cdot \{\left[d\left(a^4+b^4+b^2{c}^2\right)-c^{2}e\left(a^2-b^2\right)\right]\left(z_1+z_3\right)\hfill \\ \hfill & -\left[d\left(a^4+c^4+b^2{c}^2\right)+e\left(a^4+c^4+a^2{b}^2\right)\right]\left(z_4+z_6\right)\hfill \\ \hfill & +\left[e\left(b^4+c^4+a^2{b}^2\right)+a^{2}d\left(b^2-c^2\right)\right]\left(z_7+z_9\right)\hfill \\ \hfill & +d\left[b^4\left(z_2-3z_5\right)+c^4\left(3z_2-z5\right)+\left(a^4-2b^2{c}^2\right)\left(z_2-z_5\right)\right]\hfill \\ \hfill & +e[a^4\left(3z_8-z_5\right)+b^4\left(z_8-3z_5\right)+c^4-2a^2{b}^2\left]\left(z_8-z_5\right)\right]\hfill \\ \hfill & -2\left[a^{2}d\left(b^2-c^2\right)z_8-c^{2}e\left(a^2-b^2\right)z_2\right]\}\hfill \end{array}$$

(10)

where the parameters $a$, $b$, $c$, $d$, and $e$ are the linear distances calculated from the Haversine formula (Equations (3)–(5)), and the $z$ values are elevation values from the neighbors of a moving window (Figure 1).

2.1.4. Terrain Attributes

Local attributes, such as slope, aspect, and curvatures, are calculated from the partial derivatives of terrain [8]. The slope gradient ($G$, Equation (11)) is a flow attribute that relates to the velocity of gravity-driven flows. For measuring the direction, the slope aspect is used ($A$, Equations (12) and (13)). Additionally, one can calculate the amount that a slope is faced to the North or East, resulting in the Northernness ($A_N$, Equation (14)) and Easternness ($A_E$, Equation (15)) derived from the aspect. The remaining flux attributes that can be calculated from the first and second-order partial derivatives are the horizontal ($k_h$, Equation (16)) and vertical curvatures ($k_v$, Equation (17)). While the horizontal curvature relates if a lateral flow converges ($k_h$ < 0) or diverges ($k_h$ > 0), the vertical curvature measures the relative acceleration ($k_v$ > 0) and deceleration ($k_v$ < 0) of a gravity-driven flow:

$$G=\arctan \sqrt{p^2+q^2}$$

(11)

$$\begin{array}{ll}\hfill A=& -90\left[1-\mathrm{sign}\left(q\right)\right]\left(1-\left|\mathrm{sign}\left(p\right)\right|\right)+180\left[1+\mathrm{sign}\left(p\right)\right]\hfill \\ \hfill & -\frac{180}{\pi }\mathrm{sign}\left(p\right)\arccos \left(\frac{-q}{\sqrt{p^2+q^2}}\right)\hfill \end{array}$$

(12)

$$\mathrm{sign}\left(x\right)=\left\{\begin{array}{ll}1\hfill & \mathrm{for}x>0\hfill \\ 0\hfill & \mathrm{for}x=0\hfill \\ -1\hfill & \mathrm{for}x<0\hfill \end{array}\right.$$

(13)

$$A_N=\cos A$$

(14)

$$A_E=\sin A$$

(15)

$$k_h=-\frac{q^{2}r-2pqs+p^{2}t}{\left(p^2+q^2\right)\sqrt{1+p^2+q^2}}$$

(16)

$$k_v=-\frac{p^{2}r+2pqs+q^{2}t}{\left(p^2+q^2\right)\sqrt{{(1+p^2+q^2)}^3}}$$

(17)

Differently from flow attributes, which are gravity field-specific variables, form attributes are related to principal sections of terrain [8]. The mean curvature ($H$, Equation (18)) is a half-sum of any two orthogonal normal sections and represents two accumulation mechanisms of gravity-driven flows with equal weights: convergence and relative deceleration. Among the class of form attributes, the Gaussian curvature ($K$, Equation (19)) is a product of maximal ($k_{max}$) and minimal ($k_{min}$) curvatures. The two principal curvatures calculate the highest and lowest curvature for a given point of the topographic surface. The maximal curvature ($k_{max}$, Equation (20)) is useful for mapping rigdes ($k_{max}$ > 0) and closed depressions ($k_{max}$ < 0). Likewise, the minimal curvature ($k_{min}$, Equation (21)) is useful for identifying hills ($k_{min}$ > 0) and valleys ($k_{min}$ < 0) across the topographic surface. With the results of mean and Gaussian curvatures, a landform classification can be generated after [27] proposing the continuous form of the Gaussian classification [8,28]. Instead of providing categorical values, the shape index ($SI$, Equation (22)) ranges from –1 to 1 and map convex ($SI$ > 0) and concave ($SI$ < 0) landforms:

$$H=-\frac{\left(1+q^2\right)r-2pqs+\left(1+p^2\right)t}{2\sqrt{{(1+p^2+q^2)}^3}}$$

(18)

$$K=\frac{rt-s^2}{{(1+p^2+q^2)}^2}$$

(19)

$$k_{max}=H+\sqrt{\left(H^2-K\right)}$$

(20)

$$k_{min}=H-\sqrt{\left(H^2-K\right)}$$

(21)

$$SI=\frac{2}{\pi }\arctan \frac{H}{\sqrt{H^2-K}}$$

(22)

2.2. Package Description

Calculation methods presented in this paper were developed using the JavaScript programming interface available as the online code editor of GEE. TAGEE was developed by different modules of calculation, similarly to what was described in Methods. The first module, calculateParameters, uses convolution kernels and the Haversine formula to retrieve elevation values and the spheroidal geometries of a 3 × 3 moving window. In this module, a digital elevation model and a square polygon representing the bounding box (min. Longitude, min. Latitude, max. Longitude, and max. Latitude, in the WGS84 coordinate reference system) are required as input parameters to run. The bounding box is used both in this module and others for generating images with constant values and restrict the calculations to the study area. The first module returns an image with 14 bands, i.e., the neighbor elevation values (from $Z_1$ to $Z_9$) and the distances ($a$, $b$, $c$, $d,$ and $e$) (Figure 1).

Once the basic parameters (elevation and distances) were established, the partial derivatives of terrain are calculated with the calculateDerivatives module. This second module requires the returned parameters from calculateParameters and also the bounding box of the study region. The second module adds the partial derivatives ($r$, $t$, $s$, $p,$ and $q$) as new bands to the previous image. Then, terrain attributes are calculated by the module calculateAttributes (Figure 2).

Terrain attributes can also be calculated by a single function, without calling the intermediate modules. The final output, for both alternatives (Figure 2), is a multi band object containing the same data properties of the digital elevation model (resolution, data type, and coordinate reference system) with 13 bands (

Table 1. Attributes of terrain, with their units and description, calculated by TAGEE package.

Attribute	Unit	Description
Elevation	meter	Height of terrain above sea level
Slope	degree	Slope gradient
Aspect	degree	Compass direction
Hillshade	dimensionless	Brightness of the illuminated terrain
Northernness	dimensionless	Degree of orientation to North
Easternness	dimensionless	Degree of orientation to East
Horizontal curvature	meter	Curvature tangent to the contour line
Vertical curvature	meter	Curvature tangent to the slope line
Mean curvature	meter	Half-sum of the two orthogonal curvatures
Minimal curvature	meter	Lowest value of curvature
Maximal curvature	meter	Highest value of curvature
Gaussian curvature	meter	Product of maximal and minimal curvatures
Shape Index	dimensionless	Continuous form of the Gaussian classification

). The final attributes can be used for further modeling inside GEE or thematic mapping.

The package has an additional feature that makes easier the visualization of terrain attributes. As the range of attribute values and the pixel resolution may vary according to the visualization level (zoom), which impacts the estimated geometries and elevation neighbor values, a module called makeVisualization automatically calculates the dynamic legend defined by the 0.05 and 0.95 percentiles within the bounding box. In addition, different color palettes for making the map legend are available in TAGEE: rainbow, inferno, cubehelix, red2green, green2red, elevation, and aspect. The package code and a minimal reproducible example are available in https://github.com/zecojls/tagee (Supplementary Materials).

2.3. Statistical Evaluation

We performed the evaluation of TAGEE attributes by comparing the aspect and slope derived from two available functions of GEE (ee.Terrain.aspect and ee.Terrain.slope) on a near-global scale. For this task, we used the Pearson correlation analysis with the SRTM DEM 30m, which contains elevation in meters limited to an area between about 60° north latitude and 56° south latitude. It is important to mention that for the currently available terrain functions of GEE, the local gradient is computed using the four-connected neighbors of each pixel, differently from the proposed method of TAGEE, which uses a 3 × 3 pixel window and also considers the spheroidal geometries in its calculation. Thus, minimal differences between the calculation methods are expected to occur. This analysis was performed in GEE and, in addition to Pearson’s correlation, we calculated the relative mean absolute error (MAE) between the outputs. The relative MAE is estimated by calculating the mean absolute difference between two rasters and standardizing the result to the range (maximum minus minimum values) of the reference raster.

Similarly, we compared the results from TAGEE with terrain attributes calculated by the System for Automated Geoscientific Analyses (SAGA) GIS version 2.3.2 [12]. In this case, we downloaded from GEE the 30 m SRTM DEM together with the resulting 12 attributes calculated by TAGEE, all covering the Mount Ararat (located between 44.2° and 44.5° E, and 39.6° and 39.8° N). The Mount Ararat was selected due its high variability of landforms and the availability of published maps from previous works [8,29], allowing the visual comparison of spatial patterns. The Mount Ararat SRTM-DEM was processed in SAGA GIS using the “Slope, Aspect, Curvature” from the Morphometry module of Terrain Analysis. The calculation method was the “Evan (1979)” based on six parameters and 2nd order polynomials, similarly to the TAGEE calculation method. The comparison was performed by calculating the Pearson’s correlation coefficient (r) and the relative MAE, where the aspect, slope, horizontal curvature and vertical curvature from TAGEE were compared with aspect, slope, tangential, and profile curvature from SAGA GIS, respectively, following the equivalence described in [8].

3. Results and Discussion

The statistical analysis revealed a significant correlation (p < 0.01) of the TAGEE outputs with equivalent terrain attributes calculated from GEE and SAGA GIS (

Table 2. Comparison of TAGEE attributes with outputs from GEE and SAGA GIS algorithms.

Attribute	Region	Reference	Pearson’s r	rMAE 1
Aspect	Near global SRTM DEM 30 m	GEE	0.89 *	13%
Slope	Near global SRTM DEM 30 m	GEE	0.98 *	2%
Aspect	Mount Ararat SRTM DEM 30 m	SAGA GIS	0.96 *	4%
Slope	Mount Ararat SRTM DEM 30 m	SAGA GIS	0.98 *	3%
Horizontal curvature	Mount Ararat SRTM DEM 30 m	SAGA GIS	0.98 *	4%
Vertical curvature	Mount Ararat SRTM DEM 30 m	SAGA GIS	0.98 *	4%

* Significant for p < 0.01; ¹ relative mean absolute error.

). The slope estimated over a near-global extent reached a correlation of 0.98 (error of 2%) between TAGEE and functions of GEE, while the aspect resulted in a Pearson’s r of 0.89 (13% of error). The lower correlation of aspect can be associated to its dimension nature, i.e., a circular variable, as well as to the differences of calculation methods between TAGEE and GEE. Despite the small differences, TAGEE revealed the same spatial patterns and allowed the estimation of additional attributes at the global scale, such as the Northernness, horizontal and vertical curvatures (Figure 3A–C, respectively). The main mountain ranges of the Earth, such as the Rocky Mountains in North America, Andes in South America, Alps in Europe, Himalayas, and Tibetan plateau in Asia, etc., present the highest curvatures calculated by TAGEE. Conversely, the plains and flat surfaces had the lowest estimates for both curvatures. The degree of orientation to North (Figure 3A) also depict the main landforms of the Earth.

TAGEE was developed in GEE to take advantage of the high-performance computing of the platform. As the cloud-based interfaces have created much enthusiasm and engagement in the remote sensing and geospatial fields, many processing algorithms have been adapted to make substantive progress on global challenges involving the processing of big geospatial data [30]. In this sense, GEE is providing petabytes of publicly available remote sensing imagery and other ready-to-use products. The high-speed parallel processing of GEE servers and the libraries of operators and machine learning algorithms available by Application Programming Interfaces (APIs) in popular coding languages, such as JavaScript and Python, are enabling users to discover, analyze, and visualize geospatial big data without needing access to supercomputers [30]. Within this framework, TAGEE supports the development of customized terrain analysis with different elevation data across large geographical extents.

When TAGEE outputs were compared to those from SAGA GIS (Table 2), the statistical evaluation resulted in a significant and high correlation for the slope, horizontal and vertical curvatures of terrain (Pearson’s r of 0.98, with an error difference of 3 and 4%). Aspects from TAGEE and SAGA GIS had an inferior correlation coefficient, but the result was higher than the aspect from the algorithm of GEE. The region of Mount Ararat was also used to visually compare the slope, horizontal and vertical curvatures, calculated from both TAGEE and SAGA GIS (Figure 4). The 3D visualizations revealed a high similarity between both maps, but some small differences can be visualized by the color intensity. This is the case of the slope of the Mount Ararat calculated by TAGEE (Figure 4A), which had a higher intensity compared to the slope of SAGA GIS (Figure 4B). A slightly higher intensity for the vertical curvature calculated by SAGA GIS was also evident on an edge of the Mount Ararat (Figure 4F). Despite being small, these visual differences confirm the relative error of both methods (Table 2). In addition, the spatial patterns of aspect, slope, and curvatures from TAGEE presented a high correspondence with the terrain maps of Mount Ararat available in [8,29], reinforcing the confidence of the TAGEE calculation method.

In this work, the TAGEE algorithm was developed to consider spheroidal geometries in its calculation method. This approach diverges from the techniques available in traditional GIS, where TAGEE considers the great circle distances of the DEM defined by Latitude and Longitude positions. Common GIS software, such as SAGA GIS, requires the projection of the DEM to ensure the elevation data have the same pixel size. However, as identified by [25], some researchers continue to apply square-grid algorithms to spheroidal equal angular DEMs, which can lead to substantial computational errors in models of morphometric variables. The small relative errors between TAGEE and GEE or SAGA GIS could be linked to the differences in their calculation methods.

Finally, some limitations of TAGEE can also be noted. Only local morphometric variables can be calculated by the package, which includes flux and form attributes. Non-local attributes, such as specific catchment area, were not implemented due to the absence of a general analytical theory, which is still little developed [29], and due to the recursion processing that is still challenging within GEE [17]. Furthermore, a novel method became available to handle major problems of terrain analysis, which includes the approximation of DEM, generalization and denoising, and the computation of morphometric variables. The universal spectral analytical method based on high-order orthogonal expansions using the Chebyshev polynomials were developed by [31] to handle the aforementioned issues into an integrated framework, but was not implemented in this work.

4. Conclusions

The proposed package (TAGEE) can calculate terrain attributes using the high-performance platform of GEE with an accuracy equivalent to traditional GIS. The approach of using spheroidal geometries does not require the projection of input elevation data for terrain attributes calculation. The comparison between algorithms demonstrated that TAGEE estimates terrain slope and aspect similarly to the available functions of GEE. The advantage of TAGEE over the currently available functions is that additional outputs can be produced, such as curvatures and shape index, which can be useful for environmental mapping and modelling studies. In addition, a good agreement was also found when TAGEE was compared to equivalent outputs from SAGA GIS, reaching a Pearson’s correlation coefficient between 0.96 and 0.98, and differences between 3–4 %. Thus, TAGEE becomes a feasible tool for making terrain analysis of big geospatial data, which can be customized to any spatial resolution and scaled up to the global extent.