How to Batch Georeference Maps and Turn Pixels Into Coordinates in Python

Key Takeaways

• The task today: To georeference scanned maps and look up the earthquakes and faults inside them, an AI agent calls Stratigraphic Amenity, an open-source MCP server (a tool server that AI agents call). It calls one tool at a time, which suits a model that chooses each next step from what it reads. The same repository is also a Python package with those functions underneath, so a script can run the identical fit and query with no agent.
• The problem: A fixed batch job, such as three hundred sheets whose grid ticks are already read, has no steps left to choose. There the server adds limits: it answers in boxes rather than points, writes a new file for every knowledge query and keeps a cache you cannot switch off. Skipping it also drops checks we added after watching agents misread results in our own tests: record totals rather than a count of result items, a warning for each provider (one lookup, such as earthquake history) skipped by default, and a summary that repeats every warning.
• What the tool does: In Python, georeference_bounds returns a georeference whose pixel_to_lonlat turns any pixel into longitude and latitude. The knowledge service (the Python class that runs the lookups) returns a bundle: one result object with an item per provider, each carrying its record count and a truncated flag, plus the warnings and a trace, the log of which providers ran and which were skipped. Skipped providers appear only in that trace, with no warning, and nothing totals the records, so the script has to read the trace and the warnings before it trusts an empty answer.

TL;DR

Want to georeference scanned geologic maps and look up their earthquakes and faults from a Python script instead of through an AI agent? Stratigraphic Amenity, an open-source MCP server for AI agents, is also a Python package, and calling its functions directly adds per-pixel coordinates, a cache you can switch off and your own data sources. What a direct call loses are the checks we built into the server after watching agents misread results, so the script has to count records rather than result items, read every warning and check which lookups were skipped before it trusts an empty answer.

When a Geologic AI Job Needs No Agent

Picture a geologic AI job with nothing left to decide. There is a folder of three hundred scanned quadrangles and a spreadsheet with one row per sheet. Each row holds the coordinate system printed in the map collar, four grid ticks with their pixel positions and the pixel rectangle of the map face. Someone has already read the margins. Now you want the EPSG:4326 footprint of every sheet, and the recorded earthquakes and active faults inside each one.

Nothing in that job needs judgement while it runs. The readings are done and the steps never change. So you skip the MCP server and import Stratigraphic Amenity as a Python package. Then you write a loop: georeference_bounds on each row, then KnowledgeService.query_bounds on the bounds it returns. That is the right call. The rest of this post is about what leaves with the server when you make it.

What Does the Python SDK Share With the MCP Server?

The server and the package are not two implementations. The MCP layer describes itself as a “thin protocol-independent adapter over existing local SDK services”, and the description is accurate. geomap_georeference builds GroundControlPoint objects and calls georeference_bounds. geomap_query_knowledge builds a KnowledgeRequest and passes it to KnowledgeService.query. Same affine fit, same reprojection, same providers, same warnings at the source.

So the choice between them is not a choice of engine. It is a choice of reader. Behind the server, the reader is a model. It sees a tool result, often only its text, and decides what to do next. Behind the package, the reader is your code, and you.

Points instead of rectangles

For a reader who writes code, the package offers things the protocol does not. The clearest is points. Where the server’s tools deal in space, they deal in rectangles. Georeferencing from coordinates returns a box and six affine coefficients, and the knowledge query takes a box. Field geology, however, is full of points: a sample site marked on the scan, a drill collar, a vertex on a fault trace.

The GeoReference object that georeference_bounds returns keeps its pyproj transformer. As a result, ref.pixel_to_lonlat(x, y) turns any pixel into longitude and latitude, and ref.lonlat_to_pixel goes back. No tool on the server does this. An agent holding only the coefficients would have to apply the affine and then reproject out of UTM by itself. That is exactly the work the module’s docstring calls “the deterministic core a VLM cannot do reliably.”

Fresh data, your own providers, no bundle files

Other advantages are quieter. Provider results land in a disk cache with no expiry. A live source’s version string never changes either. So the same live query, asked again, returns the first answer, however old it is. Only Python can switch that cache off: write_cache is a field on KnowledgeConfig, not an environment variable the server reads.

KnowledgeService(providers=[...]) accepts any object with an id, a supports test and a query method. So a state survey’s web service or a company’s collar table can produce the same bundle format. That list replaces the built-in providers, though, rather than adding to them.

