Handling Raster Map Backgrounds in D3.js

One of the most popular blog post on the blog is a tutorial on making a map in d3.js. I thought it’s time to revisit the post and propose some visual improvements to it, starting with realistic backgrounds. My all time favorite interactive article from the NYT is their 2014’s How the Air Campaign Against ISIS Grew. The subtle use of terrain formation in the map’s background not only makes the graphic more interesting visually, but places the depicted countries in a real geographical context. That’s a more fitting visual vehicle for geo-political events than plain, almost cartoonish, vector-based maps. I thought my original map would benefit from this treatment too: in the following sections I will take you through generating a raster terrain background for my base vector map.

This is where we will get to:

Vector and raster data, unified

The image is enriched with a hillshade background. Note the lighter shading for Spain, the country the map is centered on. I copied that from the NYT: it differentiates the visualisation subject from its neighbours and draws the eye towards the most important part of the image.

To construct this map, you will need the following tools and source files:

  • GDAL
  • a map of administrative boundaries of the target region
  • a DEM file with hillshade information of the target region
  • a list of cities and their geocoordinates within the target region
  • pre-canned html and css files to get started
  • (optional) Photoshop

The first section guides you through the setup.

Tools & Data Sources

#1 GDAL

We will use GDAL for operations on raster and vector data. GDAL is a swiss-knife-like utility for geospatial transformations. It can be used as a command line tool or a Python package.

If you already have Python installed, getting GDAL is as simple as running:

pip install gdal; 

You can also install GDAL from source files: follow KyngChaos for directions on Mac, or UCLA’s guide for Windows installation.

#2 Administrative boundaries

I will be using a map of world administrative boundaries in 1:10M scale obtained from Eurostat. The file we are after is in the 1:10 scale and in a GeoJSON format. It will download as ref-countries-2020-10m.geojson.zip. The package includes the same map in various projections, and with various geo features. Let’s use the variant of administrative boundaries in a standard Web Mercator projection, in other words EPSG:3857: CNTR_RG_10M_2020_3857.geojson (direct download link).

#3 Elevation map

Now onto the terrain image for the background. The European Environment Agency is a great place to look as it publishes digital elevation data of Europe. Various detail levels and file types are offered: here I decided for the hillshade model with 1 km X 1 km resolution. The package will include some metadata, and the image we will use as our base: hillshade1x1.tif (direct download link).

#5 A list of cities

This one is a bit of a hack: I will reuse the list of highly-populated European cities I prepared in my 2018 post, Point-in-polygon, a mathematical cookie-cutter. The location dataset is based on the Geonames database, reduced to 50 top cities within Europe. You can follow my original article to generate the dataset yourself, or directly download the cities.csv file.

#5 Basic html and css

To make things easier and focus on the task at hand, I will reuse the map generated in the original tutorial. I reviewed the base map code and the color scheme; other than that, it’s the same map if you feel like taking a step back and following that tutorial first.

Base vector map, centered on Spain

Copy the html and css files and save them as map.html and map.css.

map.html:

