Georeferencing a Scanned Map From Coordinates Alone

Key Takeaways

• The task today: Georeferencing a scanned map means reading coordinates off its margin and fitting a transform to them, with residuals and a holdout error showing any control point that disagrees with the rest. Stratigraphic Amenity, an open-source MCP server (a tool server that AI agents call), does this from a coordinate system, four or more control points and the map’s pixel rectangle, with no map-detection model or knowledge base.
• The problem: Those checks compare the points with each other, so an error every point shares passes clean: a datum left off the CRS string (read as WGS84, around 200 m from a NAD 27 grid in the lower 48), a wrong UTM zone that moves the map six degrees of longitude, or short tick labels that squash it. In this mode the tool cannot see the map collar, so it cannot catch any of them.
• What the tool does: geomap_georeference returns EPSG:4326 bounds, the six affine coefficients, residuals in metres and warnings for flipped or rotated axes, as local math that stores nothing, so a second fit costs nothing extra. What it cannot judge is place, so the reliable check is a second witness: fit again on the latitude and longitude printed at the corners, compare the two boxes, and test the implied pixel size against the map’s scale.

TL;DR

Georeferencing a scanned geologic map from its tick coordinates needs no detector: Stratigraphic Amenity, an open-source MCP server for AI agents, turns a coordinate system, control points and the pixel rectangle into bounds, residuals and warnings. A clean result only means the control points agree with each other, so errors they all share, such as a missing datum, a wrong UTM zone or short tick labels, leave the map in the wrong place without a warning. Writing the datum out, fitting a second time on the corner coordinates and checking the pixel size against the map’s scale catch what the residuals cannot.

How AI for Geology Georeferences a Map It Never Sees

Give an agent four tick readings from the margin of a scanned geologic map. Add the coordinate system printed in the map collar and the pixel rectangle inside the neatline. Stratigraphic Amenity will then georeference the sheet without ever looking at it. Called with only crs, gcps and pixel_extent, geomap_georeference is deterministic local projection math. It needs no detector, no knowledge base and no registered map. It returns a georef/v1 payload with six affine coefficients, EPSG:4326 bounds, residuals and warnings. With no map to attach it to, it saves nothing.

Suppose that payload comes back clean. Four control points, fit_method set to affine-least-squares, a small residual_m, a holdout_error in the same range, an empty warnings list. The map can still sit a couple of hundred metres from where it belongs. Nothing in the payload will say so.

That is not a bug in the math. It is what math does when it runs on its own.

What a Residual Tells an AI for Geology Agent

Every quality number in that payload compares the control points with each other. residual is the RMS misfit of the affine in the grid’s own units. residual_m turns each point’s misfit into a geodesic distance in metres. holdout_error refits the transform without each point in turn and reports the worst miss. All three answer one question: do these points agree with one another under an affine model?

It is a good question. A mistyped easting, a tick read one grid line over, a scan that isn’t flat: each makes one point disagree. With four or more points, the disagreement shows. With two or three, the fit is exact, and the tool says so with a residual_not_diagnostic warning. That is why georeferencing from ground control points should start at four.

But agreement is not truth. An error that every point shares moves them all together, and points that move together still agree. Geologic maps offer at least three such errors, and all three are common.

Three Georeferencing Mistakes the Control Points Share

The datum. A typical USGS credit note mentions “1000-meter Universal Transverse Mercator grid ticks, zone 15” in one sentence. It names the “1927 North American Datum” in another, and that second sentence is easy to leave behind. Stratigraphic Amenity resolves a bare UTM string to WGS84, so “UTM zone 15” becomes EPSG:32615. Many older North American sheets have a NAD 27 grid. On those, every control point carries the same wrong label. According to USGS, UTM values change by around 200 metres between NAD 27 and NAD 83 in the lower 48. WGS84 sits close to NAD 83 at this scale. A NAD 27 grid read as WGS84 therefore lands about as far off. The residuals barely notice. The affine is fitted in grid coordinates before the datum enters. The datum only decides where the finished map is placed.

The zone. Write zone 16 where the map collar says 15. Every point then moves exactly six degrees of longitude east, in the grid’s own datum. A UTM zone is the same projection with its central meridian six degrees over. So the distances between the points do not change, and neither do residual_m or holdout_error. The bounds are still valid longitudes and latitudes. They are just in another place.