A loop over three hundred sheets also writes no bundle files. The server writes a new bundle file for every knowledge query and has no cleanup tool. With write_cache=False, the loop writes nothing at all.

What a Geologic AI Script Loses Without the Server

What leaves is harder to see, because it is not a feature. It is a set of lessons.

Read the MCP layer’s comments and they turn into a field log of agents getting things wrong. One, from the second round of our agent tests, reads “Round-2 agents all called this tool first and then ignored it”. Their host showed them only a provider count from the capability check, and nothing about what was ready. Another notes that counting result items “hides the yield”, reporting “1 item” for 86 mineral occurrences. A third explains why an unread legend label comes back as null rather than an empty string. Both MCP agents in that round of testing had filled the blank with invented lithology names.

Lessons the adapter turned into code

Each lesson became code. The server’s text summary repeats every warning, so a text-only host cannot reduce the warnings to a count. It adds up the records each provider found and reports total_records_found, total_records_returned and truncated. When a query runs only the default providers, the server warns about each compatible provider it did not consult. That warning names the include value that would consult it. When a provider lacks its data, the warning names the step that installs it.

What the Python SDK bundle leaves out

The package returns the bundle before any of that happens. bundle.items holds one item per provider. So a query that finds 86 mineral occurrences comes back as a single item, and len(bundle.items) is 1. bundle.warnings holds what the providers said, but the note about skipped providers is not among them. It sits in bundle.trace["providers_not_consulted"], where nothing prints it.

Call query_bounds(box) without include and the mineral, land-cover and population providers are never asked. Your log says nothing about it. Likewise, a fit from two control points carries a residual_not_diagnostic warning in ref.warnings, a tuple your loop may never read.

Reader Mistakes, Not Model Mistakes

It is tempting to file these as agent mistakes, the kind a careful developer would not make. They are reader mistakes, and a loop written late on a Friday makes them too. It counts items instead of records, skips the warning list and takes an empty result as absence. And it does so three hundred times without hesitating. A model at least looks at each result. A loop looks only at what you told it to.

We did not write the adapter because models are careless. We wrote it because results are easy to misread, and we watched it happen in our own agent tests often enough to encode the fix. Those fixes are not specific to language models. Any reader that sees counts without context will draw the same wrong conclusion. So will any reader that sees a list without the note saying what is missing from it.

Who Decides the Next Step in a Geologic AI Workflow?

So the useful question is not “MCP or Python?” It is who decides the next step. If a model chooses which tool to call, give it the server: the adapter is its reading glasses. If your code fixes the steps, use the package and accept the trade. You get points, fresh live data, your own providers and loops that write no bundle files. You also become the adapter, and the few lines it would have written are now yours:

• Sum the record counts and check truncated.
• Print both warning lists.
• Read the trace for what was skipped.
• Accept no empty answer until you know which providers were asked.

None of this calls for new tooling. Both routes ship in the same open-source package, and the line between them is short enough to read over a coffee. The server is a thin layer, and the lessons it carries are a few dozen lines of ordinary Python. That makes them easy to borrow. A script that calls the package directly can take the same counting, warning list and check for skipped providers. Someone just has to write them down.

The package also installs without the server, its detector stack or an MCP host. So the SDK route runs on a laptop with one extra and one data download for the local earthquake and fault data. Start with the method that does the counting, because it shows exactly what a direct call leaves for you to do.

The Lines the Adapter Writes for You

The counting that a geologic AI script loses when it skips the server lives in one method of the MCP adapter, _bundle_result. That method runs after KnowledgeService.query has already returned its bundle (adapter.py:502, 509):

        # Each item is a per-provider *summary*; the records it found live in
        # `record_count` (found) and `value` (returned). Summarizing by item
        # count hides the yield -- e.g. "1 item" for 86 mineral occurrences --
        # so surface record totals plus a per-provider breakdown an agent can
        # branch on without re-summing `item.value` itself.
        record_counts: dict[str, int] = {}
        returned_counts: dict[str, int] = {}
        truncated = False
        for item in structured.get("items", []):
            provider = item.get("provider", "unknown")
            returned = _returned_count(item.get("value"))
            found = item.get("record_count")
            found = int(found) if found is not None else returned
            record_counts[provider] = record_counts.get(provider, 0) + found
            returned_counts[provider] = returned_counts.get(provider, 0) + returned
            truncated = truncated or bool(item.get("truncated")) or found > returned

        total_found = sum(record_counts.values())
        total_returned = sum(returned_counts.values())
        structured["record_counts"] = record_counts
        structured["total_records_found"] = total_found
        structured["total_records_returned"] = total_returned
        structured["truncated"] = truncated

        summary = (
            f"Knowledge query found {total_found} record(s) "
            f"across {len(record_counts)} provider(s)"
        )
        if total_returned < total_found:
            summary += f" ({total_returned} returned, truncated)"
        nonempty = {provider: count for provider, count in record_counts.items() if count}
        if nonempty:
            breakdown = ", ".join(f"{provider}={count}" for provider, count in nonempty.items())
            summary += f": {breakdown}."
        else:
            summary += "; no provider returned records."
        warnings = list(structured.get("warnings", []))
        if warnings:
            summary += f" {len(warnings)} warning(s): " + " | ".join(warnings)

