30:T441,So instead I chose to overlay several density surfaces with different bandwidths to build up a multi-scale density surface, much like a sound wave is made up of oscillations at multiple frequencies. The intention was to be able to show national, regional and local trends in site concentration all in a single map. To do this, I couldn't use Plot's [density mark](https://observablehq.com/plot/marks/density) as it does not expose the density surface as a grid. Instead I used Jeff Heer's [fast KDE library](https://observablehq.com/@uwdata/fast-kde) to compute the raster grids of calculated density with different bandwidths. I used bandwidths of 1km, 2km, 4km, 8km, 16km, 32km, and 64km, accumulating them with weights of 1, 4, 9, 16, 25, 36 and 49 respectively. Increasing the weighting of the larger bandwidths meant that the very high local densities didn't dominate the cumulative surface. One further advantage of dealing with the density raster grids explicitly is they can be clipped to the coastline of Britain with a simple point-in-polygon test against the `land` GeoJSON.31:Tb91,function multiscaleKDE(pts, bounds, options = {}) { const { bandwidths = [10000] } = options; // Normalize points to the grid dimensions const normPts = pts.map((d) => [ ((d.x - bounds.x1) / (bounds.x2 - bounds.x1)) * bounds.width, ((d.y - bounds.y1) / (bounds.y2 - bounds.y1)) * bounds.height ]); // Ensure full geographic space is used for raster by including two dummy points normPts.push([0, 0]); normPts.push([bounds.width - 1, bounds.height - 1]); const values = new Array(bounds.width * bounds.height).fill(0); bandwidths.forEach((bandwidth, w) => { const scaledBwx = (bandwidth / (bounds.x2 - bounds.x1)) * bounds.width; const scaledBwy = (bandwidth / (bounds.y2 - bounds.y1)) * bounds.height; // Calculate density for current bandwidth const densityGrid = kde .density2d(normPts, { bandwidth: [scaledBwx, scaledBwy], bins: [bounds.width, bounds.height] }) .grid(); // Determine the scaling factor for this raster const scaleFactor = 1 + (3 * bandwidth) / (bounds.x2 - bounds.x1); const scaledWidth = Math.round(bounds.width * scaleFactor); const scaledHeight = Math.round(bounds.height * scaleFactor); // Scale the density grid const scaledGrid = scaleGrid( densityGrid, bounds.width, bounds.height, scaledWidth, scaledHeight ); // Determine the overlap region const xOffset = Math.floor((scaledWidth - bounds.width) / 2); const yOffset = Math.floor((scaledHeight - bounds.height) / 2); // Accumulate scaled values into the final grid for (let y = 0; y < bounds.height; y++) { for (let x = 0; x < bounds.width; x++) { const scaledX = x + xOffset; const scaledY = y + yOffset; if ( scaledX >= 0 && scaledX < scaledWidth && scaledY >= 0 && scaledY < scaledHeight ) { const scaledIndex = scaledY * scaledWidth + scaledX; const i = y * bounds.width + x; values[i] += Math.pow(w + 1, 2) * scaledGrid[scaledIndex]; } } } }); // Create a land-clipped version for the raster display. const clippedValues = []; for (let y = 0; y < bounds.height; y++) { for (let x = 0; x < bounds.width; x++) { // Convert raster cell coordinates to geographic coordinates const geoX = bounds.x1 + (x / bounds.width) * (bounds.x2 - bounds.x1); const geoY = bounds.y1 + (y / bounds.height) * (bounds.y2 - bounds.y1); // Check if the point is in the polygon const point = { type: "Point", coordinates: [geoX, geoY] }; const isInPolygon = T.pointInPolygon(point, land[0].geometry); // Set cell to null if it's not in the polygon const i = y * bounds.width + x; clippedValues.push(isInPolygon ? values[i] : null); } } return { values, clippedValues, width: bounds.width, height: bounds.height }; }32:T74a,function drawMap(ptData, options = {}) { const { title = "", scheme = "Viridis" } = options; const grid = multiscaleKDE(ptData, bounds, { bandwidths: [1000, 2000, 4000, 8000, 16000, 32000, 64000] }); return Plot.plot({ width: bounds.width * 1.5 + 40, height: bounds.height * 1.5 + 25, inset: 20, insetTop: 5, projection: { type: "reflect-y", domain: bounds.feature }, color: { scheme: scheme, reverse: true, type: "sqrt" }, marks: [ Plot.frame({ fill: "oklch(90% 0.03 180deg)" }), Plot.geo(land, { stroke: "oklch(95% 0.02 180deg)", strokeWidth: 18 }), Plot.raster(grid.clippedValues, { ...bounds }), Plot.contour(grid.values, { ...bounds, stroke: "oklch(20% 0.03 180deg)", strokeWidth: 0.5 }), Plot.geo(land, { stroke: "black", strokeWidth: 2 }), Plot.dot(ptData, { x: "x", y: "y", r: 0.1, opacity: 0.1, stroke: "black", tip: true, title: (d) => `${Object.entries(d.properties).join("\n")}` }), Plot.text(["Concentration of\narchaeological sites"], { frameAnchor: "top-left", lineAnchor: "top", fontSize: bounds.width / 10, fontFamily: "Crimson Text", fontWeight: 700, dx: 8, dy: 8 }), Plot.text(["As recorded in Open Street Map"], { frameAnchor: "top-left", lineAnchor: "top", fontSize: bounds.width / 15, fontFamily: "Crimson Text", fontWeight: 400, dx: 8, dy: bounds.width / 4 }), Plot.text([title], { frameAnchor: "top-left", lineAnchor: "top", fontSize: bounds.width / 13, fontFamily: "Crimson Text", fontWeight: 700, dx: 8, dy: bounds.width / 2.5 }), Plot.frame() ] }); }33:[{"id":"da677f67dd96fdfe3992398a85effb2e9467cffc5c58fbbd74886e8789836939191bda4bdd8674c585bce4126869f1e7fae45117ba88ec1395f16ea0db858b47","href":"https://static.observableusercontent.com/files/da677f67dd96fdfe3992398a85effb2e9467cffc5c58fbbd74886e8789836939191bda4bdd8674c585bce4126869f1e7fae45117ba88ec1395f16ea0db858b47","name":"gbLand.json","type":"application/json","lastModified":1732052807847.449,"size":10917},{"id":"398442166ea88fc9ac7cd327801d9592d1dfcc98c96fceb47d686fcc60cc3ee7587a5f8f756c9c870490d92043cc39eadfc2d6b2620b016996ea07d6a6b65929","href":"https://static.observableusercontent.com/files/398442166ea88fc9ac7cd327801d9592d1dfcc98c96fceb47d686fcc60cc3ee7587a5f8f756c9c870490d92043cc39eadfc2d6b2620b016996ea07d6a6b65929","name":"archaeology.json","type":"application/json","lastModified":1732105385686.097,"size":1850230}] 2e:{"title":"30 Day Map Challenge: 20. Data: OpenStreetMap","cells":[{"id":0,"value":"# 30 Day Map Challenge: 20. Data: OpenStreetMap","pinned":false,"hidden":false,"mode":"md"},{"id":783,"value":"When faced with the topic 'OpenStreetMap', part of the challenge is choosing which parts of the huge database to work with. After some browsing of the data, I elected to look at British archaeology, as this reflects the depth of detail available and also the challenges of inconsistent coding that is inevitable with such a large crowdsourced project.","pinned":false,"hidden":false,"mode":"md"},{"id":69,"value":"## Extracting data from OSM","pinned":false,"hidden":false,"mode":"md"},{"id":785,"value":"One of the advantages of working with OSM is the ease with which data can be extracted. Using [Overpass turbo](https://overpass-turbo.eu), a simple script can extract all the data tagged with `historic:archaeological_site` as GeoJSON files:","pinned":false,"hidden":false,"mode":"md"},{"id":67,"value":"```text\n[out:json][timeout:25];\narea[\"name\"=\"Great Britain\"]->.searchArea;\n(\n node[\"historic\"=\"archaeological_site\"](area.searchArea);\n way[\"historic\"=\"archaeological_site\"](area.searchArea);\n relation[\"historic\"=\"archaeological_site\"](area.searchArea);\n);\nout body;\n>;\nout skel qt;\n```","pinned":false,"hidden":false,"mode":"md"},{"id":787,"value":"I then used [mapShaper](https://mapshaper.org) to convert the line and polygon features into points (centroids) and merged them to form a single layer of archaeological-related point locations and properties. I removed unnecessary properties such as `@id` and `url` to leave those that best described the feature. In total there were over 11,000 such features:","pinned":false,"hidden":false,"mode":"md"},{"id":441,"value":"archaeology = TDMC.geoToOSGB(\n topojson.feature(\n await d3.json(await FileAttachment(\"archaeology.json\").url()),\n \"archaeology\"\n ).features\n).map((d) => ({\n x: d.geometry.coordinates[0],\n y: d.geometry.coordinates[1],\n properties: d.properties\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":790,"value":"For mapping context I also used mapShaper to to create a generalised map of Britain into a single merged multi-polygon TopoJSON from which the bounds for each map could be stored:","pinned":false,"hidden":false,"mode":"md"},{"id":50,"value":"land = TDMC.geoToOSGB(\n topojson.feature(\n await d3.json(await FileAttachment(\"gbLand.json\").url()),\n \"land\"\n ).features\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"bounds = {\n const bs = TDMC.boundsFromGeo(\n { type: \"FeatureCollection\", features: land },\n { inset: 0 }\n );\n\n const w = 300; // Width of raster storing density surface\n return {\n ...bs,\n width: w,\n height: Math.round((w * (bs.y2 - bs.y1)) / (bs.x2 - bs.x1))\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":184,"value":"## Creating a multiscale clipped density surface","pinned":false,"hidden":false,"mode":"md"},{"id":816,"value":"Rather than just simply display the point locations of the archaeological sites","pinned":false,"hidden":false,"mode":"md"},{"id":819,"value":"Plot.plot({\n width: 280,\n projection: { type: \"reflect-y\", domain: bounds.feature },\n marks: [\n Plot.dot(archaeology, { x: \"x\", y: \"y\", r: 0.6, fill: \"firebrick\" }),\n Plot.geo(land)\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":834,"value":"I wanted to show the concentration of sites more explicitly especially as it is hard to distinguish areas of concentrated sites at this scale. An obvious way of doing this is to calculate a density surface using [Kernel Density Estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation). Observable Plot has a [density mark](https://observablehq.com/plot/marks/density) built in to make this easy:","pinned":false,"hidden":false,"mode":"md"},{"id":116,"value":"Plot.plot({\n width: 280,\n projection: { type: \"reflect-y\", domain: bounds.feature },\n marks: [\n Plot.density(archaeology, {\n x: \"x\",\n y: \"y\",\n bandwidth,\n fill: \"density\",\n opacity: 0.2,\n stroke: \"black\"\n }),\n Plot.geo(land)\n // Plot.dot(archaeology, { x: \"x\", y: \"y\", r: 2, fill: \"red\" })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":43,"value":"viewof bandwidth = Inputs.range([1, 30], {\n label: \"Bandwidth\",\n step: 1,\n value: 20,\n transform: Math.log\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":848,"value":"The problem with any single density surface is that it is highly dependent on the _bandwidth_: A density of points varies according to how wide an area is used to count points. Moving the slider above shows how density varies with scale.","pinned":false,"hidden":false,"mode":"md"},{"id":855,"value":"$30","pinned":false,"hidden":false,"mode":"md"},{"id":569,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":571,"value":"function scaleGrid(grid, originalWidth, originalHeight, newWidth, newHeight) {\n const scaledGrid = new Array(newWidth * newHeight).fill(0);\n\n for (let y = 0; y < newHeight; y++) {\n for (let x = 0; x < newWidth; x++) {\n // Map new coordinates back to original grid\n const origX = (x / newWidth) * originalWidth;\n const origY = (y / newHeight) * originalHeight;\n\n const x0 = Math.floor(origX);\n const x1 = Math.min(x0 + 1, originalWidth - 1);\n const y0 = Math.floor(origY);\n const y1 = Math.min(y0 + 1, originalHeight - 1);\n\n const dx = origX - x0;\n const dy = origY - y0;\n\n // Bilinear interpolation\n const value =\n grid[y0 * originalWidth + x0] * (1 - dx) * (1 - dy) +\n grid[y0 * originalWidth + x1] * dx * (1 - dy) +\n grid[y1 * originalWidth + x0] * (1 - dx) * dy +\n grid[y1 * originalWidth + x1] * dx * dy;\n\n scaledGrid[y * newWidth + x] = value;\n }\n }\n return scaledGrid;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":75,"value":"## Map Design","pinned":false,"hidden":false,"mode":"md"},{"id":862,"value":"For the final maps I chose to show the full archaeological dataset, but also to select some other groups of sites according to their tags. Unfortunately, because of inconsistencies with which OSM contributors choose to tag entries, this isn't a straightforward data query. But a reasonable approach seems to be simply to search for keywords in the entirety of the feature tags. This allowed me to extract those sites associated with `Roman` and those with `tumulus`. I would have liked to extract features associated with other eras, especially as these often have distinct geographical footprints, but matching them to obvious keywords was more challenging.\n\nTo assist in exploring the data, I added an interactive tooltip to each map, displaying all the stored tags with each site (_move the pointer over any map to see details_).","pinned":false,"hidden":false,"mode":"md"},{"id":723,"value":"roman = archaeology.filter((obj) => {\n return Object.values(obj.properties).some(\n (value) => typeof value === \"string\" && /roman/i.test(value)\n );\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":778,"value":"tumuli = archaeology.filter((obj) => {\n return Object.values(obj.properties).some(\n (value) => typeof value === \"string\" && /tumulus/i.test(value)\n );\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":193,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"

\n${drawMap(archaeology, { title: \"All records\" })}\n${drawMap(roman, { title: \"Roman\", scheme: \"Inferno\" })} \n${drawMap(tumuli, { title: \"Tumuli\", scheme: \"Cividis\" })} \n

","pinned":false,"hidden":false,"mode":"html"},{"id":866,"value":"_Side note: It would appear that no sites were returned for any of the Hebrides, Orkney or Shetland Isles (all rich in archaeological sites) despite my query specifying `Great Britain`. Not sure if this was a problem with my query script or tagged area in the database. _","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"---\n\n## Appendices","pinned":false,"hidden":false,"mode":"md"},{"id":7,"value":"import { TDMC } from \"@jwolondon/30-day-map-challenge-utilities\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":872,"value":"_For point-in-polygon testing:_","pinned":false,"hidden":false,"mode":"md"},{"id":77,"value":"T = await import(\"https://cdn.skypack.dev/@turf/helpers@7?min\").then(\n async (helpers) => {\n const [pointInPolygon] = await Promise.all([\n import(\n \"https://cdn.skypack.dev/@turf/boolean-point-in-polygon@7?min\"\n ).then((d) => d.default)\n ]);\n return { pointInPolygon };\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":870,"value":"_For KDE raster calculation:_","pinned":false,"hidden":false,"mode":"md"},{"id":264,"value":"kde = require(\"fast-kde\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":868,"value":"_For map title font:_","pinned":false,"hidden":false,"mode":"md"},{"id":737,"value":"","pinned":false,"hidden":false,"mode":"html"}],"files":"$W33","stdlib":"1","titling":"auto"}