Cerne Abbas Walk is an artwork by Richard Long, in the collection of the Tate Modern in London and on display at the time of this writing. Several of Long’s works involve geographic representations of his walks, some abstract and some concrete. Cerne Abbas Walk is described by the artist as “a six-day walk over all roads, lanes and double tracks inside a six-mile-wide circle centred on the Giant of Cerne Abbas.” The Tate catalog notes that “the map shows his route, retracing and re-crossing many roads to stay within a predetermined circle.”

The Giant in question is a 180-foot-high chalk figure carved into a hill near the village of Cerne Abbas in South West England. Some archaeologists believe it to be of Iron Age pedigree, some think it to date from the Roman or subsequent Saxon periods and yet others find the bulk of evidence to indicate a 17th-century origin as a political satire. (I find the last theory to be both the most amusing and the most convincing.)

I found the geographic premise of Cerne Abbas Walk intriguing, so I decided to replicate it computationally.

Preparing the Image and Creating the Map

In order to create map imagery centered at the correct location and at the correct zoom level, we must define a few initial variables.

Here, we assign a cropped version of the previously shown map image to a variable:

Engage with the code in this post by downloading the Wolfram Notebook

&#10005 longMapBG = Import["https://wolfr.am/CerneAbbasWalk-OriginalMap"];

To create the map, we then set the center location of our walk (the Cerne Abbas Giant), the radius of the disk bounding the walk we’re going to make (our region of interest) and the radius of the maps to generate along the way (the cropping region for longMapBG was predetermined to match this radius):

&#10005 centerLocation = GeoPosition[{50.813611`, -2.474722`}]; diskRadius = Quantity[6, "Miles"]/2; mapDiskRadius = Quantity[5, "Miles"]; disk = GeoDisk[centerLocation, diskRadius]; mapRange = GeoDisk[centerLocation, mapDiskRadius];

View satellite imagery of our center location to confirm that it is correctly located on the Cerne Abbas Giant:

&#10005 ImageRotate[GeoImage[GeoCenter -> centerLocation, GeoZoomLevel -> 16], 90 \[Degree]]

Superimpose the region of interest (ROI) we’ve created atop the original Cerne Abbas Walk image to confirm that it matches the circle in the original:

&#10005 Rasterize@ GeoGraphics[{Green, Opacity[.5], disk}, GeoCenter -> centerLocation, GeoRange -> mapRange, GeoBackground -> longMapBG]

Getting Road Data

Before we can perform any computations, we have to acquire data on roads in the vicinity of our ROI. I’ve chosen to use OpenStreetMap (OSM), a user-editable collaborative mapping project, as my data source, primarily because it has a powerful and open API for querying data, and there exists an OSMImport function in the Wolfram Function Repository that makes access to that API from within the Wolfram Language a breeze. Also, OSM has excellent semantic tagging of different types of map elements, which simplifies the task of filtering out just the roads.

Call the Wolfram Function Repository OSMImport function on our ROI, which will query the OSM API for all of the map polylines (“ways” in OSM parlance) that intersect the bounding box of our ROI:

&#10005 osmData = ResourceFunction["OSMImport"][disk];

Extract the “nodes” (polyline vertices) from the OSM response and parse their coordinates to their GeoPosition:

&#10005 nodes = Lookup["Position"] /@ osmData["Nodes"]

Plot the nodes:

&#10005 Rasterize@ GeoGraphics[{Red, PointSize[Small], Point /@ nodes}, GeoCenter -> centerLocation, GeoRange -> mapRange]

Extract the “ways” (polylines of nodes) from the OSM response, filter out those that aren’t marked as “highways” (roads) and replace the node IDs within each with each node’s GeoPosition:

&#10005 ways = Lookup["Nodes"] /* (Lookup[nodes, #] &) /@ Select[osmData["Ways"], KeyExistsQ[#Tags, "highway"] &]

Plot the roads as lines on a map, overlaying our ROI:

