31:T4a3,intro = md`Let's start with the static example map below, which attempts to illustrate connectivity *and* earthquake sensitivity across the realtime station network maintained by UNAVCO. First a [Delaunay tesselation](https://en.wikipedia.org/wiki/Delaunay_triangulation) is created with **${d3.format( "," )( data.features.length )}** *realtime* station positions as its vertices. We simply take \`(longitude, latitude)\` pairs and project them onto an orthographic map (a globe) that you will be able to rotate and zoom interactively futher down this page. Then, the resulting triangles in the tesselation we created are filled in with a color representing the minimum earthquake magnitude for which the selected amount of displacement (e.g. 1 cm) would be measurable by at least one of the station/vertices making up a particular triangle. The static example below shows what happens when you hover your mouse around station **P013** in the interactive globe further down the page. The magenta radii show coverage around neighboring stations for the station closest to your mouse pointer. The blue radius shows the range of station P013 below if your mouse is closest to P013.`32:Te01,chart = { const context = DOM.context2d(width, height); context.font = "14px roboto"; const path = d3.geoPath(mutable projection, context).pointRadius(1); context.cursorStyle = "crosshair"; mutable projection.rotate([110, -45]); function render({ lofi = false }) { context.clearRect(0, 0, width, height); context.globalAlpha = 1; context.beginPath(); path(sphere); context.lineWidth = 1.5; context.fillStyle = "#cfe0ff"; context.strokeStyle = "#000"; context.fill(); context.stroke(); context.beginPath(); path(graticule); context.lineWidth = 0.5; context.strokeStyle = "#aaa"; context.stroke(); // land context.beginPath(); path(world); context.fillStyle = "#faf8f5"; context.fill(); context.stroke(); // States context.beginPath(); path(stateShapes); context.fillStyle = "#faf8f5"; context.fill(); context.lineWidth = 1.5; context.stroke(); context.lineWidth = 0.5; if (!lofi) { // Delaunay triangles triangles.features.forEach((feature) => { context.beginPath(); path(feature); context.fillStyle = d3.interpolateTurbo( magnitudeScale(d3.max(feature.properties.magnitudes)) ); context.globalAlpha = d3.max(feature.properties.magnitudes) < magnitudeCutoff ? magnitudeScale(d3.max(feature.properties.magnitudes)) : 0; context.strokeStyle = "#aaa"; context.stroke(); context.fill(); }); context.globalAlpha = 0.3; context.beginPath(); path({ type: "MultiPoint", coordinates: points }); context.fillStyle = "#000"; context.fill(); if (mutable highlightedIndex > -1) { d3.merge([ [mutable highlightedIndex], delaunay.neighbors[mutable highlightedIndex] ]).forEach((n, i) => { let radius = cutoffDistance({ magnitude: neighborMagnitude, displacement: displacement }).MelgarEtAl2015_2D; context.beginPath(); context.globalAlpha = 1; path( d3 .geoCircle() .center(points[n]) .radius(greatCircle2ArcDegrees(radius))() ); context.strokeStyle = i == 0 ? "blue" : "magenta"; context.fillStyle = i == 0 ? "blue" : "magenta"; context.setLineDash([5, 5]); context.lineWidth = 1.5; context.stroke(); context.globalAlpha = 0.1; context.fill(); context.setLineDash([]); context.fillStyle = "#000"; context.globalAlpha = 1; context.fillText( data.features[mutable highlightedIndex].properties.site_code, mutable mouse[0] + 10, mutable mouse[1] - 10 ); }); } } } let transform; function zoomed(e) { transform = e.transform; mutable projection = mutable projection; render({ lofi: true }); // Render fewer details while user moves map. } context.canvas.onmousemove = (event) => { if (transform) { let coordinate = mutable projection.invert(d3.pointer(event)); mutable mouse = d3.pointer(event); mutable highlightedIndex = delaunay.find(...coordinate); } render({ lofi: false }); }; return d3 .select(context.canvas) .call( zoom(mutable projection) .on("zoom.render", zoomed) .on("end.render", () => render({ lofi: false })) ) .call(render) .node(); }33:T514,md`We can of course consider the simple **heatmap** that results from drawing overlapping circles with a radius indicating the detection area around each station for a given displacement and magnitude. We can skip all of that triangulation business above and quickly draw our overlapping radii. Darker areas have more overlapping areas of detection from multiple stations. These darker areas can be considered more sensitive than lighter areas as we can use GNSS signals and ground stations to detect earthquakes of smaller magnitudes in these areas. However, the triangulation approach simplifies the otherwise very computationally intensive problem of determining how much overlap there is between each circle as the number of sites increases. Triangulating the network using the Delaunay tesselation allows for a topological data structure to be created and quickly queried. Overlapping circles, on the otherhand, looks cool but does not lend itself to quickly computing quantitative properties about the overlapping pile. For example, once the tesselaton is created, the time to find the neighbors around each vertex is constant time \`O(1)\` and that allows for snappy user interaction! *Try moving the sliders below to see how sensitive the network is to particular magntitude earthquakes.*`34:T98f,chart2 = { const context = DOM.context2d(width, height); context.font = "14px roboto"; const path = d3.geoPath(mutable projection, context).pointRadius(1); context.cursorStyle = "crosshair"; mutable projection.rotate([110, -45]); function render({ lofi = false }) { context.clearRect(0, 0, width, height); context.globalAlpha = 1; context.beginPath(); path(sphere); context.lineWidth = 1.5; context.fillStyle = "#cfe0ff"; context.strokeStyle = "#000"; context.fill(); context.stroke(); context.beginPath(); path(graticule); context.lineWidth = 0.5; context.strokeStyle = "#aaa"; context.stroke(); // land context.beginPath(); path(world); context.fillStyle = "#faf8f5"; context.fill(); context.stroke(); // States context.beginPath(); path(stateShapes); context.fillStyle = "#faf8f5"; context.fill(); context.lineWidth = 1.5; context.stroke(); context.lineWidth = 0.5; if (!lofi) { let radius = cutoffDistance({ magnitude: magnitudeRadius, displacement: displacement }).MelgarEtAl2015_2D; context.globalAlpha = 0.03; // Delaunay triangles data.features.forEach((feature) => { context.beginPath(); path( d3 .geoCircle() .center(feature.geometry.coordinates.slice(0, 2)) .radius(greatCircle2ArcDegrees(radius))() ); context.strokeStyle = "magenta"; context.fillStyle = "magenta"; context.fill(); context.stroke(); }); context.beginPath(); path({ type: "MultiPoint", coordinates: points }); context.fillStyle = "#000"; context.fill(); } } let transform; function zoomed(e) { transform = e.transform; mutable projection = mutable projection; render({ lofi: true }); // Render fewer details while user moves map. } context.canvas.onmousemove = (event) => { if (transform) { let coordinate = mutable projection.invert(d3.pointer(event)); mutable mouse = d3.pointer(event); mutable highlightedIndex = delaunay.find(...coordinate); } render({ lofi: false }); }; return d3 .select(context.canvas) .call( zoom(mutable projection) .on("zoom.render", zoomed) .on("end.render", () => render({ lofi: false })) ) .call(render) .node(); }35:[{"id":"b4e5eaa6af400e91c92474a94b8b503c6effe0e137a9f2b4b6b423fc800ae854e7313fef3511302931b8bfa9e834a7abc270442f662f642c2b33fa38e99112f7","href":"https://static.observableusercontent.com/files/b4e5eaa6af400e91c92474a94b8b503c6effe0e137a9f2b4b6b423fc800ae854e7313fef3511302931b8bfa9e834a7abc270442f662f642c2b33fa38e99112f7","name":"Screen Shot 2022-06-02 at 12.48.28 PM.png","type":"image/png","lastModified":1654188574676.034,"size":259622}] 2f:{"title":"EarthScope Network Sensitivity","cells":[{"id":0,"value":"title = md`# EarthScope Network Sensitivity`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1086,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1447,"value":"

\n \n

The station P013 has the radius of detection shown in purple for an earthquake magnitude 5.8 and a peak ground displacement of 1 cm. The neighboring circles show the closest neighbors to P013, where neighbors are those stations which share a cell boundary in the Vornoi diagram underlying the Delaunay Tesselation drawn as a mesh.

\n","pinned":false,"hidden":false,"mode":"html"},{"id":1475,"value":"md`## Interactive Network Sensitivity Map (with triangles)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1495,"value":"explanation = md`*Interact with the map below by clicking and dragging and zooming using your trackpad or mousewheel... or pinch and zoom with your fingers on mobile.*\n\n*Hover over the map to highlight a station and its closest neighbors (based on a [Voronoi diagram](https://en.wikipedia.org/wiki/Voronoi_diagram#:~:text=In%20mathematics%2C%20a%20Voronoi%20diagram,%2C%20sites%2C%20or%20generators).) created with stations as its vertices). Try adjusting the \\` Event Magnitude\\` slider to see the areas around each station where an earthquake of this magnitude would be detectable.*`","pinned":false,"hidden":false,"mode":"ojs"},{"id":599,"value":"viewof displacement = Inputs.range([0.005, 0.2], {\n value: 0.025,\n step: 0.001,\n label: \"Minimum Measurable Displacement (m)\" // meters\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":550,"value":"viewof magnitudeCutoff = Inputs.range([1, 9.1], {\n value: 4.2,\n step: 0.1,\n label: \"Show Triangles Up To Magnitude\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1190,"value":"viewof exponent = Inputs.range([0.1, 10], {\n value: 2,\n step: 0.1,\n label: \"Adjust Color Scaling\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1534,"value":"md`The globe below shows areas where a magnitude **${magnitudeCutoff}** earthquake could be detected above a displacement threashold of **${\n displacement * 100\n}** centimeters.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1069,"value":"legend = Legend(\n d3\n .scaleSequentialPow(\n [Math.min(triangles.max, 9.1), Math.max(triangles.min, 5.9)],\n d3.interpolateTurbo\n )\n .exponent(exponent),\n {\n title: \"Magnitude required for detection\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1478,"value":"md`## Interactive Heatmap`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1433,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1424,"value":"viewof magnitudeRadius = Inputs.range([1, 9.1], {\n value: 4.2,\n step: 0.1,\n label: \"Event Magnitude\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1409,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1389,"value":"md`---\n## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1374,"value":"mutable mouse = []","pinned":false,"hidden":false,"mode":"ojs"},{"id":1261,"value":"mutable highlightedIndex = -1","pinned":false,"hidden":false,"mode":"ojs"},{"id":1353,"value":"neighborMagnitude = magnitudeCutoff","pinned":false,"hidden":false,"mode":"ojs"},{"id":1331,"value":"highlightedNeighbors = delaunay.neighbors[highlightedIndex]","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"mutable projection = d3.geoOrthographic()\n .fitExtent([[1, 1], [width - 1, height - 1]], sphere)\n .rotate([0, -30])","pinned":false,"hidden":false,"mode":"ojs"},{"id":57,"value":"height = Math.max(640, width / 2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"world = {\n var world = await (\n await fetch(\"https://unpkg.com/world-atlas@1/world/110m.json\")\n ).json();\n return topojson.feature(world, world.objects.countries);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":59,"value":"sphere = ({type: \"Sphere\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"graticule = d3.geoGraticule10()","pinned":false,"hidden":false,"mode":"ojs"},{"id":340,"value":"data = metadataSearch({\n realtime: \"Y\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"points = data.features.map((feature) =>\n feature.geometry.coordinates.slice(0, 2)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":480,"value":"delaunay = d3.geoDelaunay(points)","pinned":false,"hidden":false,"mode":"ojs"},{"id":703,"value":"voronoi = d3.geoVoronoi()(points)","pinned":false,"hidden":false,"mode":"ojs"},{"id":540,"value":"function greatCircleDistance(radians) {\n return (radians / (2 * Math.PI)) * EARTH_CIRCUMFERENCE;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":572,"value":"function greatCircle2ArcDegrees(distance) {\n return (distance / EARTH_CIRCUMFERENCE) * 360;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":537,"value":"EARTH_CIRCUMFERENCE = 40074.275 // km","pinned":false,"hidden":false,"mode":"ojs"},{"id":706,"value":"triangles = {\n let included = voronoi.triangles(points);\n included.min = Infinity;\n included.max = 0;\n included.features.forEach((feature) => {\n let T = feature.geometry.coordinates[0];\n feature.properties.magnitudes = [];\n\n // For each edge of the triangle\n for (let i = 0; i < 3; i++) {\n let m = 1;\n let distance = greatCircleDistance(d3.geoDistance(T[i], T[i + 1]));\n let range = 0;\n\n while (distance > 2 * range) {\n m += 0.1;\n range = cutoffDistance({\n magnitude: m,\n displacement: displacement\n }).MelgarEtAl2015_2D;\n feature.properties.magnitude = m;\n }\n\n if (m < included.min) included.min = m;\n if (m > included.max) included.max = m;\n\n feature.properties.magnitudes.push(m);\n }\n });\n\n return included;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":944,"value":"magnitudeScale = d3\n .scalePow()\n .exponent(exponent)\n .domain([triangles.max, triangles.min])\n .clamp(true)\n .range([0, 1])","pinned":false,"hidden":false,"mode":"ojs"},{"id":1325,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":216,"value":"import {zoom} from \"@d3/versor-zooming\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":343,"value":"import { metadataSearch } from \"@unavco/web-services\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":385,"value":"import {\n peakGroundDisplacement,\n cutoffDistance,\n pgdVersion\n} from \"@unavco/peak-ground-displacement\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":423,"value":"import { stateShapes } from \"@enjalot/us-county-datasets\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"d3 = require(\"d3@6\", \"d3-geo-voronoi@1\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1071,"value":"import { Legend } from \"@d3/color-legend\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}