ADAPTATION OF THE GLOBAL GEOID MODEL EGM2008 ON CAMPANIA REGION (ITALY) BASED ON GEODETIC NETWORK POINTS

The knowledge of the geoid undulation, the height of the geoid relative to a given ellipsoid of reference, is fundamental to transform the ellipsoidal heights into orthometric heights. Global geoid undulation models developed from satellite gravity measurements appropriately integrated with other data, are free accessible in internet, but their accuracy may be inadequate for specific applications. Earth Gravitational Model 2008 (EGM2008) is one of those: usually available in grid form 2.5’ x 2.5’ (a geotif is developed by Agisoft with resolution 1’x1’), it defines the difference between the WGS84 ellipsoid height and the mean sea level, but in some areas the discrepancies between these geoid undulations and local correspondent measured values are on the order of various decimetres. For consequence, more accurate models are necessary. This article aims to determine a geoid undulation model suitable for Campania Region (Italy), starting from the global model EGM2008 (1’ x 1’) that is locally adjusted by using geodetic network points (GNPs) and GIS interpolation functions. Three different datasets are considered including respectively 20, 40 and 60 GNPs and three deterministic interpolators are applied in global way to generate geoid undulation grids: Inverse Distance Weight (IDW), Global Polynomial 1 order (GP1), Global Polynomial 2 order (GP2). The resultant 9 models are tested on 20 additional GNPs. The experiments demonstrate that local geoid can be produced on a little area adapting global geoid by means of GNPs: the model obtained using GP2 and 60 GNPs, the most accurate one, fits the data with ±3.2 cm root mean square error (RMSE).


INTRODUCTION
A geoid undulation model (GUM) describes the vertical separation between the WGS84 ellipsoid and a hypothetical surface corresponding to mean sea level and its imagined extension under (or over) land areas (Pugh, 1987). The orthometric height, the geometrical distance between the geoid and the point on the topographical surface, is crucial for many engineering and geoscience applications, e.g., to build Digital Terrain Models (DTMs) (Maglione et al., 2014). Spirit levelling is a very accurate method for acquiring orthometric heights (Heiskanen and Moritz, 1967). It allows to determine the elevation of points or differences in elevation to produce necessary data for mapping, engineering design, and construction (Herbert and Olatunji, 2021). The axis of the spirit level is perpendicular to the plumb line and the instrument is locked in position, while two calibrated rods are held in an upright position ahead of and behind the instrument: the difference between readings on two calibrated rods is the difference in elevation between the points (Schomaker and Berry, 1981;Nestorović and Delčev, 2014). The orthometric heights can be determined also using trigonometric levelling which involves measuring a vertical angle from a known distance with a theodolite and computing the elevation of the point. In this case, vertical angle can be measured at the same time horizontal angles are measured for triangulation: this method results more economical but less accurate than spirit levelling (National Geodetic Survey, 2007). The orthometric heights can be obtained in correct way by GPS (Global Positioning System) survey only when an accurate GUM is available. In fact GPS observations together with geoid estimates can give orthometric heights in a faster and cheaper way than using other techniques (spirit levelling survey), although with lower precision which is however sufficient in many practical applications (Barzaghi et al., 2002). Earth Gravitational Model 2008 (EGM2008) (Pavlis et al., 2008;Pavlis et al., 2012) is a global geoid height model derived from satellite gravity measurements suitably integrated with other data, but its accuracy is poor, so local accurate models are necessary. Different approaches can be used for generating a local geoid (Marti, 2007, Falchi et al., 2018. GPS/levelling surveys allow to acquire both ellipsoidal and orthometric heights and the differences between those values in a good number of points can be interpolated for determining a GUM. In alternative way, geodetic network points (GNPs), the location of each being fixed in relation to both geoid and WGS84 ellipsoid, can be used for the purpose. Some main problems afflict the determination of the local geoid starting from GNPs: the accuracy of the geoid undulation values, the characteristic of the geoid surface in the area, the availability of an adequate number of data, and the choice of the most appropriate interpolation algorithm. Ad hoc survey based on spirit levelling and GPS techniques for acquiring very accurate orthometric and ellipsoidal height values is to prefer, but the existing data concerning GNPs can be a good starting point. There are different interpolation algorithms and each of them can have different results when interpreting your data (Al-Krargy et al., 2014;Erol and Erol, 2021).There is no absolutely best interpolation method but only the optimal choice under certain circumstance (Yang et al., 2004). The availability of a global geoid model and the adaptation of it to a local area based on GNPs and interpolation process can be a good integrated approach to increase the result accuracy. This article aims to determine a GUM for the territory of the Campania Region (Italy), starting from the global model EGM2008 (1' x 1') and adapting it to the specific local reality by using GNPs and GIS interpolation functions. Three different deterministic interpolators are applied in global way: Inverse Distance Weight (IDW), Global Polynomial 1 st order (GP1), Global Polynomial 2 nd order (GP2). The paper is organised as follows. Section 2 describes the materials and methods: 3 different datasets are selected including respectively 20, 40 and 60 GNPs which are processed as Ground Control Points (GCPs) to best adapt EGM2008 global model to the Campania Region. Section 3 presents and discusses the results comparing the levels of accuracy of 9 GUMs, 3 for each interpolation algorithm in dependence of the number of the GNPs including in each dataset; particularly, the accuracy is tested using a a fourth dataset comprising 20 GNPs, different from those previously used as GCPs, now introduced as Check points (CPs). Section 4 draws out our conclusions.

Datasets
In this study we use 80 GNPs located in Campania Region or in the neighboring and surrounding areas; these points are chosen in a uniformly distributed way and organized as follows: 60 points as GCPs for data processing and local geoid production, the other 20 points as CPs for testing the results achieved. Most of them are trigonometric vertices (66) or GNSS (Global Navigation Satellite System) permanent stations (3) located in Campania, 1 is the sea level station of Gaeta (Latina) within the Italian Network managed by ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale), 10 are trigonometric vertices in other Italian Regions (1 in Lazio, 4 in Molise, 4 in Puglia and 1 in Calabria). The choice of points outside the Campania region is necessary to ensure that the interpolation process is able to generate local geoid model covering the entire regional area which extends within the following UTM/WGS84 plane coordinates -33T zone: E1 = 396,255 m, E2 = 568,569 m, N1 = 4,426,808 m, N2 = 4,595,462m. The 60 points to be used for this purpose are divided into 3 groups of 20 elements each, namely A, B and C. A constitutes subset 1 containing 20 GNPs, A + B is subset 2 containing 40 GNPs, A +B + C is subset 3 containing 60 GNPs. The remaining 20 GNPs form the check set. All GNPs are represented as vector points in ArcGIS 10.3 using UTM WGS84 Zone 33 as reference system. Figure 1 shows the 60 GNPs used for adapting EGM2008 on the Campania Region and generating local geoid models; figure 2 shows the 20 GNPs used as CPs for testing these models.  For all GNPs the orthometric height referred to the national geoid named ITALGEO2005 (Barzaghi et al., 2007) and ellipsoidal height referred to WGS84 are available as reported in the monograph, the official document that consists of a sheet of useful information concerning every point, i.e. point location, univocal number code, coordinates, etc.

