Monday, 9 November 2026

Exploring Archaeological Landscapes, Part II (long)

Fourth in a GIS 💖 AI mini-series. I was interested in possible links in Béarn, SW France where I am now, to Ley Lines I learned about when in England. Alfred Watkins reinterpreted in 1921-1925 possible alignments of monuments, tumuli & churches (BBC). Part I yielded only six recorded archaeological points: the Béarn exercise is a landscape case study rather than a statistical test of alignments; GIS can reveal shared relationships with relief, water, routes, and geology, but cannot demonstrate intentional alignment or validate ley-line claims. Here we turn to Cambridgeshire (East Anglia, where I was before) that has over 100 monuments & other data to expand this project. 

Note: as in all my AI-assisted posts and AI-guided work, my text is in italics.


Quick take

There are no apparent strong lineaments or Ley Lines to map in central East Anglia from the evidence collected. Would casting a wider net from Béarn thru Cambridgeshire & Peterborough to country-wide help? That would tax the Google Colab (scripting) & Open Street Map or BGS GeoIndex Onshore (data sources) internet use & download; these are kept free under "fair use" for personal research. It took for example over 5 min. to retrieve OSM data for Cambridgeshire & Peterborough alone, and over-traffic errors started during EU afternoons... when US started to get online!

Here we missed the whole Monte Carlo simulation of linear matches that may/not produce a Digital Twin for Ley Lines: the idea is to find alignments that can be modeled to create a "digital twin" - essentially a line (or series of lines) that exists only on paper. However, if these were to coincide with "ley lines" documented through other means, it would provide a form of scientific support. There is no proof, of course - correlation is not causation - but it certainly keeps us cartographers quite busy along with AI ​​assistants!

This is all under Creative Commons (CC BY-SA 4.O): the prompt thread beginning Part I is here for you to re-use; warning, there are not a few twists & turns, so drop a comment below if specifics interest you.


Cambridgeshire GIS Workflow

This workflow extracts OpenStreetMap (OSM) data for Cambridgeshire using a county-wide bounding box, then summarizes OSM line features within 1 km buffers around representative points.

Important: The “2 points” are not OSM results. They are representative points derived from the two polygons in the ONS “Cambridgeshire & Peterborough” layer. The buffers and OSM summary are therefore based on those two points.

1. Define the study area in QGIS

Where: QGIS Python Console.
Input: ONS Cambridgeshire & Peterborough polygon layer.
Outputs:
  • cambs_sites — two representative points
  • cambs_buffers_1km — two 1 km buffers
  • A county-wide bounding box in latitude/longitude

The script creates one point inside each polygon, buffers each point by 1 km, and then calculates a bbox from the full extent of the ONS layer, not from the buffers.

The bbox is a rectangle covering the whole county area, not just a corridor.

2. Extract OSM data in Colab

Where: Google Colab.
Input: The county-wide bbox from QGIS.
Output: A GeoPackage containing OSM roads, paths, and waterways.

The Colab script downloads OSM line features within the county-wide bbox, clips them to that rectangle, separates roads from paths, and exports the result as a GeoPackage.

It creates three layers:

  • osm_roads — roads and similar highway lines
  • osm_paths — footways, tracks, cycleways, etc.
  • osm_water — rivers, streams, canals, drains, ditches

3. Download the GeoPackage

Where: Google Colab.
Output: A .gpkg file saved to your computer.

After the extraction finishes, run the download cell. This saves the OSM GeoPackage from Colab to your local Downloads folder.

4. Load OSM data into QGIS

Where: QGIS.
Input: The downloaded GeoPackage.
Output: OSM layers available for analysis.

Add the GeoPackage to QGIS and rename the roads layer to osm_roads if needed. Script 3 expects osm_roads, osm_paths, and osm_water.

5. Run the analysis in QGIS

Where: QGIS Python Console.
Input: OSM layers plus cambs_buffers_1km.
Output: A summary table with OSM statistics per buffer.

Script 3 takes each 1 km buffer and calculates:

  • How many OSM road, path, and waterway lines intersect it.
  • Total length of those lines inside the buffer, in metres.