<!DOCTYPE html>
<html lang="en">
    <head>
        <meta charset="utf-8">
        <title>My map</title>
        <script type="text/javascript" src="https://d3js.org/d3.v5.min.js"></script>
        <link rel="stylesheet" type="text/css" href="map.css">
        <style></style>
    </head>
    <body>
        <div id="container" class="svg-container"></div>
        <script type="text/javascript">
        const w = 900;
        const h = 600;
        const svg = d3.select(selector).append("svg").attr("preserveAspectRatio", "xMinYMin meet")
        .attr("viewBox", "0 0 " + w + " " + h)
        .classed("svg-content", true);
        const projection = d3.geoMercator().translate([w/2, h/2]).scale(2200).center([0,40]);
        const path = d3.geoPath().projection(projection);

        // load data  
        const worldmap = d3.json("countries.geojson");
        const cities = d3.csv("cities.csv");

        Promise.all([worldmap, cities]).then(function(values){   
            // boundary box

            // scale

            // translate

            // position and scale

            // draw background

            // draw map
            svg.selectAll("path")
                .data(values[0].features)
                .enter()
                .append("path")
                .attr("class","continent")
                .attr("d", path),
            // draw points
            svg.selectAll("circle")
                .data(values[1])
                .enter()
                .append("circle")
                .attr("class","circles")
                .attr("cx", function(d) {return projection([d.Longitude, d.Lattitude])[0];})
                .attr("cy", function(d) {return projection([d.Longitude, d.Lattitude])[1];})
                .attr("r", "1px"),
            // add labels
            svg.selectAll("text")
                .data(values[1])
                .enter()
                .append("text")
                .text(function(d) {
                            return d.City;
                    })
                .attr("x", function(d) {return projection([d.Longitude, d.Lattitude])[0] + 5;})
                .attr("y", function(d) {return projection([d.Longitude, d.Lattitude])[1] + 15;})
                .attr("class","labels");
        });
    </script>
    </body>
</html>

map.css:

.svg-content {
    background-color: #ecf3f4;
}

.continent {
    fill: #fcf7e9;
    stroke: #434141;
    stroke-width: 0.2;
}

.circles {
    fill: #111104;
}

.labels {
    font-family: Georgia, serif;
    font-size: 14px;
    fill: #111104;
}

Merging, clipping, and projecting

In the source data preparation phase as well as in the drawing part of the tutorial, it’s necessary to keep in mind the objective of all this work. It sounds sort of obvious, but without the destination in mind we won’t pave the road as efficiently. To recap then: the goal here is to produce two maps, one vector and one raster, in the same projection and covering the same region.

How do we get about realizing the idea? Firstly, let’s establish the region we want to visualize and get its boundaries in a GeoJSON format. Next, let’s bring the two maps to a common projection, e.g. Web Mercator. Once aligned with its vector equivalent, the raster map will need to be clipped to match the target geographical area. Additional image coloring and format transformation can then take place.

#1 Extracting the target region

Since Spain has been selected as the map’s main point of interest, we will need to extract some of its neighbours so it doesn’t hang there in space on its own. Selecting a handful of countries in a close approximation to Spain is nicely achieved with ogr2ogr (part of GDAL): just specify the abbreviations of the countries you need and the tool will extract them. You can look up their FID values directly in the GeoJSON file or by using the information view in MapShaper.

To extract Portugal, Spain, Andorra, France, Switzerland, Italy, Germany, Hungary, Lichtenstein, Belgium, Luxembourg, Slovenia, Croatia, Slovakia, Bosnia and Herzegovina, Czechia, Algeria, Morocco, Tunisa, and Libya, type the following in your command line:

ogr2ogr -where "FID in ('PT','ES','AD','FR','CH','IT','AT','DE','HU','LI','BE','LU','SI','HR','SK','BA','CZ', 'DZ','MA','TN','LY')" countries.geojson CNTR_RG_10M_2020_3857.geojson

All selected polygons will be exported to the countries.geojson file.

#2 Bringing the maps to a common projection

The essential part of this visualization is ensuring that the raster image and the vector file are configured in the same spatial reference system. We know that the country boundaries have been projected in EPSG:3857 because of the impeccable naming standards applied by Eurostat. EPSG:3857 is a Mercator-based projection standard for the Web, popularized by Google and used by Mapbox, Bing, and OpenStreetMap project; it’s also supported out of the box in D3.js. As our map is Spain-centric, the Mercator projection works very well here with its clear depiction of Europe so I would suggest it for the final product, too.

We need to find out what the hillshade’s projection is and, if needed, align it with the EPSG:3857 projection. The projection can be found out with the gdalsrsinfo utility, courtresy of GDAL. Let’s see what we are dealing with here. In the terminal, run:

gdalsrsinfo hillshade1x1.tif

You should see similar output:

Gdalsrsinfo

Projection information about the hillshade file