EGM2008
EGM2008, the successor of EGM96 and EGM84, is supplied by the EGM Development Team of the U.S. National Geospatial-Intelligence Agency (NGA). This gravitational model is provided as a set of normalized, geopotential coefficients complete to spherical harmonic degree and order 2159, and contains additional harmonic coefficients extending to degree 2190 and order 2159. The harmonic_synth FORTRAN 77 software provided by EGM development team allows to compute the gravity field quantities of EGM 2008; particularly, hsynth WGS84.f is the version of harmonic synthesis program that determines the geoid undulation with respect to WGS84 ellipsoid at any WGS84 latitude/longitude coordinate pair listed in a coordinate input file (Herath Mudiyanselage, 2014). In addition, it is also provided as a 2.5-minute worldwide geoid height file, precomputed from the EGM2008. Over areas covered with high quality gravity data, the differences between EGM2008 geoid undulations and independent GPS/Leveling values are on the order of ±5 to ±10 cm (Pavlis et al., 2012). However, higher discrepancies are registered in some areas , e.g. in Saudi Arabia (RMS of 79 cm) (Alothman et al., 2014). For this study we use the EGM2008 model as geotif file presenting a cell size of 1.0 x 1.0 minute developed by Agisoft (Agisoft, 2021) and free available in web for download, resulting in a global grid of 10,801 rows x 21,600 columns covering the entire earth surface; this grid contains 4 byte IEEE float values defining the difference between the geoid surface and the WGS84 ellipsoid. In figure 3 the EGM2008 geoid height model in the study area (Campania region) is shown.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W4-2021 16th 3D GeoInfo Conference 2021, 11-14 October 2021, New York City, USA Figure 3. The EGM2008 geoid height model in the study area projected in UTM-WGS84 Zone 33N.