For each provider item, the method compares record_count (what the provider found) with the length of value (what it returned). Next, it sums both per provider and marks the result truncated if any provider cut its list. Then it writes the line a text-only host shows the model: totals, a per-provider breakdown and every warning, joined into one string. Four lines earlier (line 664), _knowledge_notice_warnings (adapter.py:1067-1111) has already turned the trace’s providers_not_consulted list into warnings. Each one names a skipped provider and the include value that would consult it.

None of this is on the object the SDK hands you. query_bounds returns the same KnowledgeBundle the adapter starts from (knowledge/service.py:174-190), with bounds, items, selected_item_ids, warnings, provider_versions, trace and schema_version, and nothing more (knowledge/types.py:96-104). That is the trade in code: the engine is shared, the reading is not.

Running geologic AI from the Python SDK

To run the SDK route, clone the Stratigraphic Amenity repository and install only the georeferencing extra. Then fetch the knowledge base, which holds the local earthquake catalogue and the GEM fault file:

git clone https://github.com/EigenformAI/stratigraphic-amenity.git
cd stratigraphic-amenity
uv sync --extra geo
uv run --no-sync stratigraphic-amenity-assets peace-knowledge-base

Then write the reading yourself. Save it in the checkout as sdk_route.py and run it with uv run --no-sync python sdk_route.py:

from stratigraphic_amenity.georef import GroundControlPoint, georeference_bounds
from stratigraphic_amenity.knowledge import KnowledgeService

ref = georeference_bounds(
    crs="EPSG:26915",
    gcps=[
        GroundControlPoint(0, 0, 660000, 5400000),
        GroundControlPoint(1000, 1000, 690000, 5370000),
    ],
    pixel_extent=(0, 0, 1000, 1000),
)
print(ref.warnings)                        # residual_not_diagnostic: two points
lon, lat = ref.pixel_to_lonlat(412, 655)   # a site marked on the scan

bundle = KnowledgeService.from_env().query_bounds(ref.bounds)
found = {item.provider: item.record_count for item in bundle.items}
truncated = any(item.truncated for item in bundle.items)
skipped = [entry["provider"] for entry in bundle.trace["providers_not_consulted"]]
for warning in bundle.warnings:
    print("warning:", warning)

ref carries the resolved CRS, EPSG:4326 bounds and fit_method set to axis-aligned-exact-2gcp. It also carries a residual_not_diagnostic warning, because two points leave nothing to check against. The bundle holds one item each for earthquake_history and active_faults, the only default providers that accept a box. Each value holds up to 50 records, and record_count holds the full total. Meanwhile, skipped lists mineral_occurrences, landcover_distribution and population_density. The warnings also carry a fallback notice from each provider, such as “earthquake_history: using legacy local asset because no source mirror was found.”, for as long as no source mirror is synced.

What to Watch For When Geologic AI Runs in Python

The server’s checks become your checks. These are the traps a geologic AI script meets first:

• Count records, not items. len(bundle.items) counts the providers that answered. Add up record_count instead, and read truncated as a sign that a list stopped at its limit. That limit is 50 records per provider unless max_records or the config’s max_records_per_provider sets another.
• Read the trace for skipped providers. On the server, each provider skipped by default comes with a warning and the include value to retry with. In Python it is only an entry in bundle.trace["providers_not_consulted"], with the reason disabled_by_default.
• Missing data looks like a short answer. A default provider whose files are missing adds only “earthquake_history: configured provider asset or dependency is unavailable.” to the warnings, and its item is simply absent. Run stratigraphic-amenity-assets peace-knowledge-base to install the files.
• Two or three control points prove nothing. The fit adds “residual_not_diagnostic: an exact fit has no independent error check; supply at least 4 well-distributed GCPs”. Also, holdout_error stays None below four points.
• Empty is not absent. A fault search that finds nothing returns “No GEM active-fault features intersected the bounds; this is not evidence that no active faults exist.” A script drops that sentence easily, so report it next to the zero: no result is not evidence of absence.
• Bundles keep local paths. The server redacts paths; the SDK does not. An earthquake item’s source is the path of the local catalogue file. Clean a bundle before you share its JSON.

