PGM2016: A new Geoid Model for the Philippines
R. Gatchalian, Geodesy Division, National Mapping and Resource Information Authority, Philippines
R. Forsberg, National Space Institute, Technical University of Denmark, Denmark
A. Olesen, National Space Institute, Technical University of Denmark, Denmark

Abstract. In 2014, a preliminary geoid model has been computed for the Philippines
i.e., Philippine Geoid Model 2014 (PGM2014), with the technical assistance of the
National Space Institute, Technical University of Denmark (DTU-Space) using datafrom land gravity, airborne gravity, marine satellite altimetry and the newest satellite
gravity data from the GOCE mission release 5. Digital terrain models used in thecomputation process was based on 15” SRTM data. The model is computed in a global
vertical reference system then fitted to the ITRF GNSS/Leveling and validated with
RMS value of 0.50m. In 2016, the PGM2014 was re- computed into PGM2016 using
the reprocessed and densified land gravity data (from 1261 to 2214 points). Significant
improvements can be seen in the reprocessed gravity data and GNSS/Leveling (RMS
= 0.040m). Further densification of the land gravity (in towns and cities) to 41,000
points will be conducted from 2017 until 2020 to refine the geoid. Re-computation will
be done for the new version of the geoid as new gravity data comes in. GRAVSOFT
system of FORTRAN routines developed by DTU-Space and Niels Bohr Institute,
University of Copenhagen was used in computing the Philippine geoid.

Introduction
Vertical coordinates (i.e. Heights) of points are referred to a coordinate surface called
Vertical Datum. The universal choice of a vertical datum is the geoid – the reference
surface for orthometric and dynamic heights (Vanicek, 1991). It is an equipotential
level surface of the oceans at equilibrium, proposed by C.F. Gauss as the“Mathematical figure of the earth” (Dr. Bernhard Hofmann-Wellenhof, 2005).
With the advent of Global Navigation Satellite Systems (GNSS), it has become mucheasier to estimate Mean Sea Level (MSL) elevations using a geoid model. Applying a
geoid model in GNSS surveys will eliminate the conduct of levelling. A geoid model isa surface (N) which describes the theoretical height of the ocean and the zero-levelsurface on land. In a modern vertical reference system, the geoid is required to obtainorthometric height H (“height above sea level”) from GPS by
𝐻 = ℎ𝐺𝑃𝑆 − 𝑁(1)

Where ℎ𝐺𝑃𝑆 is the GPS ellipsoidal height, and H the levelled (orthometric) height.
In the Philippines, determination of elevation H of points and Benchmarks (BMs) was
normally conducted thru Geodetic Levelling, a tedious process that hinders thedensification of BMs during the early times.

1

In 2007, during the Philippine Reference System 1992 (PRS92) Campaign, 22,851
BMs were established along major roads nationwide with a maximum divergence of
4.0𝑚𝑚√𝐾 between two level runs (1st Order Accuracy). These first order levelnetworks were connected to their respective reference tidal BMs to provide local MSL
elevation. Figure 1 shows the Network of Level Lines with their corresponding Tide
Gauge Stations.

Figure 1. Philippine Vertical Control Network established by Levelling

2The first attempt of computing a preliminary gravimetric geoid for the Philippines isthrough the Natural Resources Management Development Project (NRMDP) in 1991.
Land gravity data and altimetrically-derived anomalies at sea and OSU89A to degreeand order 360 (reference global model) were used. Biases between the gravimetric N
and GPS/Levelling were found ranging from 2-6 m nationwide (Kearsley, 1991).
Figures 2 and 3 shows the computed detailed geoids of the Philippines.

Figure 2. Detailed geoids of Northern (left) and Palawan with South West Philippines (right)

3

Figure 3. Detailed geoid of Visayas and Mindanao

4The making of Philippine Geoid Model 2016 (PGM2016)
In October 28, 2014 with the technical assistance of National Space Institute -Denmark
Technical University (DTU-Space) and funding from National Geospatial Intelligence
Agency (NGA), a preliminary geoid model, Philippine Geoid Model 2014 (PGM2014)
has been computed for the country using the data from land gravity, airborne gravity,
marine satellite altimetry and the newest satellite gravity data from the GOCE mission
release 5 with an accuracy of 0.30meters. In this paper, the computation of the
PGM2014 will be discussed then its re-computation into PGM2016.