It produces two rows—one per buffer—with fields such as:

  • road_count
  • road_length_m
  • path_count
  • path_length_m
  • water_count
  • water_length_m

Key distinction

Item What it is Where it comes from
cambs_sites Two representative points Derived from the two polygons in the ONS layer
cambs_buffers_1km Two 1 km buffers Created around the two representative points
osm_roads OSM road lines Downloaded from OSM within the county-wide bbox
osm_paths OSM path lines Downloaded from OSM within the county-wide bbox
osm_water OSM waterway lines Downloaded from OSM within the county-wide bbox
cambs_summary_*.gpkg Buffer-level OSM summary Created by script 3

Maps & table

Below are two maps, one for the current 2 centroids, and one for the next section for 108 monuments. 

click to enlarge, original here

The table shows the proximity analysis for 2 centroids, the one for 108 monuments is posted here.

click to enlarge, original posted here

Why the first attempt looked sparse

The intended study area covered Cambridgeshire and Peterborough, together forming the geographical footprint of the Cambridgeshire & Peterborough Combined Authority. Established in 2017, the Combined Authority brings together the area's councils rather than replacing them with a single merged council.

In the first attempt, however, the workflow used two representative centre points instead of the full administrative polygons. The initial bounding box was calculated from the combined extent of two 1 km buffers around those points. With the points arranged roughly north–south, this produced a narrow corridor rather than a rectangle covering the whole study area.

That explains why the OpenStreetMap extract initially looked sparse: the search covered only that restricted corridor. The apparent lack of features was a consequence of the search extent, not an OpenStreetMap retrieval error.

Rerunning the workflow with the full multi-polygon geometry produced the intended county-wide search extent. The two attempts therefore illustrate an important distinction: buffering representative points is not the same as using the areas those points represent.

Current setup

The updated script 1 still creates the two site points and buffers, but now calculates the bbox from the full extent of the ONS layer. This means the OSM extract covers a much larger area, while script 3 still summarizes OSM within the two 1 km buffers.


Monuments, routes and radon: mapping proximity without mixing the signals

This stage of the project brought together monument locations, mapped roads and paths, and the BGS–UKHSA Radon Indicative Atlas. The aim was to explore how monuments sit within their surrounding landscape while keeping three different questions separate:
  • How close is each monument to a mapped road or path?
  • Where do the recorded monuments cluster?
  • What radon-potential category is mapped in the surrounding grid square?

The central distinction is proximity. Measuring how near a monument is to a road, path or mapped area does not, by itself, demonstrate directional alignment, historical intent or a causal relationship.

Building the monument dataset in QGIS

The starting point was a monument point layer, together with OpenStreetMap road and path layers. The working data were reprojected to British National Grid (EPSG:27700), allowing the nearest-feature calculations to use meters.

Two successive nearest-feature joins attached information about the nearest mapped road and nearest mapped path to each monument. The resulting point layer retained the monument attributes alongside the joined route attributes and distance fields.

These distances describe the shortest mapped separation between a monument and the relevant line feature. They are not walking distances, measures of public access, or evidence that the modern route existed when the monument was placed. Their interpretation also depends on the completeness and classification of the source mapping.

Add geologic faults & aeromagnetic data?

The next obvious sources of earth disruptors that may have an effect are faults, & aeromagnetic discontinuities that map faults underground - see my thesis map revisit starting here - but BGS fault map shows nothing at the surface, and aeromagnetic data aren't for free from the same source: they would also entail a whole lot of interpretation and software or professional community I no longer have access to; even my thesis revisit above is qualitative at a continental scale.  

click to enlarge, original here

Adding radon context

The next step added the downloaded Radon Indicative Atlas as a polygon layer. Our area-of-interest extract contained 6,450 grid tiles, with values spanning classes 1 to 6 and most tiles in class 1.

Rather than calculating distance to these tiles, the workflow identified which tile intersected each monument point and attached its radon attributes. This is a containing-area comparison, not another nearest-line calculation.

The important source field is CLASS_MAX. It records the highest mapped radon-potential class occurring anywhere within a 1 km grid square. It does not necessarily describe the precise ground beneath the monument.

click to enlarge, original here

