How to present georeferencing measures: RMSE, skew, anisotropy

[mradamcox]

I've been working on the addition of a few quantitative measures that can be calculated based on a layer's georeferencing information. Specifically, the GCPs + transformation. The idea behind adding these measures that 1) during georeferencing, they can provide immediate feedback to users, and 2) they can be attached to each layer, and sorting for low values can guide targeted georeferencing improvement efforts.

The three metrics I've added to the database are:

• RMSE - Root mean square error (see [Request] Add RMSE calculation to GCPs #283)
• Skew - "Tiltedness" or degrees off square.
• Anisotropy - The ratio of the x-scale factor to the y-scale factor in the result image.

(A fourth one, distribution of GCPs, #284, is something I had thought would be useful in the past, but is somewhat superseded by the combination of these others. Still something to consider in the future though.)

Much credit to @danvk 's work in mapsnap, where he put together the calculations for all of these metrics, and I borrowed some of his functions to get started with this.

The initial addition of a layers list that displays these measures looks like this:

More info on each measure:

RMSE

• This is the average of all Root Mean Square Error values across the group of GCPs. A GCP's error is the distance (in meters) between its predicted location (where the transformation algorithm wants to place it) and its target location, where the user actually placed the point.
• To illustrate, the red lines here are the distances between the predicted locations (shown in red) and the target locations (the black X's):

When is this useful?

• A high RMSE, say over 2.5 meters, will indicate that the set of GCPs has a low internal agreement. It's likely that one GCP is way off.

When is this useless?

• With the minimum number of GCPs needed for a transformation (2 for helmert, or 3 for 1st order polynomial (our default)), RMSE will always be 0, because all GCPs will fit the model no matter where they are placed.
• Using the Thin Plate Spline transformation (RMSE will always be 0 because this transformation distorts the image to fit every GCP exactly.

When is this misleading?

RMSE has nothing to do with where the GCPs are actually placed geographically, or, indeed, with how much the page is distorted. You can make a very rotated/skewed result and have 0 RMSE:

Skew

The skew of a layer, is the amount it is "tilted", not rotated, one direction or another. It is measured in degrees, and very easy to conceptualize with these large-scale maps. For example (disregard the basemap here):

As you can see, the skew to the left results in a measurement of -18 degrees.

When is this useful?

• Almost always, because all city block should remain square after the layer is georeferenced. A skewed layer will have blocks as parallelograms.

When is this useless?

• With helmert transformation, skew will always be 0

When is this misleading?

• As above, skew has nothing to do with where the map is actually placed.
• Historical maps are not perfect, and the scanning process may have introduced warping or other distortion. So it is entirely possible that the best georeferencing fit will have some skew in the end, especially when placed within context of nearby sheets.

Anisotropy

This is the degree to which a map is stretched along its vertical or horizontal axes, and it is a calculated as the ratio of the x-scale factor / the y-scale factor. So, a map stretched wider than it is tall will have anisotropy > 1, and < 1 if it is stretched taller than it is wide. Here's an example with anisotropy ~2, meaning the x scale is about twice the size of the y scale:

When is this useful?

• Almost always, because these maps are drawn to scale

When is this useless?

• With helmert transformation, anisotropy will always be 0

When is this misleading?

• As above, anisotropy has nothing to do with where the map is actually placed.

Questions & ideas

Two main issues I see here are: 1) the names of these measures are pretty esoteric and unintuitive, and 2) these values should be normalized so they are sortable across the database as easy as possible.

• RSME
  • Renaming ideas: GCP error, Error, Fit error,
    • Fit or "GCP Fit", could work but low RMSE is good, and a low "fit" value seems unintuitive. May need to inverse/normalize?
• Skew
  • Kinda don't think this needs to be renamed, it's pretty intuitive
  • Currently the values go from negative (bad) to 0 (good) and then positive (bad). Propose using absolute values instead, so sort descending will show the most skew, in either direction, first.
  • Is there any use at all in storing the real value? Useful to know how many pages are skewed to the right vs. left? (I don't think so...)
• Anisotropy
  • Renaming ideas: Stretch, Squish??
  • Currently the values go from < 1 (larger y scale, bad) to 1 (equal scales, good) to > 1 (larger x scale, bad). Propose inverting the < 1 values so that .25 becomes 4, .5 becomes 2, etc.
    • In other words: stretch = anisotropy if anisotropy >= 1 else 1/anisotropy
  • Is there any use in storing the raw values? Useful to know whether a page is stretched vertically or horizontally?

[tannewt]

I'm not sure I'd show them all as columns. Instead, maybe a simple warning icon should show when one of these is out of range. Then clicking or hovering could give more detail.

I've definitely wished I could ignore gcps and just scale and rotate on the map. Sometimes finding specific features is hard but aligning road/feature shapes could be easier.

[mradamcox]

@tannewt Good points, and with regard to scaling/rotating, see #403, where there is some more info about how you can georeference with only two points. I will enable that on your account, and, more than likely, just make it the default soon. That will help you a lot.

As for icons on layers with error, that's an interesting idea, maybe these numbers should be aggregated into a single score and that could drive a tag like that. That would also address situations where some stats are useless and others are more helpful. If transformation is polynomial 1 and there are 3 gcps, ignore RMSE, for example.

[mradamcox]

Interesting thing I realized after deploying this on the production database (still behind a "staff" login though), storing RSME as meters, i.e. as calculated on the geo side of the process, will bias large RMSE to large geographic coverage maps, like Key Maps:

Not a show-stopper, because you can filter to only show Main Content

...but something to keep in mind. It's possible to reverse this calculation though, to show it in pixels on the source document rather than distances on the map side.

[mrkrause]

Thanks! I've been playing with this for a little while. I sometimes accidentally drag a GCP while trying to pan the map and the RMSE lines are super-helpful for fixing that. It'd be even better if there was some way to exaggerate them (a 5m error can hide pretty easily on a large-scale map) or perhaps color the points differently if their residuals are above some threshold?

A few more comments.

1. Is the RMSE calculation in utils/gcp.py correct? It looks like it's calculating the error per-coordinate (i.e., along the x- and y-axes separately), rather than pointwise. They're related, but the latter is larger; it is also what is shown via the red lines in the georeferencing interface show.

2. I'm concerned that the RMSE number may discourage people from adding GCPs. Going from 3->4 will usually jump the number up from 0 to [something] because the fit is no longer exactly determined. Less obviously, the RMSE can also jump up when you add a 5th or 6th point, because in-sample RMSE is "optimistic" (biased downwards), and strongly so when calculated from a handful of points. Thus, I suggest:

• Showing n/a until there are 4 (for Polynomial 1) GCPs or 3 (Helmert) so that the RMSE is always sort of reasonable.
• Bias-correcting the RMSE estimates when there are more points.

3. The RMSE indicates how consistent the points are with this type of transform, but the more relevant question is how well the transform generalizes to the rest of the map. This can be tested via Leave-One-Out: you re-fit the transform using N-1 points and then ask how well it predicts the held-out point. (In practice, you don't even have to do this--there's a trick). However, this tends to go pear-shaped if there are only a few points or one is driving the fit. I'm not sure if that's a bug or a feature!

Do you want a PR for any/all of this? Also, my apologies if any of this is naive! (My day job involves a lot of stats but not GIS).