Live sources, caching and the antimeridian

Three more traps appear once a script reaches beyond the local files:

• Live sources need the network extra. The shipped mineral sources, the Ontario and Quebec examples, are queried over HTTP, which needs uv sync --extra knowledge-network. Without it, a box that overlaps their coverage fails with an error that names that command. If mineral_occurrences is the only provider in include, query_bounds raises the error. Next to other providers, it becomes one more warning.
• Live answers stay cached. A live result is reused until its cache file is deleted. For fresh data, build the service as KnowledgeService(config=replace(KnowledgeConfig.from_env(), write_cache=False)), with replace from dataclasses. It then neither reads nor writes the provider cache.
• Boxes across 180° need Python, with one exception. Bounds refuses a box whose minimum longitude exceeds its maximum, so a four-coordinate query through the server cannot cross the antimeridian. In Python, query_extent splits the extent in two for the earthquake, fault and mineral providers. The Earth Engine land-cover and population providers have no split-extent method, so they fail with “Provider ’landcover_distribution’ does not support antimeridian-split extent queries.” or its population twin.

Adding your own regions or countries

The quickest route is Python-only. KnowledgeService(providers=[...]) takes any object with an id, a supports(request) test and a query(request) method that returns KnowledgeItem objects. A provider for your state survey or your own collar table therefore needs no change to the repository. The catch is that the list replaces the built-in providers, so the service runs only what you pass.

The shipped regional sources are examples: Ontario and Quebec for mineral occurrences, DISS for Italy for active faults. Adding a region to those families, so that the server sees it too, is a code change. You write a source adapter (knowledge/sources/ogs_minerals.py is the Ontario model) and register it with a coverage_bounds rectangle in knowledge/sources/registry.py. Then you bind it in knowledge/service.py, where the Ontario, Quebec and DISS bindings live. Two settings decide whether it runs:

• GEOMAP_MINERAL_OCCURRENCE_SOURCE_IDS replaces the default list rather than extending it. Name ontario_mineral_deposit_inventory and sigeom_mineral_occurrences in it too if you want to keep them.
• A fault source added through GEOMAP_ACTIVE_FAULT_SOURCE_IDS is not queried by default. The fault provider asks only the first configured source unless provider_options selects others with source or sources.

FAQs

Can I run geologic AI tools without an MCP server?

Yes. Stratigraphic Amenity is a Python package as well as a server. Install the geo extra, import georeference_bounds and KnowledgeService, and call them from any script. The MCP server is a thin adapter over the same functions, so the fit, the providers and their warnings are identical. What you leave out are its summaries, record totals and stored bundles.

How do I convert a map pixel to latitude and longitude in Python?

Fit the map with georeference_bounds(crs=..., gcps=..., pixel_extent=...), then call pixel_to_lonlat(x, y) on the result. It applies the affine and reprojects to EPSG:4326 with the transformer the result kept. lonlat_to_pixel goes the other way. No MCP tool converts a single point, so this step belongs to the Python SDK.

Why does my geologic AI query skip mineral occurrences?

The mineral provider is off by default. Without include, query_bounds runs only the earthquake and fault providers on a box, and the rest appear in bundle.trace["providers_not_consulted"], not in the warnings. To ask it, pass include=("mineral_occurrences",) and install the knowledge-network extra, because its shipped Ontario and Quebec sources are queried over HTTP.

Can I add my own region or country?

Yes. In a script, pass your own provider object to KnowledgeService(providers=[...]), which replaces the built-in list without touching the repository. For a region the server sees too, extend the shipped examples (Ontario and Quebec minerals, DISS for Italy) in code: a source adapter, a coverage rectangle and a binding. Then list its id with the default ids, and for faults select it in provider_options.

Try It With Geocluster

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

When a model should choose the steps, the full map pipeline runs the same functions through the server, one tool call at a time.