The short labels. On a USGS sheet, a tick near the corner may carry a full value, for example 249000mE. Its neighbour may read only 248, short for 248000mE. Northings are labelled the same way, with the leading digits printed small. Suppose a reader, human or model, copies the short northings from every tick. The tool then gets a map that is internally flawless and squashed onto the equator.

What the Tool Checks, and What It Does Not

This is not a complaint that the tool checks nothing. It checks the shape of the transform. If a, the x scale, is zero or negative, it flags a nonstandard axis direction. It does the same when e, the y scale, is zero or positive. Both checks assume an unrotated, north-up scan. If the rotation or shear terms exceed a tenth of the scale, it asks you to verify the map’s orientation. Those warnings catch the classic swap of easting and northing. A swapped map looks mirrored or turned ninety degrees, and the checks see that.

What they cannot catch is an error of size or of place. The map may be the right way up and in the wrong place, or squashed a thousandfold in one direction. The tool cannot see the map collar in this mode. So it cannot read the datum, the zone or the typographic convention. The tool defaults a datum-less UTM string to WGS84 instead of refusing it, so the datum is the caller’s job. That default is a guess, reasonable for a recent map and poor for a 1970s quadrangle. The defence is to write it out: “NAD27 UTM zone 15” resolves to EPSG:26715.

Checking the Datum and Zone With a Second Witness

If the fit can only cross-examine its own points, an AI for geology agent has to bring the check in from outside them. Maps are generous here, because many carry two independent coordinate systems. A quadrangle’s neatline corners are printed in latitude and longitude; its margin carries projected grid ticks. So fit twice. Fit once on the grid ticks, with the grid’s CRS. Fit again on the corner coordinates, with the geographic CRS of the same datum (EPSG:4267 for NAD 27). Then compare the two boxes.

Each fit’s mistakes are its own. A wrong zone moves only the grid fit. Short labels damage only the ticks they were copied from. A wrong datum moves the two fits by different amounts. The same USGS page gives 10 to 100 metres of shift for latitude and longitude, but around 200 for UTM values. So the boxes stop agreeing. A georeference-only call registers nothing and stores nothing. As a result, the second fit costs no more than the first.

Check the Scale the Coefficients Imply

There is also a check that needs no second fit. The coefficients a and e are the size of a pixel on the ground. A geologist can predict them from the other direction. Multiply the scale denominator by 25.4 mm and divide by the scan’s dots per inch. For a 1:24,000 sheet scanned at 300 dpi, expect about 2 metres per pixel; for a 1:250,000 sheet, about 21. A fit whose e is a thousandth the size of its a has a digit problem. No warning will mention it.

Where the Weight Sits

None of this makes georeferencing without a detector weaker than the full georeferencing pipeline. In the full pipeline the detector finds the map area; it never reads a datum either. The georeference-only route just makes plain where the weight sits. It sits on a few numbers read from a margin. It also sits on the one line in the map collar that says what they mean.

How much that matters depends on the question. For a regional earthquake search on a 1:250,000 sheet, two hundred metres is noise. For a drill collar, a claim boundary or a contact you mean to walk to, it is the whole answer. A clean residual tells you the points agree. Whether they agree about the right place is a question for a second witness.

Finding that witness is work for whoever reads the margin. The fit itself is not something you have to build. You need no projection code, no least-squares solver and no coordinate plumbing of your own. Stratigraphic Amenity is open source, and its georeference-only route installs without the detector or the knowledge base. The server an AI for geology agent would call can run on a laptop next to a stack of scanned sheets.

Before trusting a clean result from it, though, it is worth seeing exactly what it inspects. The judgement it passes on a transform is short enough to read in one sitting. What it leaves out is as easy to see as what it checks.

The Checks That Read Only Four Numbers

The warnings that judge the transform itself, rather than the number of control points, come from one short function. It runs after the affine is fitted:

def _affine_warnings(affine: AffineTransform) -> list[str]:
    warnings: list[str] = []
    if affine.a <= 0:
        warnings.append(
            "nonstandard_x_direction:a<=0; expected only for an unrotated north-up scan"
        )
    if affine.e >= 0:
        warnings.append(
            "nonstandard_y_direction:e>=0; expected only for an unrotated north-up scan"
        )
    scale = max(abs(affine.a), abs(affine.e))
    cross_term = max(abs(affine.b), abs(affine.d))
    if scale and cross_term > 0.1 * scale:
        warnings.append(
            "substantial_rotation_or_shear:cross-term ratio exceeds 0.1; verify map orientation"
        )
    return warnings

