Modelling Geoidal Undulation and Reducing Ellipsoidal Elevations

I was once faced with difficult geo-referencing problem. We had just acquired a rather extensive high resolution aerial photography and LIDAR campaign in an area that only had datum transform parameters that were acurate to 5m and no adequate geoidal model that could be used for reduction of ellipsoidal elevations to local. All data was delivered referenced WGS84 (ITRF 2008) and LIDAR data was in processed, ungridded, point cloud form (60 billion points). We needed:

1. A model for vertical reduction from ellipsoid to local.


2. An automated means of applying the vertical corrections.


3. An automated means of gridding the LIDAR data.


4. An automated means of re-projecting and applying higher precision horizontal transformations.

I determined the ellipsoid separation model by obtaining information from survey control stations where the local stations had also been occupied by GPS and those observations were available. These were few and far between and an even distribution across the survey area was highly desirable. In the end I managed to find 9 suitable primary control points and 4 secondary points that could be used for densifying the model.

Fig 1: Shows primary control points (large red cross) and 4 additional secondary control points used for densification purposes from benchmarks of various topographic surveys (smaller red cross).

A polynomial regression surface was fitted to the separation values from the control points using ArcGIS Geostatistical Analyst software. The regression fit was also measured. Fig 2: Local Polynomial Regression Parameters in Geostatistical Analyst.

Fig 3: Geostatistical Analyst showing the regression function and fit to the input control data. The model fit had an average error of 6mm.

The regression surface was output at a suitably high resolution to a raster dataset. To automate the reduction and gridding of the point cloud data I turned to my trusty Raster Engine which I modified so that the binning function could also accept a "correction" surface. I also added an interpolation method so that no "stepping" could be observed at the boundaries of the regression surface's cells on the final gridded LIDAR data. The LIDAR data had been delivered as several hundred tiles (for ease of data handling). The Raster Engine processed all tiles (60 billion points) with binning and vertical reduction in a matter of hours (without parallel processing).Fig 4: 3D Visualisation of 3 processed LIDAR tiles in Priority Area 2. No “stepping” effects can be seen even under a vertical exaggeration of 25x – again validating the success of the sub-pixel interpolation of the separation grid.

A number of secondary survey points (34) were reserved (i.e. not part of the regression model input) as a means of checking the validity of the model and the resultant reduction. During comparison these showed an average error of 2cm.

Improving the horizontal accuracy was achieved by comparison of the reprojected high resolution aerial photography with surveyed as-built features. A conformal transformation in grid space was derived for the entire survey area. Reprojection and application of the conformal transformation was automated using the ESRI re-projection engine. The entire dataset was processed in one weekend on a single workstation.