Cool! We learn that the projection is based on Lambert Azimuthal Equal-Area CRS, which is a very different approach to Mercator. We also get a precise measurement to replicate the customization applied by the European Environmental Agency: PROJ.4 : +proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +R=6378137 +units=m +no_defs.

Hillshade

The hillshade image before the re-projection

This allows us to project the file to the Web Mercator CRS. This is achieved using gdalwarp, as below:

gdalwarp -s_srs "+proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +R=6378137 +units=m +no_defs" -t_srs EPSG:3857 hillshade1x1.tif hillshade-projected.tif

0…10…20…30…40…50…60…70…80…90…100 - done! (the command line will say)

We passed the source projection, the target projection, the input and the output file names to the tool. In return we get a hillshade-projected.tif file, in EPSG:3857 SRS.

Hillshade projected

The hillshade image after the re-projection to Web Mercator

#3 Aligning the raster’s region

Now that the files exist in the same coordinate system, we can align them so that they represent the same part of the world. The raster area will be adjusted to match the extent of the globe covered by the GeoJSON collection of countries. Here’s how to achieve this:

First, we need to compute the boundary box of the GeoJSON file. In geographical terms that’s the polygon spread across the furthest north-west located point and the furthest south-east coordinate. The boundary box, or the vector’s extent, can be found out by parsing the GeoJSON with the ogrinfo utility:

ogrinfo countries.geojson -so -al | grep Extent

The following extent is displayed:

Extent: (-7022952.000000, -2438344.000000) - (6215653.000000, 7372407.000000)

The extent is returned in a (xMin,yMin) - (xMax, yMax) format. That’s the mask we will apply on the raster file to cut it to the exact area covered by the vector map. The clipping is done with the gdal_translate utility with the projwin option. You will notice that the projwin’s coordinates are passed in a different order than those given by extent - it’s a bit confusing - but that’s the convention used by projwin. extent’s output is thus reshuffled to xMin, yMax, xMax, yMin to match projwin’s ulx uly lrx lry order.

This will clip the hillshade image with the boundary box of the GeoJSON shapes:

gdal_translate -projwin -7022952 7372407 6215653 -2438344 hillshade-projected.tif hillshade-box.tif

We receive a hillshade-box.tif file from GDAL:

Hillshade clipped

The hillshade image clipped to the region covered by the vector map

Clipping is perhaps the wrong word here as the image seems to have shifted around! There is no need to worry: the output looks weirdly deformed, because the geographical area it covers includes the Canary Islands, the Azores, and Madeira (parts of Spain and Portugal!) on the South-West, and French Overseas Territories in the South-East.

#4 Coloring the map

The question on everyone’s mind by this point is why the hillshade image is all black: wasn’t it supposed to show some, um, shades?

What you’re seeing is a standard photo editor’s interpretation of a TIFF file. A typical image viewer will have a hard time making sense out of the information stored in a TIFF, especially if no red/blue/green integer bands are available.

If you run a gdalinfo on our TIFF file, somewhere mid-way in the output you will see the following color information:

Band 1 Block=4901x1 Type=Int16, ColorInterp=Gray
Description = Band_1
Min=0.000 Max=254.000 
Minimum=0.000, Maximum=254.000, Mean=169.161, StdDev=43.820
NoData Value=-32768

This is saying that our raster image is coded as grayscale and extends up to 255 values (translated as shades of grey), and that the pixels where the measurement could not be taken are represented by a dummy NoData value, set to -32768.

The good news is that can re-color the image and assign it RBG bands. You can either use Photoshop for this (check out this tutorial to learn how to create index maps for TIFF files), or continue with the tools provided by GDAL. I will stick with GDAL in this example.

All we need for re-coloring is the range values that can be used by the image (as we learnt that’s 0-254), and the color codes we want to translate them to. I played around with some colors, and came up with this palette:

0,245,245,244
160,240,240,239
200,214,213,209
254,118,117,112