Together, the three checks read four of the six coefficients and nothing else. The offsets c and f, which place the map, are never read. Two checks look at signs, one at the ratio of the rotation terms to the scale terms. No check looks at how large a or e is, or whether they match each other. Nothing checks where the map lands either: the bounds only have to be legal longitudes and latitudes (knowledge/bounds.py:37-46).

The residuals have the same blind spot for a different reason. fit_affine fits the transform from the control points alone (transform.py:129-163). The resolved CRS only enters afterwards, to reproject the corners and to express each misfit in metres. A wrong datum, a wrong zone or a column of short labels changes no sign and makes no point disagree. So _affine_warnings returns an empty list. The datum itself comes from georef/crs.py:64, where a UTM string with no datum falls back to "wgs84".

Run It From an AI for Geology Agent

To try it, clone the Stratigraphic Amenity repository and install only the MCP and georeferencing extras. The detector extras are not needed for this route:

uv sync --extra mcp --extra geo

Then give your MCP host a server entry with "command": "stratigraphic-amenity-mcp". The host launches the server and talks to it over stdio. Call geomap_georeference with no map_id, and write the datum out in full:

{
  "crs": "NAD27 UTM zone 15",
  "gcps": [
    {"pixel_x": 412, "pixel_y": 388, "world_x": 661000, "world_y": 4899000},
    {"pixel_x": 2380, "pixel_y": 388, "world_x": 665000, "world_y": 4899000},
    {"pixel_x": 412, "pixel_y": 2356, "world_x": 661000, "world_y": 4895000},
    {"pixel_x": 2380, "pixel_y": 2356, "world_x": 665000, "world_y": 4895000}
  ],
  "pixel_extent": [150, 120, 5250, 6720]
}

The reply is a georef/v1 payload, with crs resolved to EPSG:26715 and the six affine coefficients. Here a is about 2.03 and e about -2.03 metres per pixel. That is what a 1:24,000 sheet scanned at 300 dpi should give. The payload also carries EPSG:4326 bounds, residual with its residual_units, residual_m, fit_method, gcp_pixel_errors, holdout_error and warnings. There is no georef_uri, because nothing was stored. For the second witness, call it again with "crs": "EPSG:4267" and the four neatline corners in decimal degrees. Put longitude in world_x and latitude in world_y. Keep the same pixel_extent, so the two boxes cover the same part of the scan. Then compare the two bounds. In Python, georeference_bounds(crs=..., gcps=..., pixel_extent=...) does the same with only the geo extra. Its result’s pixel_to_lonlat checks a single corner.

FAQs

Can AI for geology tools georeference a map without a detector?

Yes. Stratigraphic Amenity’s geomap_georeference takes a coordinate system, at least two ground control points and the map’s pixel rectangle. It needs no registered map, detector or knowledge base. It runs local projection math, returns longitude and latitude bounds with residuals and warnings, and saves nothing. Use four or more points so the residuals mean something.

Why does AI for geology georeferencing pass with the wrong datum?

Residuals compare the control points with each other. With a wrong datum, a wrong UTM zone or short tick labels, every point shares the same error. So the points still agree, and the residual stays small. The affine is fitted in grid coordinates before the datum enters. The datum only decides where the finished map lands.

What datum does Stratigraphic Amenity assume for a UTM string?

WGS84, when the string names no datum: “UTM zone 15” becomes EPSG:32615. If the map’s grid is on NAD 27, write that out. “NAD27 UTM zone 15” resolves to EPSG:26715. According to USGS, UTM values change by around 200 metres between NAD 27 and NAD 83 in the lower 48.

How can I check that an AI for geology georeference is in the right place?

Use a second witness. Fit again on the latitude and longitude printed at the neatline corners, in the geographic CRS of the same datum. Then compare the two bounding boxes. Also check the pixel size: scale denominator times 25.4 mm, divided by the scan’s dpi. That gives about 2 metres for a 1:24,000 sheet at 300 dpi.

Try It With Geocluster

Stratigraphic Amenity is one of the Geocluster tools, MCP servers that let AI agents work with geology data and maps.

Once the bounds pass both checks, they are ready to look up earthquakes, faults and minerals for a bounding box.