Calculation of Geoid undulation.
Initially, the geoid undulation N of each GNP is calculated as the difference between the ellipsoidal height (h) and orthometric height (H) reported in the corresponding monograph, in accordance with the formula (Heiskanen and Moritz, 1967;Featherstone et al, 1998;Pepe and Prezioso, 2015;Maglione et al., 2018): The relationship between H, h and N is shown in figure 4. For each GCP, two values of geoid undulation are obtained, the one derived from the EGM2008 model and the other derived from the monograph of the point. The values of the deviation Δ between the two undulations are thus obtained, to be used subsequently for the interpolations.

Spatial interpolation
Spatial interpolation is the process of using points with known values to estimate values at other points. For consequence, spatial interpolation allows to create surface data from a limited number of sample points. It is largely used in GIS application to predict values for cells in a raster from any geographic point data, e.g. elevation, sea depth, temperature, concentrations of pollutants in the atmosphere, rainfall, etc. Many interpolation methods are present in literature and a lot of them are available in Geostatistical Analyst (Johnsto et al., 2001), an extension of ArcGIS that provides tools for optimal predictions possible by examining the relationships between all the sample points and producing a continuous surface. In this study we consider three methods of those included in Geostatistical Analyst: Inverse Distance Weight (IDW), Global Polynomial 1 st order (GP1), Global Polynomial 2 nd order (GP2). These interpolators are applied to the values of the difference of undulations previously calculated for GCPs. We report below the main characteristics of the selected interpolators.

Inverse Distance Weight (IDW)
Based on the hypothesis that closer values are more related than further values, IDW uses measured values surrounding an unmeasured location to predict its value (Lam, 1983).
The assigned values to unknown points are calculated with a weighted average of the values available at the known points. In fact, values at sampled points are weighted by an inverse function of their distance from the point of interest (Li and Heap, 2008). This is the implementation of the intuitive concept that the closest known values must carry more weight in determining the interpolated value for any unmeasured location (Aguilar et al., 2005). Therefore, as the distance increases, the weights decrease rapidly (ArcGIS Development Team, 2021a) so points that are closer to known values are more influenced than points which are distant from them. The interpolation function is represented by the following formula: (2) Where: z = value calculated in the unmeasured grid point j; = value measured at the known point i; = distance between the unmeasured grid point j and the known point i; p = power assigned to the distance d .
Normally, IDW works as an exact interpolator: since the weights assigned to the data points are fractions and the sum of all the weights is equal to 1.0, when a specific observation is coincident with a grid node, the distance between that observation and the grid node is 0.0, so a weight of 1.0 is given to that observation and weights of 0.0 to all other observations; consequently, the value calculated in the grid node is equal to the observed one (Yang et al., 2004). IDW can be applied both as a global interpolator and as a local interpolator. In .the first case, it calculates predictions using the entire dataset. In the second case, it calculates predictions from the measured points within neighbourhoods, which are smaller spatial areas within the larger study area The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W4-2021 16th 3D GeoInfo Conference 2021, 11-14 October 2021, New York City, USA

Global Polynomial Interpolators
Global polynomial interpolation (GPI) is an approximate method that fits a smooth surface defined by a mathematical function (a polynomial) to the input sample points (ArcGIS Development Team, 2021b). The user can choice the order of the polynomial that ranges from a first-order to higher order. The interpolation function can be written as: If n is the order of the equation, the following relations are valid: The values of the coefficients a i,j are determined using the known values in the sample points. Low-order polynomials generate gradually varying surfaces testifying their capability to describe some physical process; however, the more complex the polynomial, the more difficult it is to attribute physical meaning to it (Apadyn et al., 2004). GPI by estimating parameter , uses ordinary least squares to minimize the squared differences between the surface and measured points. The estimation of the coefficients permits to use the equation (3) for calculating the value of the polynomial function at any point within the map area (Wang et al., 2008). In this study, the orders 1 and 2 are considered. A first-order global polynomial (GP1) fits a single plane through the data according to the following formula: = 00 + 10 + 01 (7) A second-order global polynomial (GP2) fits a surface with a bend through the data according to the following formula: = 20 2 + 02 2 + 10 + 01 + 11 + 00 (8) Generally, GPI methods supplies surfaces which are highly susceptible to outliers (extremely high and low values), especially at the edges (ArcGIS Development Team, 2021b). However, this situation is typical of landscape which are rapidly variable alternating troughs and elevations. Nevertheless, for the undulations of the geoid in a limited area, the GPI methods may be effective.

