Transformation from a global spherical to an adjusted local rectangular harmonic model

This work presents a technique to transform a global spherical to an adjusted local rectangular harmonic model. First, the mathematical form of a global spherical harmonic model is presented. Second, the necessary conversion from global (geocentric) into local rectangular coordinates is given. Third, Laplace’s equation is solved by the method of separation of variables in local rectangular coordinates and its solutions in different functional forms are presented. Then, the estimation of the coefficients of these mathematical models by a least squares’ adjustment process is described, using as data the values of the disturbing potential of the Earth’s gravity field. The strategy for the selection of the best mathematical model for a successful transformation is described and validated in different case studies. These refer to areas in Greece, China and Germany and include comparisons with other models or methods. The results show the applicability of the presented transformation and confirm its advantages.


Introduction
For the representation of the Earth's global (geomagnetic or gravity) field, spherical harmonics are widely used, which are solutions of Laplace's equation found by the method of separation of variables. Nowadays, the functionals of the geopotential can be computed from spherical harmonic models which are mainly derived from satellite measurements and become more and more detailed and accurate (Barthelmes 2013). On the other hand, for the representation of a local or regional field, Spherical Cap Harmonic Analysis (SCHA) (e.g. Haines 1985) and Rectangular Harmonic Analysis (RHA) (e.g. Alldredge 1981) are the preferred alternative techniques over the last decades. In comparison to classical spherical harmonic models, the resulting harmonic models from the above two techniques have two important advantages: First, the local or regional field represented by an analytical model needs much less parameters (coefficients) and has less memory requirements in computation and storage than a spherical harmonic model. Second, these models can be easily enhanced using new data of different kinds (e.g. gravity anomalies, geoid undulations), leading to a further refinement of the local field. For these reasons, the two techniques are now being used more and more in numerous problems of the geoscientific community.
A comprehensive and systematic review of the SCHA literature is presented by Torta (2020), underlining the respective merits and weaknesses of the ways in which the technique has been used since it was first proposed in the context of a regional geomagnetic field modeling by Haines (1985). Several works followed, where the main problem of this approach was the computation of the associated Legendre functions with non-integer degree (e.g. De Santis 1992;Younis 2013). To overcome this problem, a new approach, called Adjusted Spherical Cap Harmonics (ASCH), was proposed by De Santis (1992) using the well-known integer order and degree Legendre functions. In order to avoid the same problem, another new approach, called Virtual Spherical Harmonics (VSH), was proposed by Wang and Wu (2019) in regional geoid modeling. An application of modeling regional gravity fields in an integrated approach, using ASCH by a least squares' estimation of the coefficients of the mathematical model, was presented by Younis (2013). Also, the transformation from a global spherical to an adjusted local spherical cap harmonic model was developed by Younis et al. (2013).
RHA was first introduced by Alldredge (1981) and it has been used in geomagnetic and gravity field modeling (e.g. Alldredge 1981Alldredge , 1982Alldredge , 1983Haines 1989Haines , 1990aHaines , 1990bMalin et al. 1996;Jiang et al. 2014). In the works of Alldredge (1981Alldredge ( , 1982Alldredge ( , 1983,  and , RH is applied only in geomagnetic fields. In this approach, Laplace's equation is solved by the method of separation of variables in local rectangular coordinates. As mentioned in Jiang et al. (2014), the computation of the trigonometric basis functions of the RH is more stable and efficient at high degree expansion than that of the associated Legendre functions of the SCH. The mathematical properties of rectangular harmonic expansions are studied in the works of Haines (1989Haines ( , 1990aHaines ( , 1990b. Malin et al. (1996) re-examined Alldredge's method of rectangular harmonic analysis and modified the functional form of the mathematical model. In the work of Jiang et al. (2014), a mathematical model of the disturbing potential of the gravity field and its functionals (gravity anomaly and disturbance, geoid undulation, deflection of the vertical) was derived and a strategy was proposed in order to reduce the edge effects caused by periodic continuation in RH. Additional recent applications of RH include an iterative procedure of the weighted least squares' prediction, which was used by Zhang et al. (2017) for merging the land, altimetric, shipborne and airborne gravity data. Also, Feizi et al. (2021) investigated the problem of local gravity field modeling over Tanzania mountainous regions using RH of vector airborne gravity data.
In this study, a global spherical geopotential model is transformed to a rectangular harmonic expansion, over a relatively small area, by a least squares' adjustment process, using as input data values of the disturbing potential of the Earth's gravity field. In contrast to previous works on RH, here we use the most general solution of Laplace's equation, including the term corresponding to the case in which the constants of the separation of variables are zero. Then, we apply different conditions and derive additional mathematical models which can be validated through a numerical process. The strategy for the selection of the best mathematical model for such an application and, specifically, the selection of expansion degree and extension parameter of a model, is described in different case studies.
Although this work resembles the work of Jiang et al. (2014), there are two important differences between the two works: First, here we provide the different functional forms of the solution of Laplace's equation and we develop a strategy for the selection of the best mathematical model for a successful transformation. Second, the final models of the two works (Eq. (8) in the work of Jiang et al. (2014) and Eq. (17) in this work) are formally different and for this reason the models are validated through a numerical process for the same data (Sect. 6.2). Comparing the results, the model used here achieves better performance, a result that is also justified by the different number of parameters between the two models. Specifically, the model of Eq. (8) in the work of Jiang et al. (2014) consists of 1 + 2N + 2M + 4NM terms, whereas Eq. (17) of this work contains 8 + 4N + 4M + 4NM . The 7 + 2N + 2M non-periodic extra terms are included for two reasons: (i) The general functional form is essentially an interpolating function and the non-trigonometric terms (such as linear trends) or the exponential term are retained in order to keep down the problems associated with the rectangular approximation of the spherical geometry of the problem, as pointed out in Haines (1990a) (ii) These terms limit the problems of an ill-posed system of normal equations, as proved through the presented numerical tests and other tests that may be worked out using the uploaded code.
This work is organized as follows: Sect. 2 describes the global spherical harmonic models used to compute the required data (values of the disturbing potential). Section 3 presents the necessary transformation from global (geocentric) into local rectangular coordinates, in which Laplace's equation is expressed. The solutions of this partial differential equation are determined by the usual method of separation of variables, leading to a combination of functions and thus in different mathematical models. Section 4 describes the adjustment process for the estimation of the coefficients of these functional models. Section 5 evaluates the proposed adjusted local rectangular harmonic models. Section 6 presents various case studies, demonstrating and validating the performance of this methodology. Finally, Sect. 7 presents our conclusions.

Global spherical harmonic models
The Earth's gravitational potential V at any point, with spherical coordinates (r, , ) on or above the Earth's surface, is conveniently expressed on a global scale by summing up over degree n and order m of a spherical harmonic expansion as follows: where GM is the product of gravitational constant and mass of the Earth, R is the reference radius, P nm are the Legendre functions and C nm , S nm are Stokes' coefficients. Hence, a global spherical harmonic model of the gravitational potential consists of n max + 1 2 coefficients and the two values for GM and R to which the coefficients relate (Barthelmes 2013).
On the other hand, the normal gravitational potential V at any point with oblate spheroidal coordinates (u, , ) outside the level spheroid, is given by the following formula (Heiskanen and Moritz 1967): is the linear eccentricity of the spheroid with semi-axes a and b and is its angular velocity. Also, the transformations between the global (geocentric) rectangular coordinates (X, Y, Z) and spherical (r, , ) , spheroidal (u, , ) and geodetic , , h coordinates can be found e.g. in Heiskanen and Moritz (1967).
Finally, the disturbing potential T at any point with global rectangular coordinates (X, Y, Z) outside the masses is given by which does not contain the centrifugal potential and is a harmonic function.

Laplace's equation in local rectangular coordinates and its solutions
The disturbing potential may also be expanded into a series of rectangular harmonics. To find them, we introduce local rectangular coordinates x , y and z . The z-axis is taken upward and the y-axis and x-axis northward and eastward, respectively, on the plane normal to a line which is perpendicular to an oblate spheroid (e.g. WGS84) at the center 0 , 0 of the area under investigation. The global (geocentric) rectangular coordinates (X, Y, Z) of a point are transformed into these local rectangular coordinates (x, y, z) by the following relations: where X 0 , Y 0 , Z 0 are the global coordinates and 0 , 0 , h 0 are the corresponding geodetic coordinates of the origin of the local coordinate system.
It can be shown that Laplace's equation for the disturbing potential T , in local rectangular coordinates, is given by We shall attempt to solve this equation by separating the variables x , y , z by means of the trial substitution where f is a function only of x , g is a function only of y and h is a function only of z . Performing this substitution in Eq. (7) and dividing by fgh , we obtain which leads to the following three equations and where p and q could be any numbers (including zero).
Solutions of Eqs. (10) and (11) may be given by the functions On the other hand, solution of Eq. (12) may be given by which, for p or q positive, tends to zero when z tends to infinity. Since the potential cannot diverge at large values of z , the positive exponent will be absent. Moreover, negative values of p and q are meaningless due to the periodicity in Eqs. (13) and (14) and the presence of squares in Eq. (15), therefore we only distinguish cases of zero and positive values. Also, to make the above solutions meet the needs of the problem under consideration, we replace the values p and q with (Alldredge 1981) where n (degree) and m (order) are integers and D x and D y are the overall lengths in the x and y directions, respectively, that cover the entire area where the data are located. Thus, the most general solution for the disturbing potential, as linear combination of the solutions (13), (14) and (15), in terms of truncated series at maximum degree N and order M , is expressed in the form This general solution (say model M0) has n RHC = 8 + 4N + 4M + 4NM Rectangular Harmonic Coefficients (RHC). However, we can restrict the terms of the infinite series in Eq. (17) and obtain five additional models (M1 to M5), which are presented in Table 1, along with the number of the corresponding RHC.
The spatial resolution represented by the above models, which is characterized in terms of half wavelength, can be computed from N and M as (Jiang et al. 2014): Finally, the RHC in each of the previous models can be estimated by the method of least squares, as described in the next section.

Adjustment process of the transformation
Using the method of least squares (see e.g. Ghilani and Wolf (2006)) we can estimate the RHC of the models of Table 1 and, therefore, to transform a global spherical to an adjusted local rectangular harmonic model, which is practically an interpolating function. To do this, a 3D grid of k points is generated over the area under investigation, where the number of points ( k ) is larger than the number of RHC ( n RHC ). Then, the values of the disturbing potential, at the grid points, are computed by Eq. (5) using a global spherical harmonic model and the function of the normal gravitational potential, for the same oblate spheroid of reference. These values are used as observations for an arbitrary point i , where i = 1, … , k . Then, we set up the linear observation equation e.g. for the model M0: CC nm cos p n x cos q m y + CS nm cos p n x sin q m y and y = D y

2M
where n = 1, 2, ... , N and m = 1, 2, ... , M. Similar equations can be written for the other models, although here we describe the process for the general model M0. Hence, the system of k linear equations is represented in matrix form as where the k × n RHC (design) matrix contains the derivatives with respect to the elements of the n RHC × 1 vector of unknowns The k × 1 vector of observations contains the values of the disturbing potential and is a k × 1 vector of residuals. The solution of the adjustment is given by where the weight matrix , which is diagonal if the observational errors are uncorrelated, is computed from the variances of the observations. Furthermore, the variance factor ̂ 2 0 is computed by while the vector of residuals from Eq. (24). Also, from the residuals, we can compute the statistical index Root Mean Square (RMS) Finally, the variance-covariance matrix of the estimated parameters is given by In this work, the values of disturbing potential (observations) are assumed unweighted and uncorrelated and thus the weight matrix is equal to the identity matrix, i.e. = . Also, following the suggestion of Jiang et al. (2014), in order to avoid the edge effects, we extend the sizes D x and D y of the data domain (only the sizes, not the data domain) to D x + d x and D y + d y , respectively, although numerical instabilities may occur in the adjustment process. Then, Eq. (16) is modified as Clearly, in order to minimize the RMS, which is influenced for a specific model (e.g. M0) by the maximum degree N and order M (expansion degree and order) and the sizes of d x and d y (extension sizes), we have to examine this dependence for various values of these parameters. However, for simplicity reasons, we assume N = M , d x = d and d y = d , where is a known real positive number, so that we only have to examine the values of the two parameters N and d , which minimize the RMS of the disturbing potential, as we show in the following case studies. As an upper bound of the extension parameter d , for and q m = m 2 D y + d y the restriction of numerical problems, should be used the size of the data domain. Furthermore, for the inversion of the symmetric matrix in Eq. (26), we use the method of Cholesky, which has been widely accepted in adjustments. Initially, we have to select the best mathematical model among the different functional forms presented in Table 1.

Validation of different adjusted local rectangular harmonic models
In order to select the best mathematical model (M0 to M5) which minimizes the RMS of the disturbing potential, we conduct numerical experiments based on data generated from the geopotential model EGM2008 with maximum degree 2190 (Pavlis et al. 2012). A 3D grid of k = 6724 points in Greece, with spatial resolution 3 � × 3 � , was generated over a 2 • × 2 • area from ( , ) = (37 • , 22 • ) to ( , ) = (39 • , 24 • ) . The values of the disturbing potential were computed at four ellipsoidal height levels (0 m, 100 m, 1000 m and 2000 m) via the ICGEM online platform (Ince et al. 2019). The WGS84 reference ellipsoid was used to compute the normal gravity field parameters while maintaining the model's defined tidal system. Gaussian filtering was not applied and the zero-degree term was taken into consideration. Then, the global (geocentric) rectangular coordinates (X, Y, Z) of these points were transformed into the corresponding local rectangular coordinates (x, y, z) through Eqs. (6), using as origin of the local coordinate system the point 0 , 0 , h 0 = (38 • , 23 • , 0) . Also, the parameters of Eq. (30) were D x = 180 km and d x = 0.80d , while D y = 220 km and d y = d.
Starting with model M0, we examined the performance for different values of expansion degree or order ( N = M ) and extension parameter d . The resulting RMS values and the execution time for each adjustment process are shown in Table 2. The minimum RMS value for the disturbing potential was produced with degree N = 20 and extension parameter d = 125 km.
We repeated the same process for the other models (M1 to M5) and different values of expansion degree and extension parameter. The results of minimum RMS values are shown in Table 3. The experiment reveals that the RMS value of model M0 is smaller compared to those of the other models. Concluding, the numerical results confirm that the functional form M0 is able to provide a reliable and precise adjusted local model for the representation of a local gravity field, as we demonstrate and in the following case studies.
The execution times reported in the tables correspond to a single adjustment (one combination of N and d values). The computations were performed on a personal computer running a 64-bit Linux Debian 5.9 operating system. The main characteristics of the  hardware were: i7-5820 K CPU (with 6 cores-12 threads, clocked at 4.0 GHz) and 16 GB of RAM. All algorithms were implemented in C ++ and compiled using the GNU Compiler Collection version 10.2 at optimization level 3. We also employed the open-source "libquad-math", the GCC Quad-Precision Math Library, which provides a precision of 33 decimal digits.
6 Case studies and numerical results

Local rectangular harmonic model in Greece
Having selected the model M0, we extended the previous area in Greece and a 3D grid of k = 26244 points, with spatial resolution 3 � × 3 � , was generated over a 4  Table 4, where the symbol "-" indicates that the corresponding adjustment was omitted. The minimum RMS value for the disturbing potential was produced with degree N = 35 and extension parameter d = 125 km and from these values we produced Fig. 1. The values of the residuals plotted against the distance from the system origin are shown in Fig. 2 (Left). In order to validate the resulting local model, a new grid of 6400 evaluation points was created in the same area. These were located intermediate to the adjustment points and had a single height value, corresponding to the actual height of the topography (for points in the sea, we used the corresponding height of the geoid). The differences of the values of the disturbing potential between the reference global and the local model are presented in Fig. 2 (Right).
The histogram of residuals of the disturbing potential, presented in Fig. 1, and the maximum values of residuals ( Fig. 2-Left) correspond to ±0.1 mm precision of geoidal heights, which indicates that the model is successfully transformed. As it is known, the geoid undulation N is related to the disturbing potential T through the Bruns formula N = T∕ , where is the normal gravity (≈ 10m∕s 2 ) . Furthermore, the differences of the disturbing potential between the values from the global and local model, at the grid of validation points (Fig. 2-Right), reflect the reliability and accuracy of the local model. Comparing the above results with those achieved for the smaller ( 2 • × 2 • ) area of Greece, with the same origin of the local coordinate system, we observe that the

Local rectangular harmonic model in China
For a comparison with another mathematical model of RH, in this case study we generated the grid described in the work of Jiang et al. (2014). The 5 • × 5 • data area in China is limited by the parallels of 29° and 34° northern latitude and by the meridians of 108° and 113° eastern longitude, with D x = 470 km, D y = 551 km and extension parameter d x = d y = d . The data area was sampled with 14,641 points, spaced by 2.5 � × 2.5 � in latitude and longitude and at the corresponding height of the geoid.
The RMS values for the disturbing potential from the adjustment of model M0 with different degrees and extension parameters are given in Table 5. From these values, we select the expansion degree N = M = 40 and extension parameter d = 75 km. The results are shown in Fig. 3.
We note that, as in the previous case study, we validated the model using 5184 points, located at intermediate positions of the adjustment points, in the area between parallels of 30° and 33° northern latitude and meridians of 109° and 112° eastern longitude.
It can be seen from Fig. 3 (Left) that the residuals of the disturbing potential (corresponding to ±10 μm precision of geoidal heights) and the differences of the disturbing potential between the values from the global and local model at the validation points ( Fig. 3-Right), both reflect the reliability and accuracy of the successfully transformed local model. Comparing the above results (Table 5) with those achieved by Jiang et al. (2014), presented in Table 1 of mentioned paper, we obtain better results for the disturbing potential because of the use of a different general mathematical model (M0). In the case of the 5°×5° area, the RMS values related to the values of extension parameter d (last column of Table 1 in Jiang et al. (2014)) are significantly greater than the RMS values computed for different expansion degree N and the same values of parameter d in our Table 5, converted in terms of geoid undulation using Bruns formula. In addition, the RMS values in Table 5, in terms of geoid undulation, are the same order of magnitude as the RMS values referring to the smaller area 3°×3° in Table 1 presented in Jiang et al. (2014). Finally, it is worth noting that the number of coefficients in the model used by Jiang et al. (2014) equals to (2N+1)(2M+1).

Local rectangular harmonic model in Germany
For a comparison with the method of ASCH by a least squares' estimation of coefficients of the mathematical model, in this case study we generated almost the same grid which was used in the work of Younis (2013). Thus, a 3D grid of k = 11163 points, with spatial resolution 3 � × 3 � , was generated over a 3 The RMS values for the disturbing potential from the adjustment of model M0 with different degrees and extension parameters are given in Table 6. From these values, we select the expansion degree N = M = 25 and extension parameter d = 125 km and create Fig. 4.
As in the previous cases, we validated the model using 3600 points, located at intermediate positions of the adjustment points and at the actual height of the topography.
It can be seen from Fig. 4 (Left) that the residuals of the disturbing potential correspond to ±0.1 mm precision of geoidal heights and the differences of the disturbing potential at the validation points ( Fig. 4-Right) are generally small and increase with the distance from the origin. However, from Table 6 we ascertain some problems of numerical   Jiang et al. (2014). In this work, we did not apply such methods because, as presented in Tables 2, 4, 5 and 6, the RMS value is smallest within the range of the parameters where the proposed transformation is still applicable and the numerical problems have not appeared.
Comparing the above results with those achieved by Younis (2013), from the same data, we obtain better results for the potential, mainly because of the use of rectangular harmonics instead of spherical cap harmonics. Furthermore, since the gravitational potential of the normal gravity field can be given analytically by Eq. (2), we find preferable to transform the disturbing potential rather than the gravitational potential. As already stated in the Introduction, spherical cap harmonics have their own problems and, additionally, they demand a gravitational constant of the Earth GM and a selected reference radius R.

Conclusions
A global spherical harmonic model requires a very high number of coefficients to represent the gravitational potential, even in a local area. In contrast, a local rectangular harmonic model requires less coefficients and hence memory requirements in computation and storage. For example, the global model EGM2008, with maximum degree 2190, includes 4,800,481 coefficients, while the corresponding local model produced in this work, for instance in a 4 • × 4 • area in Greece with maximum degree 35, includes only 5188 coefficients and the disturbing potential is accurately reproduced. Therefore, the transformation from a global spherical to an adjusted local rectangular harmonic model, as presented in this work, offers an advantage. This transformation is based on a mathematical model (M0), derived in this work, and uses a least squares' adjustment process for the estimation of the rectangular harmonic coefficients of the model.
The results in the various case studies demonstrate the reliability and accuracy of the successfully transformed local model. The derived model of the potential is simpler to Fig. 4 From the adjustment of model M0 with expansion degree N = M = 25 and extension parameter d = 125 km in 3 • × 3 • area in Germany: Residuals of the disturbing potential (Left) and differences of disturbing potential between the values from the global and local model at the validation points (Right), as a function of the distance from the origin implement than other analyses (e.g. SCH, VSH) and, in addition, one can also calculate its functionals, such as gravity anomaly and disturbance, height anomaly, geoid undulation and deflection of the vertical. On the other hand, a limitation of the rectangular harmonic model is that it cannot be extrapolated, as evidenced by the increased errors in the edge of the region constrained by the data. However, the main drawback of the method for a larger area, known as edge effects (Jiang et al. 2014), may be alleviated using advanced least squares' adjustment techniques.
The examination of the behavior of the residuals and their distributions, such as the statistical significance of the estimated parameters, should be an object of further study that could include measurements of the same or of a different kind, aiming at the refinement of the local model. On the other hand, the regional improvement of a global geopotential model using modern data such as GNSS/Leveling data leads to instabilities, such as the leakage effect (for details see the work of (Mosayebzadeh et al. (2019)). In any case, the first step is the transformation which has been presented in this study.