Using QGIS and Mapbox to Map Tornado Hotspots on bunker.land

I’ve spent the last couple weekends putting together bunker.land, which I posted about a couple weeks ago (the tl,dr is, “mapping out the riskiest places in the US to live”). 

My focus has been on targets during a nuclear war, but I thought it would be fun to expand the project to include natural disaster likelihoods.  I didn’t have too much trouble finding datasets for elevation (to model sea level rise), earthquakes, and wildfires. Tornado risk seemed like a good next step.

I wanted a map that looked something like this (source):

(aside:  I grew up in the middle of Missouri, which according to this map, is literally the armpit of tornado alley.  And yet in the 15 years I lived there, I never saw even one tornado, and I am pretty salty about this.  Now I know why. More on this later.) 

However, I needed the tornado hazard data formatted as GIS shapefiles or rasters so I could render it via Mapbox GL JS, the library I use to display webmaps.  Sadly, I had a lot of trouble finding a GIS-formatted risk map for tornadoes  The closest thing I found was a set of tornado starting points from 1950-2017.  This is a comprehensive dataset, but when I pulled it into QGIS and mapped it out, the raw data was a bit… noisy:

Since I couldn’t find a map that worked out of the box, I had no choice but to learn something new.  Luckily, I found a guide for making heatmaps in QGIS, which gave me a really good starting point. Less fortunately, the guide is for an old version of QGIS, and as a result I hit a number of obstacles I had to Google around.  

I’m pretty happy with the result, and spent a fair amount of time learning how to use QGIS, so I figured I’d write up how I filtered this data into vector datasets and got it into Mapbox, where I display it as a layer on bunker.land.

Making a heatmap

Starting from the very beginning, we’ll want a new QGIS project.  So we have some visual context when playing around with these shapes, I added an OpenStreetMap base layer.  The tornado dataset we want to work with is available as shapefiles, and we can add that to QGIS via Layer → Add Layer → Add Vector Layer:

Our lives will be easier layer on if this data is all reprojected into EPSG:3857 – WGS 84 before we do any editing.  We can just do that first. Right click on the layer → Export → Save Features As:

Set the CRS to WGS 84, save, and we can work with the new reprojected layer from now on.

So our first problem is that this dataset is huge and noisy.  While I don’t recommend ignoring any tornadoes, I would not personally get off my couch for anything less than an F2 tornado, so that’s what I’m going to filter for.  Since this data is a shapefile, I can filter on the fields of the objects; right click on the new layer → Filter.   

We’ll just filter on the “mag” column, looking for F2+ tornadoes:

This is a bit less cluttered, but still not super actionable.  From here, our goals are:

• turn these points into a heatmap
• extract discrete layers from the heatmap
• save the extracted layers as shapefiles

Luckily for us, QGIS has a nifty heatmap tool which lets us turn our points into a heatmap raster.  Click on Processing → Toolbox → Interpolation → Heatmap:

Iterating on the settings here took a while; I had to experiment before I found settings that looked good.  I went with a 150km radius on the heatmap points, 4000 rows, and 10009 columns (once you select the number of rows, the columns auto-populate).  I played around with the colors on the resulting heatmap for a bit and ended up with this:

Intensity bands

While this is an improvement, it’s still kind of a hot mess (pun intended).  Heatmaps based on spotty data like this probably over-exaggerate the hotspots (there’s likely reporting bias, and we don’t want to overweight individual data points).  I’d prefer to get discrete intensity bands. To get those, we can use the raster calculator: Raster → Raster Calculator:

Since our heatmap values are no longer really connected to any actual unit, choosing units was a bit of guesswork.  Frankly, I just chose numbers that lined up with the risk areas I saw on other maps; the lowest gradient, 10, gives us this:

This is the kind of gradient band I’m interested in.  Unfortunately, this is still a raster image. We really want shapefiles — we can do more interesting things with them on Mapbox and in the browser, and the data is dramatically smaller.  Luckily, QGIS has tool to turn raster images into shapefiles: “polygonalize”. We can go Raster → Conversion → Raster to Vector:

We can select whatever we’ve named our filtered raster.  This gives us the previous image broken into two chunks:

We want to filter out the part that falls below our heatmap threshold.  Right click the layer → Properties → Filter:

Filter for where the feature value is equal to 1.  Now we’re down to the shapes we care about:

Of course we can play around with the layer styling to get it to look like whatever we want:

To capture the gradients we care about, we can repeat this process at a few thresholds to capture distinct bands.  These don’t correspond to any particular intensity, they are just intended to demarcate more and less intense risk areas.  

Fast-forwarding the repetitive bits, I’ve repeated these steps with four raster calculator thresholds (with this dataset, I ended up using thresholds of 10, 25, 40, and 65).  By setting a different color on each layer I’ve produced and decreasing opacity to 50%, I got this:

This captures what I want; distinct gradient bands without overly-weighting hotspots.  If your goal is just to generate a static raster image, you can stop here and export this image directly. 

Mapbox

My goal however is to import these layers into Mapbox so I can attach them to an existing interactive web map.  Mapbox is a platform for hosting customized maps and embedding them in apps or webapps; I use Mapbox, plus the corresponding Mapbox GL JS library, to host maps for bunker.land.  To get this data into Mapbox, we want to upload the data as a Tileset and use the data within a Style as feature layers.

I learned the hard way that there is a good way to do this, and a bad way to do this.  The simple way is to export each of the 4 bands as a GeoJSON file, upload it to Mapbox, and add it as a layer.  This is a mistake. Mapbox has a limit of 15 data “Sources” per Style, so saving each layer as a separate GeoJSON file and uploading them separately quickly caps out how many layers we can have per Style.

Luckily, Mapbox has released a nice tool called tippecanoe which lets us combine GeoJSON files into a single mbtiles file (it can do a ton of other things too; this is just what I’ve used it for).  An mbtiles file can have as many layers as we want, as long as it is under 25 GB.

First we want to extract each layer as a GeoJSON file; right click the layer → Export → Save Features As.

Choose GeoJSON and repeat for each layer.  This gives us four geojson files:

$ ls -lh *.geojson
-rw-r--r--  1 bpodgursky  640K Aug  1 22:46 tornado10.geojson
-rw-r--r--  1 bpodgursky  590K Aug  1 22:45 tornado25.geojson
-rw-r--r--  1 bpodgursky  579K Aug  1 22:45 tornado40.geojson
-rw-r--r--  1 bpodgursky  367K Aug  1 22:44 tornado65.geojson

We can use tippecanoe to combine these into a single, small, mbtiles file:

$ tippecanoe  -zg -o tornado.mbtiles — extend-zooms-if-still-dropping *.geojson
$ ls -lh tornado2.mbtiles
-rw-r--r--  1 bpodgursky  128K Aug  1 22:54 tornado.mbtiles

This gives us a single tornado.mbtiles file.  

In practice I added these layers to an existing map for bunker.land; for simplicity, here I’m going to set up a new empty Style.  After setting up a Mapbox account, navigate to Studio → Styles → New Style.  I use a blank background, but you can also choose an OpenStreetMap background.

We can add these layers directly to the Style.  Navigate through Add layer → Select data → Upload to upload the mbtiles file we just generated.  These features are small and should upload pretty quickly.  Once that’s available (you may need to refresh), we see that there are four layers in the new source:

We’ll create four new layers from this source.  We’ll just use the Mapbox studio to recreate the styling we want, and set the opacity so the overlay is visible but doesn’t obscure anything:

All I needed to do now was get this into a website.

Embedding on bunker.land

Mapbox GL JS has great examples about how to get a Style in a map, so I won’t dig into the code too much; the important part is just loading a map from this style:

mapboxgl.accessToken = YOUR_TOKEN;

var map = new mapboxgl.Map({
  container: 'map', // the div we want to attach the map to
  style: 'mapbox://styles/bpodgursky/cjxw0v4fr7hd81cp6s0230lcw', // the ID of our style
  center: [-98, 40], // starting position [lng, lat] -- this is about the middle of the US
  zoom: 4 // starting zoom level
});

We can see the final result here, overlaid against OpenMapTiles on bunker.land:

Since our layer is just a simple vector tile layer, it’s easy to detect these features on-click for a particular point, along with any other enabled layers:

Wrapping up

It’s now pretty clear why I missed all the tornadoes as a kid — Tornado Alley (kind of) skips right over central Missouri, where I grew up!  My only explanation for this is, “weather is complicated”.

On the technical side, I was surprised how easy it was to generate a decent-looking map; Mapbox and QGIS made it stupidly easy to turn raw data into a clean visualization (and I’ve only been using QGIS for a couple weeks, so I’m sure I missed a few nice shortcuts.) 

Now that I know how to turn ugly data into nice heatmaps or gradient data, I’ll probably work on adding hurricanes and flooding over the next couple weeks.  Stay tuned.