&#10005 Rasterize@ GeoGraphics[{Red, Line /@ ways, Green, Opacity[.5], disk}, GeoCenter -> centerLocation, GeoRange -> mapRange, ImageSize -> Full]

Trimming Road Data

The OpenStreetMap API takes only a rectangular bounding box as a region specification, so we have to crop it to our circular region of interest.

Here we create a 2D geometric region from our ROI, converting it to an ellipse proportional to latitude and longitude:

&#10005 diskRegion = Region@Disk[centerLocation[[1]], Abs[Subtract @@@ GeoBounds[disk]]/2]

Then we convert each road to a 1D geometric region (embedded in 2D space):

&#10005 wayRegions = Region[Line[#[[All, 1]]]] & /@ ways;

Trim each road by intersecting it with our ROI ellipse, discarding any empty regions (that is, roads that are entirely outside the ROI—the RegionDimension of an EmptyRegion is –∞):

&#10005 wayRegionsTrimmed = Select[RegionIntersection[#, diskRegion] & /@ wayRegions, RegionDimension[#] > 0 &];

Convert the trimmed road regions back to lists of GeoPosition objects:

&#10005 waysTrimmed = Map[GeoPosition, DiscretizeGraphics /* MeshCoordinates /@ wayRegionsTrimmed, {2}];

Plot the trimmed roads as lines on a map, overlaying our ROI:

&#10005 GeoGraphics[{Red, Line /@ waysTrimmed, Green, Opacity[.5], disk}, GeoCenter -> centerLocation, GeoRange -> mapRange]

Plot the trimmed roads as lines on a satellite image:

&#10005 GeoGraphics[{Red, Line /@ waysTrimmed}, GeoCenter -> centerLocation, GeoRange -> mapRange, GeoServer -> "DigitalGlobe"]

Plot the trimmed roads as lines on the original Cerne Abbas Walk image, overlaying our ROI (black lines are Long’s walk on the original image; red lines are the trimmed roads we created):

&#10005 Rasterize@ GeoGraphics[{Red, Thick, Line /@ waysTrimmed, Green, Opacity[.5], disk}, GeoCenter -> centerLocation, GeoRange -> mapRange, GeoBackground -> longMapBG, ImageSize -> Full]

Our trimmed roads match up pretty well to the black lines on the original piece, although there are some inconsistencies. These can probably be attributed to changes in the road layout over the past 40+ years, as well as to missing data for some walking paths on OpenStreetMap—in rural England, the distinction between “public footway” and “dirt path on someone’s farm” can become blurred.

Making a Graph

By modeling the network of interconnected roads as a graph, we can use the graph theoretic functionality of the Wolfram Language to perform computations on the network.

Each road is described as a polyline consisting of multiple vertices, so we convert each pair of consecutive vertices to a graph edge:

&#10005 roadEdges = Flatten[UndirectedEdge @@@ Partition[#, 2, 1] & /@ Values[waysTrimmed]];

Combine the edges into a graph (which may contain disconnected components representing roads that cannot be reached from the main component):

&#10005 roadsDisconnectedGraph = Graph[roadEdges]

Since our final walk cannot traverse roads unreachable from the main network, we must remove such “orphaned” components. We do so by extracting the largest connected component (the main road network) and, separately, all smaller orphaned components from the graph:

&#10005 connectedComponents = ConnectedGraphComponents[roadsDisconnectedGraph]; Column[{roadsGraph, roadsOrphanedGraphs} = Through[{First, Rest}[connectedComponents]]]

See which vertices are orphaned:

&#10005 GeoGraphics[{Red, Thickness[.008], VertexList /* Line /@ roadsOrphanedGraphs, Green, Opacity[.5], disk}, GeoCenter -> centerLocation]

These two components happen to represent, respectively, a dead-end road jutting in from just outside the ROI and an isolated walkway in an Anglican friary—nothing of consequence.

We can find a more optimal tour by weighting each edge of the graph by the geographical distance between its two vertices. Here, we use EdgeWeight to assign a numerical weight to each edge, and the "SpringElectricalEmbedding" GraphLayout, which takes weights into account when laying out the graph visually:

&#10005 weightedRoadsGraph = Graph[roadsGraph, EdgeWeight -> QuantityMagnitude[ EdgeList[roadsGraph] /. UndirectedEdge -> GeoDistance, "Feet"], GraphLayout -> {"SpringElectricalEmbedding", "EdgeWeighted" -> True}]

Finding a Tour

Now that we’ve created a graph representing the road network within our ROI, we can find a walk that traverses every road at least once. In graph theory, this is the so-called Chinese postman problem—finding the shortest path (“tour”) that visits each edge of a graph at least once. The FindPostmanTour function automatically takes edge weights into account, which results in a postman tour minimizing real geographical distance.

Find a Chinese postman tour of the graph, get the first vertex of each edge in the result and append the second vertex of the last edge to close the path, yielding the sequence of vertices to visit:

&#10005 edgeTour = FindPostmanTour[weightedRoadsGraph][[1]]; roadsTour = Append[edgeTour[[All, 1]], edgeTour[[-1, 2]]]

Plot the tour as a single line on a map, overlaying our ROI:

&#10005 GeoGraphics[{Red, Line[roadsTour], Green, Opacity[.5], disk}, GeoCenter -> centerLocation, GeoRange -> mapRange]

Find the total distance of the tour:

&#10005 tourDistance = GeoDistance[roadsTour]

Calculate the time it would take to walk the tour, assuming a walking speed of five kilometers per hour:

&#10005 walkingSpeed = Quantity[5, ("Kilometers")/("Hours")]; walkingDuration = UnitConvert[tourDistance/walkingSpeed, MixedUnit[{"Days", "Hours", "Minutes", "Seconds"}]]

(The OpenStreetMap database is constantly being updated with new data, so don’t worry if you get different figures from these calculations—they’re probably more accurate than mine!)

This is well within Richard Long’s six-day figure, allowing plenty of time for rest, eating and general artistic contemplation. (Of course, Long didn’t have the Wolfram Language back in 1975 to help him plan an optimal route!)

Making an Animation

Create an animation playing at 10,000x walking speed that shows the tour at each step, along with the time elapsed:

&#10005 tourDistanceList = GeoDistanceList[roadsTour]; animationSpeedMultiplier = 10000; animation = Animate[GeoGraphics[ {Red, Line[roadsTour[[;; t + 1]]], Green, Opacity[.5]}, ImageSize -> 500, GeoCenter -> centerLocation, GeoRange -> mapRange, GeoServer -> "DigitalGlobe", Epilog -> Inset[Framed[Style[ ToString@ QuantityForm[ UnitConvert[ Floor[Total[tourDistanceList[[;; t]]]/walkingSpeed, Quantity[1, "Minutes"]], MixedUnit[{"Days", "Hours", "Minutes"}]], "LongForm"], Orange, 20 ], Sequence[ FrameStyle -> GrayLevel[0.9], Background -> GrayLevel[1], RoundingRadius -> 5] ], Scaled[{.025, .025}], {Left, Bottom}] ], {{t, 1, "Step"}, 1, Length[tourDistanceList], 1}, Sequence[ AnimationRunning -> False, AnimationRepetitions -> 1, SaveDefinitions -> True, DefaultDuration -> QuantityMagnitude[ animationSpeedMultiplier^(-1) walkingDuration, "Seconds"]] ]

Following the Road Further

This code attempts to reproduce the path of Long’s original walk as closely as possible, but there are many possible permutations on the idea that beg further exploration. What about walks within a polygonal region? A walk within a region shaped like the Giant it’s centered on? Walks centered on other landmarks? How does the length of the walk vary as a function of the enclosing circle’s diameter? How does the diameter of the circle needed to enclose a walk of a fixed distance vary as one moves from urban to rural areas? If you pursue any of these questions, or perform your own explorations based on this or other works of art, share your results on Wolfram Community for myself and others to see!