0 stands for the lowest altitude, 254 for the highest. My mapping colors the lowest elevation regions in light gray and the highest in darker, brownish-gray. I also put some additional color points in the middle to break out the scale. Save the file as colormap.txt, and run the following command to apply the new color scheme to the map:

gdaldem color-relief hillshade-box.tif colormap.txt hillshade-color.tif

gdaldem allows for specifying other cool things, such as the location of Sun over the globe, which would influence how the shades are appearing. Our bare minimum changes will result in the following image (hillshade-color.tiff.):

Hillshade color

The hillshade image after recoloring

You will notice that the areas previously drawn in pitch black are milky white now. This is because ALL values are translated to the new palette, also the NoData values. This has its benefits, as it colors those spots within the map that were marked as missing measurement, and would have stayed transparent or black otherwise. Especially by coloring the missing areas inland, it provides a smooth visual experience for the viewer.

Running gdalinfo on the new image will now show 3 bands: one for red, one for green, and one for blue. This or adjusting the color map to a more vivid palette than my conservative choice should be proof enough that we are now in the RGB realm.

#5 Cutting out the target area

Now that the image areas have the same proportions, and we like the new color scheme, we can proceed to cutting out only the country shapes that we want to keep in the final visualization. It’s like clipping cute dog pictures from magazines and leaving all the text behind. Yes, that’s the best comparison I came up with. Here we will leave out all countries we don’t need and all water regions. You bet this can be easily done with gdalwarp too!

gdalwarp -cutline countries.geojson hillshade-color.tif hillshade-cutout.tif

Neat! We get hillshade-cutout.tif after the command executes. The area matching the GeoJSON shape is distilled from the original image.

Hillshade cutout

Selected countries are cutout from the hillshade image

#6 Saving the raster as a png image

We reach the stage where its safe to take off the spatial engineer hat, and don that of a webpage designer. Here we bastardize the raster by saving it as a png and thus drop all geographical information associated with it. Now that the image has been reprojected, clipped, and colored to nicely accompany the GeoJSON shapes, there is no benefit of keeping the additional metadata.

gdal_translate -of PNG hillshade-cutout.tif europe.png

This produces a europe.png image that can be referenced by D3 in the next section.

Drawing the map

With two files spanning over the exact same geographical area, we have almost everything needed to build a coherent visualization. The secret sauce that will put it all together is the correct scaling and centering. The strategy here is to see how the files relate to the svg, and then use this information to transform them accordingly. Currently on our to-do list:

  • Calculate the boundary box of the vector shape stored in the GeoJSON file
  • Calculate the ratio of the vector’s dimensions against the svg
  • Calculate the central point of the scaled vector
  • Resize the raster image

#1 The raster meets the vector

Our starting point is the map.html and it’s supporting style sheet, map.css. Let’s make place for the raster by setting the continent color to transparent. Adjust the css file to:

.continent {
    fill: none;
    stroke: #434141;
    stroke-width: 0.2;
}

It will produce the following underwater visualization of Europe:

The continent fill is set to transparent

Time to add the hillshade map to the background. First, we’ll hard code its width and height so it can be referenced later. Place this under the load data section:

const raster_w = 8234/10
const raster_h = 6102/10

Note that the raster is resized only for the purpose of fitting it nicely on the svg - it’s too large otherwise - there is no technical reason to do it. Later we will drop the hard coded values.

Let’s append it to the svg. Paste this within the draw background section:

svg.append("image")
        .attr('id', 'europe')
        .attr("xlink:href", "europe.png")
        .attr("width", raster_w)
        .attr("height", raster_h);

That’s it. After refreshing the screen we should see the map emerge in the background:

The raster image is appended to the svg

The two files are not aligned, and that’s what we will be fixing in the next section.

#2 Scale and Transform

As announced earlier, we will be scaling the vector against the available svg estate. To achieve this, remove the set projection completely, so that we get undistorted vector lenghts.

Replace this line:

const projection = d3.geoMercator().translate([w/2, h/2]).scale(2200).center([0,40]);

With the following:

const projection = d3.geoMercator().translate([0,0]).scale(1);

The next step is to compute the boundary of the GeoJSON shape. This will be used to scale it accordingly. D3 comes with a method to do just that, d3.path.bounds(). The offical documentation gives the following description to path.bounds: “returns the projected planar bounding box (typically in pixels) for the specified GeoJSON object. The bounding box is represented by a two-dimensional array: [[x₀, y₀], [x₁, y₁]], where x₀ is the minimum x-coordinate, y₀ is the minimum y-coordinate, x₁ is maximum x-coordinate, and y₁ is the maximum y-coordinate.”

Boundingbox

An illustration of path.bounds(), coutresy of StackOverflow

Let’s compute the boundary of the GeoJSON. Paste the following lines under the boundary box section:

const b = path.bounds(values[0]);
console.log(b);

Save the file, and observe what happens.

The vector’s projection is temporarily removed

The vector is placed in the left hand side corner of the svg. Since we didn’t feed it any geographical information it cannot place itself on the plane. If you open up the console log you will see how the boundary was calculated:

Boundary box

The boundary box of the vector

Now onto bulding a scale.

The idea here is to calculate the ratio of how the GeoJSON translates to the svg: this will give us a scale we can then apply to the raster image. The svg size is 600x900, so we need to calculate how the vector’s height compares to the svg’s height, then do the same for the width. We then take the max of the returned values to see how many times would the vector fit in the svg. This is represented by s in the script.

Variable t centers the image to the svg’s center. See how I’m using the scale on the vector’s total height and width to project it on the svg plane. The totals are halved to return the center coordinate of the GeoJSON.

Finally, we pass the scale and the offset to the projection.

Paste the following snippet in the indicated spots:

// scale
const s = 0.99 / Math.max((b[1][0] - b[0][0]) / w, (b[1][1] - b[0][1]) / h); 

// transform
const t = [(w - s * (b[1][0] + b[0][0])) / 2, (h - s * (b[1][1] + b[0][1])) / 2];

// update projection
projection
    .scale(s)
    .translate(t);

Refresh the page. The vector boundaries are now set to their target spot.

The vector is scaled and shifted

The raster image will get a slightly different treatment. The image has been cut to cover the exact same area as the GeoJSON, hence scaling it is as simple as assigning it the dimensions of the scaled vector. The image doesn’t use any projection (it’s a basic png), so we’ll replace the width and height attributes directly and shift the image on the screen usingthe transform method. The following rows should go under scale and position section:

// scale and position
const raster_width = (b[1][0] - b[0][0]) * s;
const raster_height = (b[1][1] - b[0][1]) * s;

const rtranslate_x = (w - raster_width) / 2;
const rtranslate_y = (h - raster_height) / 2;

And this will replace the old width and height attributes:

.attr("width", raster_width)
.attr("height", raster_height)
.attr("transform", "translate(" + rtranslate_x + ", " + rtranslate_y + ")");

Save and refresh:

The vector and the raster are fixed to their target locations on the svg

That’s right, the vector and the image are aligned! The visualization clearly needs some final touches, but we are almost there.

Prettifying the map

First things first: let’s zoom in to center on Spain. I tried various attributes and ended up with the following viewBox setup. Replace the original values with the following:

.attr("viewBox", "370 70 " + w/5 + " " + h/5)

We are catapulted to the Mediterranean!

Adjusting the viewBox property the map is zoomed in

Some minor css adjustments, and we’re almost there.

Now that the image is zooomed in 5 times, we need to make some adjustments. Reduce the point radius to 0.2px, and the distance of the labels to the points to 1px up, and 1px right.

.attr("r", "0.4px")

.attr("x", function(d) {return projection([d.Longitude, d.Lattitude])[0] + 1;})
.attr("y", function(d) {return projection([d.Longitude, d.Lattitude])[1] - 1;})

