Browse by Topic

# Taking the Cerne Abbas Walk: From Conceptual Art to Computational Art

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:

 ✕ `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):

 ✕ ```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:

 ✕ ```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:

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

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:

 ✕ `osmData = ResourceFunction["OSMImport"][disk];`

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

 ✕ `nodes = Lookup["Position"] /@ osmData["Nodes"]`

Plot the nodes:

 ✕ ```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:

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

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

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

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:

 ✕ ```diskRegion = Region@Disk[centerLocation[[1]], Abs[Subtract @@@ GeoBounds[disk]]/2]```

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

 ✕ `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 –):

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

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

 ✕ ```waysTrimmed = Map[GeoPosition, DiscretizeGraphics /* MeshCoordinates /@ wayRegionsTrimmed, {2}];```

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

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

Plot the trimmed roads as lines on a satellite image:

 ✕ ```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):

 ✕ ```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:

 ✕ ```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):

 ✕ `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:

 ✕ ```connectedComponents = ConnectedGraphComponents[roadsDisconnectedGraph]; Column[{roadsGraph, roadsOrphanedGraphs} = Through[{First, Rest}[connectedComponents]]]```

See which vertices are orphaned:

 ✕ ```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:

 ✕ ```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:

 ✕ ```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:

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

Find the total distance of the tour:

 ✕ `tourDistance = GeoDistance[roadsTour]`

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

 ✕ ```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:

 ✕ ```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"]] ]```