qgSurf, a plugin for QuantumGIS, an open-source and free GIS software, adds this functionality in its new version (0.2.1, installable from QGIS plugin manager or downloadable from QGis plugin repository). Previously its aim was to calculate the intersections between a geological plane and a DEM. The calculation of the best-fit plane given points on a DEM is a somewhat inverse procedure so that they complement each other.
To see a short help for the plugin, look at: https://bitbucket.org/mauroalberti/qgsurf/wiki/Help
The basis of the algorithm is the application of singular value decomposition to derive the eigenvectors of a set of measures. See for instance the discussion in Best fit plane algorithms why different results?
Obviously we have to dispose of a DEM, e.g., Aster or SRTM or others.
A very interesting QGis plugin, OpenLayers, allows to load GoogleEarth or Bing maps in a project, in order to use them as backgrounds for our analysis. However, this requires the on-the-fly reprojection of layers while qgSurf does not support it for DEMs, so we need to set as the project projection that of the DEM: UTM, Lambert, or whatsoever.
The processing sequence is:
a) load in the QGis project the required DEM(s) layers and whatsoever vector or image layers needed for your analysis
b) use "Get current raster layers" in qgSurf plugin: this will allow the plugin to know which raster layers are currently loaded
c) from the "Use DEM" combo box, choose the required DEM and make sure that the QGis project and the DEM have the same projection
d) from the "Best-fit-plane calculation", press "Define points in map": this will allow you to define in the canvas at least three, and possibly more, points, whose coordinates will be listed in the plugin widget.
e) with at least three points defined, you can calculate the best-fit plane by pressing "Calculate best-fit-plane": a message box will report the dip direction and dip angle of the calculated plane
f) you can add even more points and again calculate the best-fit plane; otherwise, if you want to start a new analysis on the same DEM, go to e), or if you want to use another DEM, go to c) if it is already loaded in the project, or load it in the project and then go to b)
Example
Mt. Raparo, Basilicata, limestones of the Panormide Complex.
Source DEM is Aster, projected in WGS84-UTM33N. Base map is Bing, loaded in QGis via openLayers plugin. Project CRS is set to the same of the source DEM, i.e. WGS84-UTM33N, otherwise the plugin wouldn't work.
Two level are recognizable in the aerial image, and a few points are drawn to derive their attitude using the new plugin functionality. The results are 336.2 / 8.7 and 337.8 / 13.0 (dip direction nd dip angle).
The results are checked with field measurement made by myself and M.C. Lapenta during field work in 1991: 320 / 5, 290 / 15 and 345 / 15. The geological plane attitudes are congruent both with themselves and with field measurements (note that while the field measurements refer to magnetic North, GIS-derived values refer to map top for the current projection).
