October 23, 2019

Mapping Swiss coordinates (LV95) with d3: Part 2

This is the second of a two-part series about creating a map in SVG based on Swiss geographical data. Each part contains a certain step of the process:

  1. Converting shapefiles to GeoJSON/TopoJSON with node
  2. Creating a bivariate choropleth map with d3-geo and LV95 coordinates

What we are going to end up with at the end of the second part is this:

result

A bivariate choropleth map of Switzerland showing how income is distributed over Swiss municipalities.

Where we left off at part one: We finished creating a TopoJSON file that was at reasonable file size. We now take that TopoJSON and create an svg map out of it.

Working with d3-geo

I won’t go into great detail about the basics of d3-geo. There are plenty of resources (here, here or here). Instead, we’ll focus on how Swiss coordinates are different to work with.

The Swiss coordinate system has one key difference to others: it is already flat. What do I mean by that? While lat and long describe degrees on a globe, the Swiss coordinates X and Y describe a position on a normal plane:

cartesian

So with most coordinate systems, you need a projection like geoMercator to map coordinates from a sphere to a flat screen. This is not needed if you work with Swiss coordinates.

To draw our municipalities onto our svg we need the following ingredients: First, we convert the TopoJSON back to geojson. This is done using the topojson-client package.

import { feature } from 'topojson-client'

// convert topojson back to geojson
const geoJson = feature(topoJson, topoJson.objects.municipalities)

TopoJSON can contain multiple objects (we could for example also store cantonal boundaries in the same file). When converting back to GeoJSON, we need to specify which of the objects we want to read. The name municipalities was introduced when we created the TopoJSON on this line here: topojson.topology({ municipalities: geojson } (read more about that process in Part 1).

Next, we need to set up two classic d3 linear scales for x and y. To get the extent of the two directions, we use a small helper that returns us the bounding box of a GeoJSON object: @turf/bbox. @turf has also some other neat helpers, check them out when you work with geo data.

import bbox from '@turf/bbox'

// get extent of switzerland as x, y coordinates
const [minX, minY, maxX, maxY] = bbox(geoJson)

// calculate aspect ratio and derive height
const height = ((maxY - minY) / (maxX - minX)) * width

const x = scaleLinear()
  .range([0, width])
  .domain([minX, maxX])

const y = scaleLinear()
  .range([0, height])
  .domain([maxY, minY])

// this is where the magic happens
// read more about the background here:
// https://bl.ocks.org/mbostock/6216797
const projection = geoTransform({
  point: function (px, py) {
    this.stream.point(x(px), y(py))
  }
})

// this function will create svg paths when
// passed a geojson feature
const path = geoPath().projection(projection)

The most important part here is what happens inside of geoTransform. The function is a helper that lets us construct our own projection function. (Another example by mbostock here).

In there we just return for every coordinate px and py the corresponding pixel position generated by our linearScales x and y (Thx to @lucguillemot for pointing this out to me).

To plot the map, we could also just use plain d3, the code for that would look something like this:

svg
  .selectAll('path')
  .data(geoJson.features)
  .enter()
  .append('path')
  .attr('class', 'feature')
  .style('fill', 'steelblue')
  .attr('d', path)

But as we needed both some npm packages as well as a way to import JSON files, I decided to set up a small react-based repository.

In react, the code to draw our municipalities looks as follows:

<svg
  width={width}
  height={height}
>
  {geoJson.features.map(feature =>
    <path
      key={`path-${feature.properties.bfsId}`}
      stroke='white'
      strokeWidth={0.25}
      d={path(feature)}
      fill='steelblue'
    />
  )}
</svg>

The color of each municipality is just steelblue at the moment but yay 😊 we have a map – and it uses LV95 coordinates!

one colored map

Meet the data

From the website of the Swiss Federal Tax Administration, I downloaded a table that contains both a Gini coefficient and a mean income for every municipality in Switzerland.

We want to create a map that shows both where the income is high or low and where the income is distributed more or less equal amongst the inhabitants of a municipality.

The bivariate color scale

A bivariate color scale will help us do exactly that: it can show two different things at the same time. In our case: Whether a municipality has a low / medium / high inequality in income distribution (red axis) and whether it has a low / medium / high average income (blue axis). When both of these are low, the municipality will be plotted in grey, when both are high, it will be plotted in violet:

bivariate scale

You can find more information on how to create bivariate color scales in this blog post.

Using the “Sketch” software we combined two scales with fill “blue” (base color #1E8CE3) and “red” (base color #C91024), each with three different opacities (20%, 50%, 80%), and blend mode “multiply”, which resulted in the 3 * 3 = 9 hex values.

This GIF from Joshua’s blog linked above nicely shows the process of blending the two scales:

To replicate this process now on our map with d3 we need to create not one scale, but three. The first two are two d3 threshold scales that return either low, medium or high when given a mean income per municipality or a Gini coefficient (a number used to measure inequality).

const meanScale = scaleThreshold()
  .domain([35000, 45000])
  .range(['low', 'medium', 'high'])

const giniScale = scaleThreshold()
  .domain([0.4, 0.5])
  .range(['low', 'medium', 'high'])

Example: meanScale(20000) will return low as the number 20’000 is below the first threshold 35’000. And giniScale(0.6) will return high. I did not chose the thresholds randomly but I fist calculated three quantiles and then rounded the numbers so they’re more human friendly to read.

The third scale will be responsible for deriving a real hex color value for the nine possible categories that a municipality can be in:

const bivariateColorScale = scaleOrdinal()
  .domain([
    'high high',
    'high medium',
    'high low',
    'medium high',
    'medium medium',
    'medium low',
    'low high',
    'low medium',
    'low low'
  ])
  .range([
    '#3F2949', // high income, high inequality
    '#435786',
    '#4885C1', // high income, low inequality
    '#77324C',
    '#806A8A', // medium income, medium inequality
    '#89A1C8',
    '#AE3A4E', // low income, high inequality
    '#BC7C8F',
    '#CABED0' // low income, low inequality
  ])

We use a d3 ordinal scale for that. When called with the argument high high for example, the function will return violet: #3F2949.

So we map over our data array and add the hex color code to every item:

const dataWithFill = data.map(d => {
  const meanGroup = meanScale(d.mean)
  const giniGroup = giniScale(d.gini)
  const fill = bivariateColorScale(`${meanGroup} ${giniGroup}`)
  return { ...d, fill }
})

When we use these color codes for our <path /> elements we get the following result:

result no relief

That’s already pretty cool. But there’s a lot of red! A lot of municipalities in the Swiss Alps have a quite high income inequality. That is interesting, of course, but as only some hundred people are living in most of these municipalities, it’s better if we can give them a little less visual weight.

We’ll do that by only displaying the so-called “productive areas” of each municipality. Thx to the Swiss Federal Statistical Office who provided us these boundaries that you’ll find in the files called gde-1-1-15.shp in the public github repo.

We will underly our map with a relief of the swiss alps so the “unproductive” areas are not left white. The creation of the relief image file was done in R. You can read more about that process in the blog post that my colleague Timo Grossenbacher and I wrote where we create this very map using R, sf and ggplot.

And tadaa:

result

We end up with the map you saw at the very top. I also added a small legend at the top left corner to explain to the viewer what the colors mean.

This map was published on srf.ch where I added some more details to it like annotations or tooltips.

If you want to know how I did this exactly, let me know via twitter at @angelozehr. I could cover this in a third part of this tutorial then.


Angelo Zehr

Written by Angelo Zehr, data journalist at SRF Data and teacher.


Further reading