31:T68c,md`I got an early morning email from our science editor, [Andrew Freedman](https://www.axios.com/authors/AFreedman), asking for a map to accompany a story he was writing about how the frequency of tornadoes had changed in different parts of the U.S. He needed it within a few hours so he could run it with his [Science newsletter](https://www.axios.com/newsletters/axios-science). The story would be based on [this article in Nature](https://rdcu.be/9rLd) (DOI: 10.1038/s41612-018-0048-2). I emailed the article’s authors asking for the data used to produce [their map](https://www.nature.com/articles/s41612-018-0048-2/figures/5). Usually, when scientists send me gridded geospatial data, they send a GeoTIFF or netCDF — formats you can simply drop into QGIS, style, and export. Instead, they sent me three .txt files. (This was my fault. One of the authors emailed me right away asking what format I'd like, but I didn't see his email because I was working on something else. About 45 minutes later, he emailed me again saying he'd be out of pocket for the rest of the day and attached the .txt files.) The three files were named lats.txt, lons.txt, and trend.txt. A quick look at the files showed what appeared to be thousands of space-separated values. I assumed (correctly, it turned out) that I could match the values from one file to the corresponding values in the other files to produce a dataset of thousands of observations, each one with the properties latitude, longitude, and trend (the value upon which I would base the color scale).`32:T619,md`This confirmed that I had processed the data correctly: each trend value had been assigned to the proper coordinate. However, I still had some work to do. In this dataset, the coordinates represent not specific locations, but rather centroids of contiguous cells. To fill those cells, I would have to calculate their boundaries. I decided to do this with a Voronoi diagram. I thought about using [d3-geo-voronoi](https://github.com/Fil/d3-geo-voronoi) for this job. But since I was only displaying a part of the globe, and since each cell was so small, I felt that the computational cost of calculating a spherical Voronoi tessellation was not worth the benefit. Instead, I used [d3-voronoi](https://github.com/d3/d3-voronoi) to calculate the planar Voronoi tesselation. Mike Bostock now recommends that you use [d3-delaunay](https://github.com/d3/d3-delaunay) for this sort of job, but I'm familiar with how d3-voronoi works — and I was still on a tight deadline. Author's note, May 18, 2019: I was under the false impression that d3-geo-voronoi was relatively slow because I had only attempted using an older version of it. As d3-geo-voronoi's creator Philippe Rivière [generously pointed out to me](https://observablehq.com/@fil/d3-geo-voronoi-and-gridded-data), d3-geo-voronoi now uses d3-delaunay internally, so it's probably even faster than d3-voronoi. Moreover, in its current verison, 1.2.1, d3-geo-voronoi can handle invalid input, so it could handle my dataset (which has NaN coordinate values) like a champ.`33:[{"id":"f3338f3780b75ee12650159e8c2af33ab0a8abd39dd5584fcb99a3fd4633b9a949053a93bb9929314bbe02109625135a58466af3b593cb12d0c553b2281ea155","href":"https://static.observableusercontent.com/files/f3338f3780b75ee12650159e8c2af33ab0a8abd39dd5584fcb99a3fd4633b9a949053a93bb9929314bbe02109625135a58466af3b593cb12d0c553b2281ea155","name":"states.json","type":"application/json","lastModified":1613593775362.75,"size":101495},{"id":"61d83ae8bd664874b9bb7f8b1a6da500c418818f01533b4055f30e3576d61897b43a6d7ef87059c07fa20a986e295aa9042be05604209079afb8a8d652ae0f02","href":"https://static.observableusercontent.com/files/61d83ae8bd664874b9bb7f8b1a6da500c418818f01533b4055f30e3576d61897b43a6d7ef87059c07fa20a986e295aa9042be05604209079afb8a8d652ae0f02","name":"lats.txt","type":"text/plain","lastModified":1628024621636.462,"size":750000},{"id":"1345779d1436772f15d4253d21dfd37d36dd8ce348237d0804b5cc9c80df35cc76228fbfcf6086ecfe28f3a8ce9849bac243e638079618c0d2ff919093b59278","href":"https://static.observableusercontent.com/files/1345779d1436772f15d4253d21dfd37d36dd8ce348237d0804b5cc9c80df35cc76228fbfcf6086ecfe28f3a8ce9849bac243e638079618c0d2ff919093b59278","name":"lons.txt","type":"text/plain","lastModified":1628024657429.113,"size":765000},{"id":"a6b11b9ce749353d91f4378599d3c2e0f3ad68721d442adaccbfaa4aaaa23e26fd47d6e41623f6e5bf9b902582af8c9ce7059163229bd018d2aa399ed207049d","href":"https://static.observableusercontent.com/files/a6b11b9ce749353d91f4378599d3c2e0f3ad68721d442adaccbfaa4aaaa23e26fd47d6e41623f6e5bf9b902582af8c9ce7059163229bd018d2aa399ed207049d","name":"trend.txt","type":"text/plain","lastModified":1628024687904.606,"size":750842}] 2f:{"title":"Mapping gridded data with a Voronoi diagram","cells":[{"id":0,"value":"md`# Mapping gridded data with a Voronoi diagram\n\nHere's the story of how I made this map for Axios, which was published on May 16, 2019, under the headline [“Climate change may subtly shift Tornado Alley.”](https://www.axios.com/climate-change-subtly-shift-tornado-alley-cda4c6f3-fdde-40d8-a9a1-c93c911695b4.html)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"html`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"rawdata = ({\n lats: parse(lats), lons: parse(lons), trend: parse(trend)\n});","pinned":false,"hidden":false,"mode":"ojs"},{"id":48,"value":"md`I was going to be displaying the data on a map, so I formatted my output as [GeoJSON](https://geojson.org/) rather than as an array of objects.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"geojson = rawdata.lats.map((lat, i) => {\n return {\n \"type\": \"Feature\",\n \"geometry\": {\n \"type\": \"Point\",\n \"coordinates\": [+rawdata.lons[i], +lat]\n },\n \"properties\": {\n \"trend\": +rawdata.trend[i],\n \"index\": i\n }\n } \n});","pinned":true,"hidden":false,"mode":"ojs"},{"id":86,"value":"md`I then placed each point on a U.S. map to see if the processing went okay.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n \n svg.selectAll(\".point\")\n .data(geojson)\n .enter().append(\"path\")\n .attr(\"d\", path);\n \n svg.append(\"path\")\n .datum(statesInner)\n .attr(\"class\", \"inner\")\n .attr(\"d\", path);\n \n svg.append(\"path\")\n .datum(statesOuter)\n .attr(\"class\", \"outer\")\n .attr(\"d\", path);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":166,"value":"md`So far so good: I had one datum per grid cell, and the data appeared to be positioned along lines of latitude and longitude. But to ensure that my data processing had assigned the correct trend value to each coordinate, I made a color scale and filled the dots based on the data.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":173,"value":"values = [-.63, -.56, -.49, -.42, -.35, -.28, -.21, -.14, -.07, 0, .07, .14, .21, .28, .35, .42, .49, .56, .63];","pinned":false,"hidden":false,"mode":"ojs"},{"id":170,"value":"colors = [\"#023858\", \"#045a8d\", \"#0570b0\", \"#3690c0\", \"#74a9cf\", \"#a6bddb\", \"#d0d1e6\", \"#ece7f2\", \"#fff7fb\", \"#fff\", \"#ffffcc\", \"#ffeda0\", \"#fed976\", \"#feb24c\", \"#fd8d3c\", \"#fc4e2a\", \"#e31a1c\", \"#bd0026\", \"#800026\"]","pinned":false,"hidden":false,"mode":"ojs"},{"id":176,"value":"colorScale = d3.scaleLinear()\n .domain(values)\n .range(colors)\n .interpolate(d3.interpolateLab);","pinned":true,"hidden":false,"mode":"ojs"},{"id":180,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n \n svg.selectAll(\".point\")\n .data(geojson)\n .enter().append(\"path\")\n .attr(\"d\", path)\n .style(\"stroke\", \"#000\")\n .style(\"stroke-width\", .1)\n .style(\"fill\", d => colorScale(d.properties.trend)); \n \n svg.append(\"path\")\n .datum(statesInner)\n .attr(\"class\", \"inner\")\n .attr(\"d\", path);\n \n svg.append(\"path\")\n .datum(statesOuter)\n .attr(\"class\", \"outer\")\n .attr(\"d\", path);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":184,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":191,"value":"voronoiGenerator = d3.voronoi()\n .x(d => d.x)\n .y(d => d.y)\n .size([width, height]);","pinned":true,"hidden":false,"mode":"ojs"},{"id":203,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n \n geojson.forEach(d => {\n const p = projection(d.geometry.coordinates);\n d.x = p ? p[0] : null;\n d.y = p ? p[1] : null;\n return d;\n });\n \n svg.selectAll(\".voronoi\")\n .data(voronoiGenerator(geojson).polygons())\n .enter().append(\"path\")\n .attr(\"class\", \"voronoi\")\n .style(\"stroke\", \"steelblue\")\n .style(\"stroke-width\", .5)\n .style(\"fill\", \"none\")\n .attr(\"d\", d => d ? `M${d.join(\"L\")}Z` : null);\n \n svg.selectAll(\".point\")\n .data(geojson)\n .enter().append(\"path\")\n .attr(\"d\", path)\n .style(\"stroke\", \"#000\")\n .style(\"stroke-width\", .1)\n .style(\"fill\", d => colorScale(d.properties.trend));\n \n svg.append(\"path\")\n .datum(statesInner)\n .attr(\"class\", \"inner\")\n .attr(\"d\", path);\n \n svg.append(\"path\")\n .datum(statesOuter)\n .attr(\"class\", \"outer\")\n .attr(\"d\", path);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":229,"value":"md`I now no longer needed the points because I could simply fill each Voronoi polygon with the color corresponding to its trend value.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":232,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n \n svg.selectAll(\".voronoi\")\n .data(voronoiGenerator(geojson).polygons())\n .enter().append(\"path\")\n .style(\"fill\", d => d ? colorScale(d.data.properties.trend) : null)\n .style(\"stroke\", d => d ? colorScale(d.data.properties.trend) : null)\n .attr(\"d\", d => d ? `M${d.join(\"L\")}Z` : null);\n \n svg.append(\"path\")\n .datum(statesInner)\n .attr(\"class\", \"inner\")\n .attr(\"d\", path);\n \n svg.append(\"path\")\n .datum(statesOuter)\n .attr(\"class\", \"outer\")\n .attr(\"d\", path);\n \n yield svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":237,"value":"md`That’s about it. At this point, I used [SVG Crowbar 2](https://nytimes.github.io/svg-crowbar/) so I could finish the rest of the graphic — cropping the map, creating a color legend, and adding a headline — in Adobe Illustrator.\n\nAll of that could have been done with code, of course, but finishing it up in AI let me export my map as an image so it would load as fast as possible. (Actually, I rendered it with [ai2html](http://ai2html.org/), which also ensured the text would scale with different screen sizes.)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":80,"value":"md`## Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":268,"value":"md`Below, you’ll find more information about data processing and styling.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":99,"value":"md`### The raw data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":106,"value":"md`The three .txt files I got from the scientists, lats.txt, lons.txt, and trend.txt:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"lats = FileAttachment(\"lats.txt\").text()","pinned":false,"hidden":false,"mode":"ojs"},{"id":50,"value":"lons = FileAttachment(\"lons.txt\").text()","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"trend = FileAttachment(\"trend.txt\").text()","pinned":false,"hidden":false,"mode":"ojs"},{"id":103,"value":"md`I wrote a function to parse each .txt file into a flat array.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":64,"value":"parse = text => text.split(\" \").map(d => d.trim())","pinned":false,"hidden":false,"mode":"ojs"},{"id":126,"value":"md`### Geospatial data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":242,"value":"md`A [TopoJSON](https://bost.ocks.org/mike/topology/) file of the U.S. states.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":98,"value":"states = FileAttachment(\"states.json\").json();","pinned":false,"hidden":false,"mode":"ojs"},{"id":245,"value":"md`I used topojson.mesh to make separate MultiLineStrings for the states’ inner and outer boundaries so that I could style them separately.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":129,"value":"statesInner = topojson.mesh(states, states.objects.states, (a, b) => a !== b);","pinned":false,"hidden":false,"mode":"ojs"},{"id":132,"value":"statesOuter = topojson.mesh(states, states.objects.states, (a, b) => a === b);","pinned":false,"hidden":false,"mode":"ojs"},{"id":248,"value":"md`I used [d3-geo](https://github.com/d3/d3-geo) to set up an [Albers USA projection](https://github.com/d3/d3-geo#geoAlbersUsa) and a [path generator](https://github.com/d3/d3-geo#geoPath) to render my geospatial data to the DOM.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":136,"value":"projection = d3.geoAlbersUsa()\n .fitSize([width, height], statesOuter);","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"path = d3.geoPath()\n .projection(projection)\n .pointRadius(2);","pinned":false,"hidden":false,"mode":"ojs"},{"id":149,"value":"md`### Styling`","pinned":false,"hidden":false,"mode":"ojs"},{"id":300,"value":"md`I styled the inner and outer boundaries separately.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":151,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":297,"value":"md`For U.S. maps on an Albers projection, I've found that the height is about .58 times the width.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":120,"value":"height = width * .58;","pinned":false,"hidden":false,"mode":"ojs"},{"id":198,"value":"md`### Libraries`","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"d3 = require(\"d3-geo@3\", \"d3-interpolate@3\", \"d3-scale@4\", \"d3-selection@3\", \"d3-voronoi@1\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}