RESULTS AND DISCUSSION
The accuracy of the EGM2008 is analysed considering the dataset of the 20 CPs. Particularly, the difference between the geoid undulation obtained from the global model and that obtained from the monograph is calculated in each CP.  Table 4. Statistical values of the residuals obtained in 20 CPs in the case of GUMs based on 60 GNPs.
By analyzing table 2 you can see how, using 20 GCPs, the best results are provided by GP1 as the RMSE value measured on the CPs is equal to 4.6 cm (5.3 cm by IDW and 4.9 cm by GP2). The table 3 shows how using 40 GCPs the best interpolator is GP2 which provides an RMSE equal to 3.7 cm (4.1 cm by IDW and 3.9 cm by GP1). The table 4 confirms the good performance of GP2 which provides, in the case of 60 GCPs, a RMSE equal to 3.2 cm (3.7 cm by IDW and 3.5 cm by GP1).
From the values obtained it is also possible to comment on the improvements or worsening of errors according to the number of points used. For example, the IDW interpolator, while presenting the worst performance among the three interpolators considered, records a marked improvement in results as the number of points used increases. In fact, by increasing this number from 40 to 60, the RMSE decreases by 1.6 cm. If we take into account the difference between the EGM2008 and the undulation of the Check Points derived by the monographs (table 1), we note that even with the IDW a clear improvement is introduced. In fact, the RMSE value with 60 points is 3.7 cm compared to 29.9 cm found for the EGM2008.
In the results, it is seen that the global deterministic interpolation methods are applicable for adapting EGM2008 on the study area and generating local geoid. The increasing number of GNPs homogeneously distributed in the study area allows to achieve more performing models. However, GP2 fits the data better than IDW as well as GP1 methods according to test criteria. In particular, the model obtained by the application of GP2 interpolator and 60 GNPs is the most accurate GUM fitting the data with ±3.2 cm (RMSE). Figure 5 shows the most accurate local geoid model resulting from the adaptation of EGM2008 on the Campania Region based on 60 GNPs and the GP2. Figure 6 shows the same model in 3D visualization: since they are much smaller than the size of the area under examination and are extending in a limited range, geoid undulation values are multiplied by an amplifying factor equal to 10.000, so the enhance their variability in threedimensional view.

CONCLUSIONS AND FUTURE WORKS
The results demonstrate that a local geoid can be produced on a little area adapting global geoid by means of GNPs. Usually these points present both orthometric and ellipsoidal heights, so geoid undulations can be easily calculated. The experiments carried out on Campania Region also prove that a limited number of GNPs can contribute to define in rapid way a local geoid when they are appropriately chosen and used in combination with a global model, i.e. EGM2008.
Particularly, these points must be uniformly distributed over the area affected by the model and also cover the neighbouring areas in order to avoid extrapolation. The adopted approach is based on two steps: firstly, the differences between the GCPs geoidal heights and the corresponding EGM2008 1' x 1' ones are interpolated using a fixed interpolation method (i.e. IDW, GP1 or GP2); then the resulting model is summed to the EGM2008 1' x 1' to obtain the GUM. All operations are carried out using UTM-WGS84 Zone 33 N (cell size: 150 m) The results obtained show how, by varying both the type of spatial interpolator and the number of GCPs used, it is possible to obtain different models characterized by different accuracy values. In fact, fixed an interpolator, the greater the number of points, the greater the accuracy of the model. Taking for example GP1, if 20 points are used to interpolate the geoid undulation values, the final RMSE is equal to 4.6 cm, while with 40 points it will improve up to 4.0 cm. With GP2, the passage from 40 to 60 GCP determines a decrease in the RMSE value which goes from 3.8 cm to 3.2 cm. GP1 always provides better results than IDW. GP2 is always more performing than the other two interpolators, with the exception of the case in which only 20 GCPs are considered, but in this situation the difference in terms of RMSE with GP1 is only 3 mm. The most accurate model is the one obtained with 60 GCPs and with GP2.
Concerning the future developments of this work, further experiments will be carried out integrating datasets with other GNPs in order to assess the impact of the increasing number of points on model accuracy. In addition, other interpolation methods will be considered, e.g. local deterministic and geostatistical ones.
ArcGIS Development Team, 2021b. How global polynomial interpolation works. In ArcGIS Pro help, ESRI, Redlands, CA, USA.