Replace the css body completely with the following snippet:

.svg-content {
    background-color: #ecf3f4;
}

.continent {
    fill:#fafafa;
    stroke: #434141;
    stroke-width: 0.2;
    opacity: 0.6;
    stroke-opacity: 0.3;
}

.circles {
    fill: #111104;
}

.labels {
    font-family: Georgia;
    font-size: 15%;
    fill: #555555;
}

#ES {
    fill: #fafafa;
    opacity: 0.8;
    stroke: #434141;
    stroke-width: 0.2;
    stroke-opacity: 0.3;
}

It will calm down the fonts and add a transparent mask to the whole continent. Spain will be additionally distinguished by a slightly lighter color of its area.

To enable that, we need to be able to access the country names (or ids). Update the draw map section as below:

svg.selectAll("path")
        .data(values[0].features)
        .enter()
        .append("path")
        .attr("class","continent")
        .attr("d", path)
        // assign ids
        .attr("id", function(d) {
            return d.properties.FID
        }),

Save all your changes and refresh the map. This is where our journey ends. Thank you for reading!

The final visualization

Further reading

To build this tutorial I have consulted the following articles and examples:

Appendix

If you followed the tutorial, at the end of it your html page will look like this:

map.html:

function draw_final(selector){
    const w = 900;
    const h = 600;

    const svg = d3.select(selector).append("svg").attr("preserveAspectRatio", "xMinYMin meet")
    //.attr("viewBox", "0 0 " + w + " " + h)
    .attr("class","svg")
    .attr("viewBox", "370 70 " + w/5 + " " + h/5)
    .classed("svg-content", true);

    const projection = d3.geoMercator().translate([0,0]).scale(1);
    const path = d3.geoPath().projection(projection);

    const worldmap = d3.json("countries.geojson");
    const cities = d3.csv("cities.csv");

    Promise.all([worldmap, cities]).then(function(values){
        // boundary box
        const b = path.bounds(values[0]);

        // scale
        const s = 0.99 / Math.min((b[1][0] - b[0][0]) / w, (b[1][1] - b[0][1]) / h);  

        // transform
        const t = [(w - s * (b[1][0] + b[0][0])) / 2, (h - s * (b[1][1] + b[0][1])) / 2];

        // update projection
        projection
            .scale(s)
            .translate(t);

        // scale and postion
        const raster_width = (b[1][0] - b[0][0]) * s;
        const raster_height = (b[1][1] - b[0][1]) * s;

        const rtranslate_x = (w - raster_width) / 2;
        const rtranslate_y = (h - raster_height) / 2;

        // draw map
        svg.append("image")
            .attr('id', 'Raster')
            .attr("clip-path", "url(#europe_color.png)")
            .attr("xlink:href", "europe.png")
            .attr("class", "raster")
            .attr("width", raster_width)
            .attr("height", raster_height)
            .attr("transform", "translate(" + rtranslate_x + ", " + rtranslate_y + ")");

        svg.selectAll("path")
            .data(values[0].features)
            .enter()
            .append("path")
            .attr("class","continent")
            .attr("d", path)
            // assign ids
            .attr("id", function(d) {
                return d.properties.FID
            }),

        // draw points
        svg.selectAll("circle")
            .data(values[1])
            .enter()
            .append("circle")
            .attr("class","circles")
            .attr("cx", function(d) {return projection([d.Longitude, d.Lattitude])[0];})
            .attr("cy", function(d) {return projection([d.Longitude, d.Lattitude])[1];})
            .attr("r", "0.4px"),
        // add labels
        svg.selectAll("text")
            .data(values[1])
            .enter()
            .append("text")
            .text(function(d) {
                return d.City;
            })
            .attr("x", function(d) {return projection([d.Longitude, d.Lattitude])[0] + 1;})
            .attr("y", function(d) {return projection([d.Longitude, d.Lattitude])[1] - 1;})
            .attr("class","labels");
    }); 
}