# How to (vertically) register a DEM to GPS readings?

I don't know if what I ask is possible or not but let's try...

I have two datasets:

• a DEM as a point cloud in a relative coordinate system, constructed with 3D algorithm structure-from-motion.

• 20 ground control points with absolute elevations values (between 130 and 140 meters) taken with a GPS.

• 20 points of the point cloud representing the GCPs in the relative coordinate system (georeferencing targets visible on the point cloud).

I would like to constrain the elevation of my point cloud DEM with the GPS ground truth elevation, in order to have a DEM with absolute elevation calibrated to my study area.

At the moment, the study area is about 1 hectare and with slopes around max. 20% but the aim is to export the methodology in river catchment with a size of some kilometers-square and with max. slopes around 60%.

Any help is welcome, preferably in ArcGIS or R but I can deal with other softwares if necessary!

Thanks for the help.

• What is the extent of your study area and what are the typical slopes in the DEM? Commented Jun 26, 2013 at 16:08
• Thanks. Assuming typical GPS location error around 1 - 5 meters and with slopes of 20+%, it becomes important to account for the GPS location accuracy. (A 20% slope translates to a potential 0.2 - 1.0 meter error in the elevation associated with each GPS reading, with the error varying with slope: thus GPS readings at low-slope areas provide better information than readings in high-slope areas.) With lower slopes, the location accuracy is less of an issue. Commented Jun 26, 2013 at 20:51

Finally, I found the solution to my question...

So basically, to georeference a point cloud based on GPS measurements, we need to find a linear transformation (rotation matrix, translation vector and scaling factor) between the relative coordinates of the GCPs visible on the point cloud and the world coordinates of the GPS measurements. And apply this transformation to each point of the point cloud...

To do that, the procrustes() method, available both in R and MATLAB, can be used. It returns the 3 parameters needed for the linear transformation to be done. Then the coordinates of each point of the point cloud can be transformed into world coordinates based on this equation:

World Coordinates of the point = Scaling * Rotation * Relative Coordinates of the point + Translation

Hope this is clear!

• This otherwise interesting solution is deficient in two respects. First, it is usually the case that better registration is obtained by considering an arbitrary affine transformation: you cannot expect the GPS and DEM coordinates to agree perfectly up to just a rotation and translation. Second, this solution is inferior to those that explicitly account for the heteroscedasticity (variable positional error) introduced by the varying slopes. If you're referring to a three dimensional solution, that's even worse, because the GPS error is much greater vertically than horizontally. Commented Sep 2, 2013 at 15:07
• Could you please explain a bit what do you mean by "a solution that account for the heteroscedasticity introduced by varying slopes"? Commented Sep 11, 2013 at 12:31
• When the slope is steep, a small change in horizontal position leads to a relatively large change in elevation. When the slope is zero, changing the horizontal position leaves the elevation unchanged. Therefore steep areas need to be registered accurately whereas horizontal flat areas provide no useful information at all. There are statistical methods to handle this varying sensitivity, such as weighted least squares, and they are straightforward to apply with statistical software like `R` or with suitable Python packages. Commented Sep 11, 2013 at 14:06
• If you have the data, applying the DEM corrections to the GPS pseudo-ranges (rather than the output positions), as is usual for DGPS, might give a better result. Commented Nov 8, 2013 at 23:09

Using Kriging there are various methods to constrain your interpolation based on other data:Regression Kriging/Kriging with external drift.

http://spatial-analyst.net/book/Regressionkriging

on the diffrence between methods: http://www.itc.nl/library/Papers_2003/misca/hengl_comparison.pdf

For a similar problem I have used the Kriging with External Drift method using the g-stat package in R. In R using the gstat package, something like:

``````Krig_res <- krige(val~constraint,locations,data,model=variogram)
``````

where the constraint is a collumn in your 'data' dataframe.

hope it helps Annette

• This is an interesting idea. It's not quite clear how it would be used to change the values of an existing DEM. Could you clarify that point? Commented Jun 26, 2013 at 17:48