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

: 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 3x3 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.


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 package code and a minimal reproducible example are available in https://github.com/zecojls/tagee (Supplementary Materials).
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 globalscale 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.

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.

Topographic Surface
The land topography can be approximated by a topographic surface defined by a continuous, single-valued bivariate function (Equation (1)) [8]: where is elevation (meters), and and 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 = ( , ) is expressed as the second-order bivariate Taylor polynomial (Equation (2)): where , , , , and are the partial derivatives, and 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 3x3 moving window must retrieve both the geometry elements and the elevation values of the window nodes ( Figure 1).

Terrain Parameters: Neighbor Elevations and Geometries
The elevation values of a 3x3 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), (4), (5)): where is latitude for the first given node in radians, is the latitude for the second given node in radians, is the longitude for the first given node in radians, is the longitude for the second given node in radians, Δ and Δ are the respective differences of latitude and longitude between the given nodes, and is the mean radius of Earth equals to 6,371,000 meters. The linear distance 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 to ) to calculate the partial derivatives of terrain.

Terrain Attributes
Local attributes, such as slope, aspect, and curvatures, are calculated from the partial derivatives of terrain [8]. The slope gradient ( , Equation (11)) is a flow attribute that relates to the velocity of gravity-driven flows. For measuring the direction, the slope aspect is used ( , Equations (12) and (13)). Additionally, one can calculate the amount that a slope is faced to the North or East, resulting in the Northernness ( , Equation (14)) and Easternness ( , 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 ( , Equation (16)) and vertical curvatures ( , Equation (17)). While the horizontal curvature relates if a lateral flow converges ( < 0) or diverges ( > 0), the vertical curvature measures the relative acceleration ( > 0) and deceleration ( < 0) of a gravity-driven flow: Differently from flow attributes, which are gravity field-specific variables, form attributes are related to principal sections of terrain [8]. The mean curvature ( , 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 ( , Equation (19)) is a product of maximal ( ) and minimal ( ) curvatures. The two principal curvatures calculate the highest and lowest curvature for a given point of the topographic surface. The maximal curvature ( , Equation (20)) is useful for mapping rigdes ( > 0) and closed depressions ( < 0). Likewise, the minimal curvature ( , Equation (21)) is useful for identifying hills ( > 0) and valleys ( < 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 ( , Equation (22)) ranges from -1 to 1 and map convex ( > 0) and concave ( < 0) landforms:

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 3x3 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 to ) and the distances ( , , , , and ) ( 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 ( , , , , and ) 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). 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.

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 3x3 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 2 nd 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].

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). 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 ( Figures 3A, B and 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-touse 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 highorder 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.

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.
Supplementary Materials: The package code and a minimal reproducible example are available in https://github.com/zecojls/tagee.