These bands concern the estimated percentage of homes exceeding the indoor radon Action Level. They are not outdoor gas measurements, percentages of radioactive ground, or direct measurements at the monument locations. The atlas supplies regional context; it is not a fault map or an independent measurement of ground permeability.

What the new columns record

  • radon_TILE: the reference of the selected grid square.
  • radon_CLASS_MAX: the maximum mapped radon-potential class within that square.
  • radon_VERSION: the version recorded in the source dataset.
  • radon_band: a readable percentage-band label.
  • radon_matches: the number of tiles intersecting the point.
  • radon_status: a flag showing whether the point had a single match, multiple matches, no match or another issue.

The updated script preserved one output row for every input alignment row. Where a point intersected more than one tile, it selected the highest valid class and flagged that choice. Unmatched points were retained rather than silently discarded.

A new GeoPackage stored the enriched point layer, leaving the earlier results untouched. Because it already contains the monument geometry and inherited alignment attributes, a further join back to the original monument layer is optional.

Keeping the maps interpretable

A heatmap weighted by radon class initially offered an eye-catching view, but it combined two different things: monument clustering and the numeric class assigned to each point. We therefore separated the displays.

1. Radon context with monument locations

The radon polygons are best shown with categorized colors using CLASS_MAX, following the six atlas bands. Ordinary monument points sit above them. This keeps the grid boundaries visible and makes clear that the radon information belongs to an area rather than to a measured point.

2. Monument density

A separate heatmap gives every monument equal weight. It answers the simpler question: where do the recorded monuments cluster?

The heatmap radius controls how far each point's influence spreads. It should be documented and kept consistent when comparing maps. Any duplicate monument records must also be checked before interpreting the density pattern.

3. Distance to roads and paths

For individual distances, graduated point symbols are clearer than a distance-weighted heatmap. Each monument can be colored or sized according to its own verified road-distance or path-distance field, with the relevant network visible underneath.

A distance-weighted heatmap would combine distance with clustering: several moderately distant monuments could create a stronger hotspot than one very remote monument. Graduated points avoid that ambiguity.

Here, “isolation” must always be qualified. Distance from a mapped road is not the same as distance from a mapped path, and neither measures isolation from other monuments.

Presenting the comparison

The radon-context map and the equally weighted monument-density map are placed side by side. Note:

  • Radon classes describe containing grid squares, not measurements at monument locations. 
  • Monument density does not demonstrate association with radon.

click to enlarge, original here


Checks before drawing conclusions

Matching input and output row counts is necessary, but it is not enough. Nearest-feature joins can return multiple equally close matches, so the number of unique monument identifiers should also be checked. The radon script preserves input rows; it does not automatically resolve duplicates inherited from earlier joins.

  • Confirm the working coordinate reference system and distance units.
  • Verify which distance field belongs to roads and which belongs to paths.
  • Check unique monument identifiers as well as total row counts.
  • Inspect unmatched, multiple-tile and invalid-class radon records.
  • Remove temporary filters before calculating full-dataset totals.
  • Use a persistent monument identifier for any later joins, rather than the GeoPackage row key or radon tile reference.

The next analytical question

Finding monuments in higher-class tiles is not enough to demonstrate an association. The relevant comparison is whether those tiles contain a larger share of monuments than their share of the study area would suggest.

The next step is therefore to summarize unique monuments by radon class, then compare their percentages with the percentage of the exact area of interest occupied by each class. Tiles crossing the study boundary should be clipped and compared by area rather than simply counted.

Even a difference between monument percentages and area percentages would remain exploratory. Settlement patterns, land use, survival, accessibility and recording effort could all influence the observed distribution. A formal test would need an appropriate comparison design.

For now, the result is a reusable spatial dataset and a clearer mapping method: show radon context as areas, monument clustering as density, and individual route proximity as distances. Keeping those signals separate makes the patterns easier to examine without turning an attractive map into an unsupported explanation.

Data and software

  • QGIS for reprojection, nearest-feature joins, spatial matching, symbology and print layouts.
  • OpenStreetMap source layers for the mapped monuments, roads and paths. © OpenStreetMap contributors.
  • BGS–UKHSA Radon Indicative Atlas for the radon grid and class definitions. The source dataset's license and attribution requirements should accompany any published maps.