The airborne gravity survey
The success of the first long-range airborne gravity survey in Greenland 1991-92 by
the group of US Naval research Laboratory, in cooperation with NOAA, NIMA andDanish National Survey paved the way for the use of airborne gravity in filling theintermediate wavelength bands between satellite gravity, e.g. GRACE and GOCE (R.
Forsberg & Olesen, 2010) for accurate geoid modelling. There is now full operational
capability of collecting seamless airborne gravity data across land and marine areas
at 1-2 mGal r.m.s. accuracy and at around 4-6 km resolution due to the new gravity
acceleration sensors and improved GPS processing methods.
After their successful airborne gravity campaigns in the other regions of the world (e.g.,
Malaysia, Mongolia, Ethiopia, South Korea, Nepal), DTU-Space conducted theairborne gravity survey in the Philippines on March until May of 2014 using a Cessna
Caravan aircraft. This is part of the project to improve the global gravity field model
EGM2008 under the umbrella of the NGA - Danish Geodata Agency BECA agreement.
Figures 4 and 5 shows the aircraft used in the Philippine airborne gravity survey andits cabin layout respectively.
The following instruments were used:- LaCoste and Romberg Air/Sea gravimeter S-38
- Chekan AM gravimeter #24- Javad Lexon GPS receiver- Javad Delta GPS receivers- Novatel dual frequency aircraft GPS antenna- Power rack with DC/AC inverters, UPS etc.
GPS reference stations operated in all airports together with active geodetic stationsof NAMRIA were used as base stations in computing the position of the aircraft.
AUSPOS online positioning provides the coordinates of the reference stations (inITRF2008) with an estimated accuracy of 0.50 to 2cm in vertical. Aircraft trajectories
were computed with waypoint software package from Novatel (Calgary, Canada) using
precise ephemeris from International GNSS Service (http://igscb.jpl.nasa.gov/).
Parking spots were tied to NAMRIA first order gravity stations using Scintrex CG-5
gravimeter.

5Mean altitude for all flights was 3185m with a terrain clearance of 545m above
mountains and 3760m in lowlands. Figure 6 shows the color-coded flight track
elevations.

Figure 4. The Cessna Grand Caravan aircraft used in the airborne survey

Figure 5. Cabin layout in the Cessna Caravan aircraft. To the right, the LCR gravimeter andbehind it the Chekan meter6

Figure 6. Flight track elevations of the airborne gravity survey

Free –air gravity anomalies at aircraft level are obtained from:
𝜕𝛾

𝜕2 𝛾

∆𝑔 = 𝑓𝑧 − 𝑓𝑧0 − ℎ" + 𝛿𝑔𝑒𝑜𝑡𝑣𝑜𝑠 + 𝛿𝑔𝑡𝑖𝑙𝑡 + 𝑔0 − 𝛾0 − (𝜕ℎ (ℎ − 𝑁) + 𝜕ℎ2 (ℎ − 𝑁)2 )(2)Where 𝒇𝒛 is the gravimeter observation, 𝒇𝒛𝟎 the apron base reading, 𝒉" the GPS
vertical acceleration, 𝒈𝒆𝒐𝒕𝒗𝒐𝒔 the Eotvos correction computed by the formulas ofHarlan(Harlan, 1968), 𝒈𝟎 the apron gravity value, 𝛾0 normal gravity, h the GPS
ellipsoidal height and N the geoid undulation from EGM08 2.5x2.5minutes, 𝛾0 and thesecond order height correction is based on GRS80 definitions (Moritz, 1980). All
altitude dependent atmospheric correction has been applied, see e.g. Hintze, et al.
2005. The platform off-level correction 𝑔𝑡𝑖𝑙𝑡 is based on a platform modelling approachdescribed in (Olesen, 2002). All data were filtered with a symmetric second order
Butterworth filter with a half power point at 170 seconds, corresponding to a resolutionof 6km (half-wavelength).

7

Apron base readings were performed each day after the flight and approximately every
second day before the flight to monitor drift of the gravimeter and for proper connectionof airborne readings to the gravity network.
Figures 7 and 8 shows the gridded and acquired free air anomalies at flight altitude.
Colour agreement at cross lines indicate consistent data. Closer examination of themisfit in the 289 line intersections shows a 3.7 mGal RMS error indicating 2.6 mGal
average noise level.

Figure 7. Free air anomalies gridded and with artificial shade highlighting short wavelength features8

Figure 8. Free air anomalies at altitude

Comparison to EGM08 shows significant differences in many places amounting tomore than 130mGal over SE Mindanao, see Table 1 and figures 9 and 10.
Table 1. Data statistics and residuals

Unit: mGal
Airborne
GOCE R5-Direct residuals
EGM08 residuals

Mean
66.3
3.2
1.1

Std. Dev.
65.7
41.6
15.6

9

Minimum
-168
-126
-56

Maximum
293
191
135

Figure 9. EGM08 residuals at altitude

Figure 10. GOCE R5-Direct residuals at altitude.
10

Gravimetric geoid computation - principles
The PGM2014 is computed by the GRAVSOFT system, a set of Fortran routinesdeveloped by DTU-Space and Niels Bohr Institute, University of Copenhagen(Forsberg R, 2008). It forms the base of major recent geoid computation projects, suchas the joint Nordic “NKG” geoid models, undertaken as joint geoid model computationsof the Nordic and Baltic countries (R. Forsberg, D. Solheim & J. Kaminskis, 1996)
under the auspices of the Nordic Commission for Geodesy (NKG), as well as the
OSGM02 geoid model of the UK and Ireland (R. Forsberg et al., 2002), and several
national geoid models done from airborne surveys in recent years (Malaysia,
Mongolia, Indonesia and others).
The “remove-restore” technique was used in computing the geoid, where a sphericalharmonic earth geopotential model (EGM/GOCE combination) is used as a base. The
geoid is divided into three parts namely: the global contribution 𝑁𝑒𝑔𝑚 , a local gravity
derived component 𝑁2 , and a terrain part 𝑁3 .
𝑁𝑔𝑟𝑎𝑣 = 𝑁𝑒𝑔𝑚 + 𝑁2 + 𝑁3

(3)

The Remove – Restore Steps as illustrated by PAHLEVI, PANGASTUTI, SOFIA, and
KASENDA (2015) are:Remove Steps: Gravity measurement data is subtracted with global gravity anomaly
and surface correction, ∆𝑔𝑟𝑒𝑠 = ∆𝑔𝐹𝐴𝐴 − ∆𝑔𝑒𝑔𝑚/𝑔𝑜𝑐𝑒 − ∆𝑔𝑡𝑒𝑟𝑟𝑎𝑖𝑛
1. Extract gravity anomaly from the spherical harmonic model to produce∆𝑔𝑒𝑔𝑚/𝑔𝑜𝑐𝑒 (Long wavelength)
2. Extract gravity anomaly from SRTM data to produce ∆𝑔𝑡𝑒𝑟𝑟𝑎𝑖𝑛 (short –medium wavelength)
3. Subtract gravity anomalies ∆𝑔𝑒𝑔𝑚/𝑔𝑜𝑐𝑒 and ∆𝑔𝑡𝑒𝑟𝑟𝑎𝑖𝑛 from ∆𝑔𝐹𝐴𝐴 (measured
gravity) to get ∆𝑔𝑟𝑒𝑠
4. Apply FFT to ∆𝑔𝑟𝑒𝑠 to obtain geoid residuals ∆𝑁𝑟𝑒𝑠
Restore Steps: Geoid residuals summed with global geoid undulation and indirect
effect, resulting to geoid height (Undulation), 𝑁𝑔𝑒𝑜𝑖𝑑 = ∆𝑁𝑟𝑒𝑠 + ∆𝑁𝑒𝑔𝑚/𝑔𝑜𝑐𝑒 + ∆𝑁𝑡𝑒𝑟𝑟𝑎𝑖𝑛
1. Compute global geoid undulation from egm/goce ∆𝑁𝑒𝑔𝑚/𝑔𝑜𝑐𝑒
2. Compute the terrain part ∆𝑁𝑡𝑒𝑟𝑟𝑎𝑖𝑛 from SRTM data
3. Sum up ∆𝑁𝑟𝑒𝑠 (from remove steps), ∆𝑁𝑒𝑔𝑚/𝑔𝑜𝑐𝑒 , and ∆𝑁𝑡𝑒𝑟𝑟𝑎𝑖𝑛 to obtain
𝑁𝑔𝑒𝑜𝑖𝑑11

The spherical harmonic expression as a function of latitude, longitude and height isof the form:
𝐺𝑀

𝑅 𝑛

𝑛
𝑁( , , 𝑟) = 𝑅𝛾 ∑𝑁
𝑛=2 ( 𝑟 ) ∑𝑚=0(𝐶𝑛𝑚 𝑐𝑜𝑠𝑚 + 𝑆𝑛𝑚 𝑠𝑖𝑛𝑚 )𝑃𝑛𝑚 (𝑠𝑖𝑛 )(4)where 𝐺, 𝑀 𝑎𝑛𝑑 𝛾 are earth parameters. The EGM08/GOCE combination model used
involves more than 4 million coefficients 𝐶𝑛𝑚 𝑎𝑛𝑑 𝑆𝑛𝑚 derived from very large set ofglobal satellite data and regional (average) gravity data from all available sources,
both
open-fileandclassified,fordetailssee
http://earthinfo.nga.mil/GandG/wgs84/gravitymod/egm2008/index.html.
The EGM08 model is incorporating GRACE satellite data, which determines the error
spectrum of the EGM08 up to spherical harmonic degree 80 or so. New satellite datafrom the GOCE mission have recently been made available by the European Space
Agency, for details see www.esa.int/goce. The latest GOCE spherical harmonic model(“Direct” Release 5 model), complete to degree and order 280 was used to update
EGM08 field in the following way:-

EGM08 used unchanged in spherical harmonic orders 2-80, and from 200 up
GOCE R5 direct model used in band 90-180
A linear blending of the two models done in bands 80-90 and 180-200

The mixed spherical harmonic model (termed EGM08/GOCE) has been used tospherical harmonic degree N = 720, corresponding to a resolution of 15’ orapproximately 28 km. From the recent DTU geoid projects in France (Auvergne),
Malaysia, and Nepal, this resolution appears to be a good ”trade off” between the full
resolution of EGM08 (degree 2160) and the local gravity data. Because the fullresolution gravity data used in the construction of EGM08 is classified, there is no
good information on the quality of the errors in EGM08 at the high wavelengths in SE
Asia, and only 15’ mean gravity data are assumed to be underlying EGM08. All
spherical harmonic computations were done in a grid using the geocol17 program in
grid mode.
The terrain part of the computations were based on the RTM method, where
topography is referred to a mean elevation level, and only residuals relative to thislevel is taken into account. The mean elevation surface were derived from the SRTM
15” detailed model through a moving average filter with a resolution of approximately
20’ (37 km; slightly longer than the 15’ data resolution implied by spherical harmonic
degree 720, in order to have a more smooth residual gravity signal ∆𝑔2). The
difference in resolution between reference field and RTM is not a theoretical issue, as
the remove-restore method takes any “double accounted” topography into accountfully.
The method for the gravimetric geoid determination is spherical FFT with optimized
kernels. This is a variant of the classical geoid integral (“Stokes integral”), in which
there is a proper weighting of the long wavelengths from EGM08 and the shorter
wavelengths from the local gravity data. Mathematically it involves evaluating
convolution expressions of form
12𝑁2 = 𝑆𝑟𝑒𝑓 (∆𝜑, ∆ ) ∗ (∆𝑔2 (𝜑, )𝑠𝑖𝑛𝜑) = 𝐹 −1 (𝐹(𝑆𝑟𝑒𝑓 )𝐹(∆𝑔𝑠𝑖𝑛𝜑))(5)

Here 𝑆𝑟𝑒𝑓 is a modified “Stokes” kernel, ∆𝑔2 = ∆𝑔 − ∆𝑔𝑒𝑔𝑚 is the EGM08/GOCEreduced free-air gravity anomalies, and F is the 2-dimensional Fourier transform
operator. For details see references (R. Forsberg, D. Solheim & J. Kaminskis, 1996),(R. Forsberg et al., 2002) and (R. Forsberg & Olesen, 2010).
The geoid is computed on a grid of 0.025 x 0.025 resolution (corresponding toroughly 2.7 x 2.5 km grid). The area of computation is 04-22N and 112-128E,
covering the Kalayaan Islands of West Philippine Sea. Computations was based onleast squares collocation and Fast Fourier Transformation methods which involve
1440 x 1280 grid points corresponding to 100% zero padding. The data are griddedand downward continued by least squares collocation using the planar logarithmic
model. GRAVSOFT programs such as gpcol1, spfour, gcomb, geoip are involved in
the process. Figure 11 shows the PGM 2014 at 5m contour.
The final gravimetric geoid solution was computed by the following steps:- Subtraction of EGM08GOCE spatial reference field (in a 3-D “sandwich mode”)- RTM terrain reduction of surface gravimetry, after editing for outliers- RTM terrain reduction of airborne gravimetry- Reduction of DTU-10 satellite altimetry in ocean areas away from airborne data
- Downward continuation to the terrain level and gridding of all data by least-squares
collocation using a 1 x 1 moving-block scheme with 0.6 overlap borders- Spherical Fourier Transformation from gravity to geoid- Restore of RTM and EGM08GOCE effects on the geoid- Correction for the difference between quasigeoid and geoid (using a Bouguer anomaly
grid)- Shifting of the computed geoid by +80 cm to approximately fit to Manila tide gauge
datum

Figure 11. The Preliminary Philippine Geoid 2014 (PGM2014). Contour interval 5m
13

Data used and Quality Control for the geoid computation
The PGM 2014 is based on the following data:- Airborne gravity data
- Land gravity from NAMRIA, reformatted to GRAVSOFT and mildly edited- DTU10 global gravity anomalies from multi-mission satellite altimetry(Selected only in the open ocean area, away from the airborne gravity)- SRTM 15” DEM data for the region- EGM08 and GOCE RL5 satellite data

Some plots of the used and processed data are shown in the Figs. 12-14 below.

Figure 12. SRTM DEM data (left); low-pass filtered mean elevation surface (right), used as referencein RTM terrain reductions. Elevations in meters. Ocean depths are not used in the Philippines geoid
computations.

14

Figure 13. NAMRIA land gravity data and the airborne gravity data after terrain and EGM-reduction.
Some outliers were deleted in the geoid processing.

The available data from the airborne and surface sources were quality controlled
through plotting of the EGM08/GOCE and terrain reduction residuals, showing a few(< 1%) obvious surface gravity outliers, which were deleted in the final geoid
processing. The overall EGM/GOCE and RTM terrain “reduce” statistics for the data
are shown in Table 2. Overall this statistics is good, with relatively small bias andstandard deviation for all data sets.
Table 2. Statistics of remove steps in the PGM2014 computation (mGal)

Unit: mGal
NAMRIA edited land gravity data (1261pts)- above minus EGM08GOCE and RTM
Airborne gravity data (58515 pts)
- above minus EGM08GOCE and RTM
DTU10 altimetry gravity (28841 selected pts)
- above minus EGM08GOCE and RTM

15

Mean
8.7
0.9
1.5
1.6
-9.6
-0.3

Std.dev.
24.7
21.4
20.7
18.6
12.3
10.8

Figure 14. The used DTU-10 satellite gravity (left), and the collocation downward continued merged grid(right)Geoid processing results
The plots in the sequel shows the intermediate results of the final remove-restore
geoid processing (“geoid.gri”), computed with full 3-dimensional modelling, going via
the quasigeoid to the classical, final geoid.
The final geoid covers the region 4-22N, 112-128E, and has a resolution of 0.025 x
0.025. The airborne and surface gravity data were gridded by spatial least squares
collocation (gpcol1, using covariance parameters √C0 = 18 mgal, D = 6 km, T = 30km). A priori errors assumed were 2 mGal for both the airborne data and the surface
data (averaged in 0.025 blocks), and 5 mGal for DTU-10. The collocation downward
continuation was done in 1 x 1 blocks, with 0.6 overlaps.
For the spherical FFT transformation of gravity to geoid, 3 reference bands were used.
Fig. 15-16 below shows the primary data grids done in connection with the geoid
processing. The final geoid “restore” statistics is shown in Table 3.
Table 3. Statistics of the restore quantities on the geoid
Unit: meters
Mean
Reduced geoid (after spherical FFT)
0.00
RTM restore effects (computed by FFT)
0.00
Final gravimetric geoid statistics
39.06

16

Std.dev.
0.25
0.04
18.36

Min.
-1.61
-0.23
-9.02

Max.
2.88
0.74
76.43

Figure 15. Left: Reduced geoid (after spherical FFT transformation); right: RTM terrain effect on the geoid

Figure 16. Bouguer anomaly grid, derived from the reduced data. Used for the geoid-quasigeoid estimation
17

GPS/Levelling data comparison and fitted geoid
A set of 190 GPS data in ITRF2005 levelling benchmarks was used to compare withthe final geoid. These GPS data showed large errors relative to the geoid, with large
outliers in some regions, likely due to a combination of geodynamic effects, levellingand GPS errors. The rms fit is 0.5 m; it is therefore not possible to use these data forvalidation of the geoid. Fig. 17 shows the offset values, and the geoid correction
surface for a fitted geoid ph_geoid_fit (corrector surface gridded with 80 km correlation
length, and GPS-levelling apriori error of 10 cm). Fig. 18 shows the PGM2014
comparison to EGM2008; large improvements are seen, especially in the south. Fig.
19 shows a comparison to the DTU10-MSS, i.e. the MDT, shifted by -91 cm to take
into account the different reference systems.

Figure 17. Location of GPS/Levelling data. Colour show the correction surface for the fitted geoid18

Figure 18. Differences between the PGM2014 and EGM2008

19

Figure 19. Mean Dynamic Ocean Topography (MDT) from DTU-10 mean sea surface (shifted by
91cm)

The PGM2014 files are given as GRAVSOFT grids, and can be interpolated by the
GUI program grid_int (or the command line program “geoid”), provided to NAMRIA aspart of the computations, along with the general software and geoid job setups(Appendix 1).
Re-computation of the geoid
To further improve a geoid model, Professor Forsberg in his paper “Towards a cmgeoid in Malaysia” (R. Forsberg, 2003) recommends the following:-

Levelling networks must be carefully analysed for adjustment errors
Connections and antenna height errors of GPS data on benchmarks
must also be revisited and re-analysed
Erroneous points (geoid outliers) must be resurveyed by Levelling and
GPS
20

-

New GPS-fitted version of the geoid must be computed as newbatches of GPS-Levelling data, additional gravity surveys in major
cities and GPS users height problem reports comes in.

In 2015 NAMRIA started the re-computation of the PGM2014 with the help ofProfessor Forsberg. The densification of land gravity stations was conducted in somemajor cities of the country. Also, the GPS-levelling data was re-analysed, corrected
and outliers deleted.
Land Gravity Data
In this re-computation, the original airborne and satellite data processing results wereused, only the densified land gravity data (2214 points to date) were reprocessed andquality controlled. Figure 20 shows the new plots of the land and airborne gravity data.
Significant improvements can be seen in the land data (thicker dots). Most dots are ingreen, some yellows and light blue (i.e., 25-50 mGals difference in mountainous
areas).
GPS/Levelling Data
101 out of the 190 GPS/Levelling data on BMs remain (after cleaning up for outliers)
and used as validation points. After fitting the new GPS/Levelling, the RMS now is
0.054m with a minimum and maximum offset value of -0.124m and 0.169m
respectively. This improvement is mainly due to the removal of erroneous levelling
points. More points will be added to the GPS/Levelling data as the re-adjustment of
the levelling network progresses. Figure 21 shows the offset values and the new geoid
correction surface for the ITRF-fitted PGM2016.
Philippine Geoid Model 2016
The PGM2014 was re-computed to the new PGM2016 with an accuracy of 0.022m
using additional land gravity stations combined with the same airborne and satellite
gravity data. More land gravity data (up to 41,000) will be added from 2017 until 2020
in order to re-compute a new version and further refine the Philippine geoid. Figure 22
shows the new PGM2016 and its plotted differences with PGM2014 in figure 23. There
are differences in most parts of the country as big as 0.30m.

21Figure 20. Land gravity data after reprocessing and densification (2214points), plotted with theairborne data. Most differences are below 25mGals, although there are some points that exceeds
35mGals in mountainous regions.

22

Figure 21. New correction surface of the fitted geoid (fitgeoid-itrf2016)23

Figure 22. The new PGM2016. Contour interval 1m.

24

Figure 23. Differences between 2014 and the new 2016 geoid (ph_geoid files)

25

References
Dr. Bernhard Hofmann-Wellenhof, D. H. M. (2005). Physical Geodesy (pp. 46).
Forsberg, R. (2003). Towards a cm-geoid for Malaysia. Paper presented at the geoid computation
workshop in Kuala Lumpur.
Forsberg R, C. C. T. (2008). Overview manual for the GRAVSOFT Geodetic Gravity Field Modelling
Programs (Technical Report, 2nd ed.): DTU-Space.
Forsberg, R., D. Solheim & J. Kaminskis. (1996). Geoid of the Nordic and Baltic area from gravimetry
and satellite altimetry. Paper presented at the Symposium on Gravity, Geoid and Marine
Geodesy, Tokyo, 1996.
Forsberg, R., & Olesen, A. V. (2010). Airborne gravity field determination Sciences of Geodesy-I (pp.
83-104): Springer.
Forsberg, R., Strykowski, G., Iliffe, J., Ziebart, M., Cross, P., Tscherning, C. C., . . . Finch, O. (2002).
OSGM02: A new geoid model of the British Isles. Gravity and geoid, 3.
Harlan, R. B. (1968). Eotvos corrections for airborne gravimetry. Journal of Geophysical Research,73(14), 4675-4679.
http://igscb.jpl.nasa.gov/.
Kearsley, A. H. W. (1991). Evaluation of the geoid of the Philippines (Vol. 1): SAGRIC International.
Moritz, H. (1980). Advanced physical geodesy. Advances in Planetary Geology, 1.
Olesen, A. V. (2002). Improved airborne scalar gravimetry for regional gravity field mapping andgeoid determination. Faculty of Science, University of Copenhagen.
PAHLEVI, A., PANGASTUTI, D., SOFIA, N., & KASENDA, A. (2015). Determination of Gravimetric Geoid
Model in Sulawesi–Indonesia. Paper presented at the FIG Working Week.
Vanicek, P. (1991). Vertical datum and NAVD88. Surveying and Land Information Systems, 51(2), 8386.

26

APPENDIX 1 – Directory structure and GRAVSOFT geoid jobs
Directory Files

Comments

DATA

dtu-10_edited.fa
namria_edited.fa
airborne.fa

DTU-10 free-air anomaly data
NAMRIA edited free-air anomalies
Airborne free-air data

NAMRIAnamria.fa
gps-lev.dat
gps-lev-wgs84.dat

Original NAMRIA gravity data, reformatted from xls- files
Files with N from GPS-levellingDEM

ph_dem005.gri
ph_dem025.gri
ph_demref.gri

SRTM 15” file (30” used in Malaysia)
Averaged file for remote zone computations + geoid
RTM mean elevation surface, made with dem_ref.job

RD-TC

Remove steps (run “jobs” in sequence)
egm08goce5_720
Spherical harmonic coefficients for EGM08GOCE model
geocol17.job
Program to compute EGM08GOCE in grids for N and grd_surface.job
EGM and terrain reduction of NAMRIA data
rd_air.job
Reduction of airborne data, incl. terrain effect filteringrd_dtu10.job
Reduction of DTU-10
qc.jobGeogrid job to make difference NAMRIA minus airborne g

DOWNCONT Downward continuation and gridding
gpcol1.job
Collocation job, blockwise solutions
gpfit.job
Simple covariance estimation, use only as rough guide
GEOID

FFT and restore steps for geoid
bouguer.job
Job for making land Bouguer anomaly from downward
continued, reduced gravity grid file “downrd.gri”
n-zeta.job
Difference of geoid and quasigeoid, from DEM & Bouguer
n_rtm.job
RTM geoid terrain effect by prism integration (slow task!)
geoid.job
Composite geoid job, doing FFT, add terrain effect, add
EGM, apply corrections for quasigeoid -> “ph_geoid.gri”
geoip.job
Difference to GPS-levelling data (statistics and dif-file)

GPS-LEV Geoid fit files
fitgeoid.job

Auxillary:- fitgeoid_dif.dat:- fitgeoid_dif.gri:

Fitting of geoid to GPS-lev file with collocation. Always
plot the following files to judge fit of “ph_geoid_fit”:
point list with differences GPS-geoid
corrector surface for fitted geoid

grid_int_ph.exe
se-asia.bna

User-friendly interpolation program
Coastline file for “surfer” graphics

Note: to run GRAVSOFT on 64-bit Windows machines, check “job.bat” in GRAVSOFT
directory for the following statement in line 2, it should read: jobb64 <%1.job (job <%1.job for
Win32)

27