Tag Archives: georectification

The title says it all: how to georectify archaeological field drawings in a GIS system when you only have two ground control points (GCP).

This is a problem that I can imagine surfaces also in other projects than only those I’m involved with. In general, when doing field drawings on an archaeological excavation, two known points, or GCPs, are enough to position the drawing exactly within the local coordinate system. The extent covered by a single drawing is usually some square meters, and even in the very unlikely case the local coordinate system would be in some other projection than purely rectangular, the errors caused by importing the plan to a GIS using only two points are minimal compared to the errors made while drawing anyway, that to use more than two points is only necessitated by the needs of a large drawing to have other fixed references to cover the whole drawing area.

On the other hand, a GIS system, which otherwise is very good in organizing the measurement data and the drawings of an excavation, often has a georeferencing system that is based on a minimum of three GCP’s. This is (probably) because the most commonly used forms of affine transformations in principle call for three GCP’s, and in many cases, the GIS includes options to provide way much more GCP’s than three. The reason for this is twofold: firstly, the use of more than the necessary amount of GCP’s allows for the use of statistical means to estimate errors; secondly, the use of more points allows the use of other than linear transformations, ranging from second grade polynomials to “rubber sheeting”, very useful for historical and otherwise distorted maps, that can be georeferenced using an array of points distributed evenly over the whole map area.

However, the problem discussed in this post is at the other end of the scale: exact drawings covering a small area, often in scales between 1:10 and 1:50 – although 1:1 drawings of especially interesting details are not unknown.

Let us assume the we have such a drawing, with two known control points with the coordinates in the local coordinate system. The problem is, how to import this drawing into a GIS? A scanned drawing can be imported, and a georectifying interface, like g.gui.gcp in GRASS, can be used to mark the GCP’s. In GRASS, this can be made using the imported drawing so, that the GCP’s are marked on the original drawing, and the corresponding coordinates for the local coordinate system entered on the list of GCP’s visible on the window. However, to perform the rectification, one needs the enter at least three coordinates. The third coordinate can be simulated, for example drawing additional elements both in the original drawing location as well as the target location, i.e. the location where the imported drawing is to be placed. However, this process is rather complicated and slow. (If the drawing capabilities were at the level of AutoCAD, it would be much speedier, as one could add elements based on the geography of the existing elements, i.e. the GCP’s.)

One, simple solution is to use an Octave/MATLAB script like below to calculate a third common point between the source and target location:

For this script, you need the Octave image library, which contains the cp2tform and tformfwd functions.

The first line puts the coordinates of the GCP’s on the original drawing into the variable movingpts. The second line puts the GCP coordinates in the local, target coordinate system to targetpts. The third line adds a third point in the drawing to newpoints– this coordinate can be anything, as it does not refer to any existing point on the drawing itself. The fourth line created the transformation matrix; important here is the third parameter 'nonreflective similarity', which defines the transformation type. The ‘nonreflective similarity’ in practice means, that the only allowed operations are rotation, scaling and translation, which is exactly what we need for the archaeological field drawings.

On the fifth line, the coordinate for the third point is transformed to the target coordinates using the defined transformation matrix. The result in this case is 5150.766081868117 5006.89610151436.

And not we have the necessary coordinates for defining the third GCP: 1500 1500; 5150.77 5006.90.

A careful person integrates this whole process as a part of an Emacs Org-mode file using Babel for reproducible results. It is always good to document where you got the numbers you were using.