Reading the wider landscape

Visually, the radon-class pattern appears to pick out parts of the northern and central lowland landscape extending towards the Wash. This is a qualitative observation, not a calculated relationship. The radon atlas provides geological context rather than a substitute for a bedrock or elevation map.

The flood-risk and sea-level-rise maps below add separate perspectives on the landscape - see here for more on this topic - they are presented here for descriptive comparison; no additional proximity calculations or statistical tests have been applied to these layers.

click to enlarge, original here

Each map should be read according to its own source, date, legend and assumptions. Flood-risk mapping and sea-level-rise scenarios should not be treated as interchangeable, or as reconstructions of the landscape when the monuments were established.

Concluding remark

As mentioned b re: faults & aeromagnetic data, the maps above most definitely show congruences at a regional scale. In the prompt thread, AI mentioned:

Call radon geological context, rather than a “bedrock substitute”. BGS derives radon potential from both bedrock and superficial geology, combined with indoor measurements. So the pattern you’re noticing may be useful, but it cannot independently establish bedrock boundaries or low-lying terrain.

So Radon definitely outlines the high ground running N-S to the E, as well as lesser one running SW-NE in the Center. This Space Shuttle view will give you an appreciation of the relief in East Anglia. Note that views are looking south and Cambridgeshire is always at top right.

click to enlarge, original here

This is the view southward from The Wash (North Sea embayment) toward the Chalk Hills south of Cambridge. High Norfolk at left is the high ground to the East mentioned above, and "Isle of Ely" is the lesser high ground in the Center also mentioned above: this is a cleaned up frame ~ 15 s. in my Sea Level Rise video here; AI allowed to remove the YouTube onscreen display.

L is map location, R is picture perspective above, The Wash is blue

Saturday, 7 November 2026

A little AI arty romp

And now for something completely different again... look Ma, only the pics are AI here!

AI and humor

Taking off of  Medium article "Does AI have a sense of humor?" here. This after a British Library exhibit on 100th anniversary of the Russian Revolution here & specifically Lenin's library card here.

This is the prompt Perplexity AI suggested when I asked it to illustrate "Draw Lenin's gleaming scalp in the diffuse library lanterns as he leans over poring over books, as a pretty librarian brings him a stack of books, while a v severe looking wife sizes them up in the distance at the doorway come to fetch him home":

Friday, 6 November 2026

Exploring Archaeological Landscapes with QGIS and Python

Third of what became a GIS 💖 AI mini-series starting here and continuing here. 

Update: #4 in this mini-series here showed a little more success... though no Digital Twin for Ley Lines!

And now for something completely different... This gets technical, so here's the synopsis: A local Bio-energetician uses GeoBiology (think outdoors Feng Shui) to diagnose clients with energy fields used for alt. medicine. She is also a Water Douser to find underground waterways & springs. That got me interested in possible links to Ley Lines I learned about in England, Alfred Watkins' 1921-1925 modern reinterpretation of possible alignments of monuments, tumuli & churches (BBC). 

AI quickly provided the areal extent of Béarn (where I am in southwestern-most France), the "Atlas des Patrimoines" for historic features  and OpenStreetMap for culture data to use in the following workflow. Reading lines in the Béarn landscape, for example, the GR 653 pilgrimage route and the Aspe valley reveal how relief, water, and historic paths shape the region’s enduring corridors (ffrandonnee).

Château de Pau and Sauveterre-de-Béarn

Monday, 2 November 2026

Karl Bushby: 28 years, 31,000 miles, and the long walk home

Update 2: second of what became a GIS 💖 AI mini-series starting here and continuing here.

Update 1: added a cool 3D globe view of same from a GIS Educator at bottom.

Following on map trajectory posts in the main blog here and this blog here, here's another AI-assisted one from this prompt. See Author's GIS note at bottom.


On 2 November 2026, British adventurer Karl Bushby completed the Goliath Expedition: a 28-year, roughly 31,000-mile (50,000 km) journey on foot from Punta Arenas in southern Chile to his mother’s home in Hull, England.

His rule was brutally simple: no transport could advance the journey. He could fly or sail only to wait out bureaucracy, weather or funding problems—but he had to return to the exact point where he had stopped.