Cartographic palettes and colour harmonies

This story begins one day when I was assembling a map of the city of Edmonton, Alberta from OpenStreetMap data. It was going to be a big map, a 42″ (106 cm) wide poster for a wall.

OpenStreetMap colour palette

The data was good, but the standard OSM colours were not. They would work fine for a street map but, for a piece of art being seen all day in someone’s living space, there was an unpleasant amount of grey. Plus, the accents of pinky red and ochre, and the cold green parks, gave it a lack of colour harmony.

If you analyze the colour profile of this map, counting the number of pixels of each colour and assuming there are 12 colours being used, you can see that white and the first two greys occupy about 70% of the map surface. You’d not be wrong if you said, “That’s a grey map.”

Moreover, this colour palette does have any of the classic “harmonies” that we expect when we look for how colours work together in an image. (E.g., see here.) There’s no dominant point on the colour wheel, no complementary colour, no triadic distribution or analogous set. It’s not OpenStreetMap’s fault: this is just the colour mix that you get when you look at this specific framing of this city.

If we want the map to look harmonious, we need to think about two things. What colours are used? And what proportions are they used in?

For example, look at the 1986 map, below, of Central Asia from the Royal Geographic Society. My sense here is that the cartographers thought about the colours, their relationships to one another, and the percentages of the map they would occupy. Most of the colours are clustered around the green-yellow-orange side of the colour wheel, but there are a few complementary blues. There is enough green in India to showcase the yellows, oranges and browns of Tibet. The adjacent blobs of orange Tibetan highland and grey-brown Tarim Basin work well together. There’s a nice balance of values from very dark to white. There’s no particularly bright colour that draws the attention from other areas. It says “reference map,” rather than “map with a point to make.” Overall, it’s quite pleasant to look at.

Image Color Summarizer

It turns out that there’s an interesting tool for doing an initial analysis of a colour palette and the proportions of an image the different colours occupy. It’s called Image Color Summarizer and it comes to us courtesy of Martin Krzywinski at the BC Genome Sciences Centre. You go to http://mkweb.bcgsc.ca/color-summarizer/?home and go to the Analyze tab. Upload your image and choose, say, 12 colour clusters and very high precision.

The software then reduces the number of colours in the image to (in this case) 12 , using an algorithm that shifts colours a minimum to fit them into 12 categories. It then gives you all sorts of great information. So here we have first an analysis of colours in the RGS map of Central Asia…

…and then my OSM map of Edmonton:

This is very useful for a quick look at colour proportions and interrelatedness between colours. You get each colour in hex, RGB, HSV and LCH, and the percentage of the image that it occupies.

You might have noticed, however, that Image Color Summarizer gave different results on my OSM map from the ones I showed at the beginning of this page. It ranked the two greys first, and then white, with the three totalling only 60% of the image. There are a number of reason for this, but one is that this software does have a shortcoming: even at the “very high” precision setting it reduces your image to 200 pixels wide before doing its analysis. This means that if you have tiny features, like small labels or trim lines running alongside roads, their contribution to the overall colour mix can be lost.

If you want to use Image Color Summarizer fairly often, I recommend you download the software (on the Download tab) and run it locally at the command line. It’s a PERL script. This both reduces the load on the BCGSC server that hosts it, and allows you to analyze images larger than 200 px wide. You can just use your own computing resources and just let it run all night! (Think about it, though: an image 1500 px wide will take 100 times longer to analyze than one which is 150 px wide.)

At this point you might think the problem is basically solved: we extracted a 12-colour palette and their proportions from a map whose colour harmony I like, and all we have to do is to apply those to my OSM map.

But it’s not so simple. First of all, the OSM map did not have a 12-colour palette. In fact, according to GIMP (Color>Info>Colorcube Analysis), it had 213,078 colours. And this shouldn’t be surprising: OSM Downloader generated scores of rules for styling its many lines and polygon types. And then when I exported the map out of QGIS, antialiasing produced even more colours.

Just a few of the polygon styling rules from OSM Downlaoder

Image Color Summarizer suggested a palette, but before it did that it had scaled the map image to a mere 200 pixels across, and then simplified it to 12 colours.

Working with the actual map (some 12600 pixels across) there are two possible ways we can go:

  1. We could take the finished map with its OSM colours, and a “model” image, simplify them both to a reasonable number of colours (say 12), and then do a 1-to-1 substitution based on the order of colours when they are ranked by percentage of image. Or…
  2. We could simplify the styling of the map within QGIS to, say, 12 colours, and then apply the colours we get from simplifying the model.

Let’s look at each of these in turn.

Replacing the colour palette of a finished map, using GIMP

This is certainly the easier way, but it may not give you the level of control you want.

GIMP has an excellent tool for simplifying an image to a set number of colours at Image > Mode > Indexed. It allows you to decide how many colours you want (we’ll use 12), and whether these can be dithered or not. Generally I do not use dithered, except sometimes I do.

Once you have reduced your map to 12 colours, you can see those colours by opening the Colormap dialogue (Windows > Dockable Dialogs > Colormap).

Unfortunately they are not ordered by how much of the image they represent. You have to determine this yourself, and it’s a bit clumsy. Have the Histogram dialogue open (Windows > Dockable Dialogs > Histogram) and right click each colormap colour in turn, taking “Select this colour.” The Histogram window will tell you how many pixels in the image are selected. I take notes into a spreadsheet and that gives me the table I showed above:

The top line, “Index” tells me the position of the colour in the original colormap, from 0 to 11. I then re-order the colourmap in descending order of “real estate” (that is, what percentage of the image is occupied by each colour), by right-clicking any colour and choosing “Rearrange colormap.” Now my colourmap looks like this:

At this point, white, the most common colour, is index 0; the two greys are indices 1 and 2; two greens are indices 3 and 4; and so on.

I do the same tedious process with my model image. I reduce it to 12 colours, I make notes on how much of the image each colour represents…

…and re-order its colormap correspondingly:

In order to substitute the new palette colours for the original colours I have to do one more preparatory step. I have to convert this GIMP colourmap into a GIMP palette. I do this by opening the Palette dialogue (Windows > Dockable Dialogs > Palettes), right-clicking it and taking “Import Palette.” I import the palette using the model image as source.

Finally, I can apply that palette to the original map by going to Colors > Map > Set Colormap, and choosing the palette I just imported.

And there it is.

RGS “Mountains of Central Asia” colour palette

Now, you can see right away that this is both successful and not exactly what we expected. The map has acquired the palette of the original image, and by a stroke of luck the North Saskatchewan River became a pale green, which actually works for a river. It looks pretty good! But, to be picky, it doesn’t exactly have the original colour proportions. As noted before, the RGS map has quite a bit of green in it (anchoring the lower left), but we didn’t get that much green.

The reason for this becomes apparent when we compare the original map image (left) with its indexed version (right)…

Yes, there were a lot of green pixels there, but the indexing process turned much of it into patches of yellow. This is a hazard of reducing the number of colours in an image.

Still, that was successful enough that I want to try extracting a colour palette and proportions from some non-map image, and seeing if I can apply those to my map. Let’s try—I don’t know—this Georgia O’Keefe painting of an iris:

Once reduced to 12 colours (I did use dithering here) the palette in descending order of pixels looks like this:

And the result looks like this:

Georgia O’Keefe “Light Iris” colour palette

Again, it’s basically right, although the overall colour signature isn’t quite as light as the original image. The palette is centred around many varieties of purple and white, with a bit of nearby blue, opposed by a small proportion of complementary yellows. Looking at the original Georgia O’Keefe painting, I’d expect a bit more of the bright yellow, but I see that this colour went into the pixels of small features like riverside bike paths, only visible when you get up close. I’d also expect a bit more black, but the black got shunted into the railways lines so is quite fragmented and not visible when you look at the map as a whole.

This points out another factor: the extracted colour palette and proportions only tell part of the story of how an image presents itself. Another important factor is whether a given colour is fragmented into many small pixel groups, or collected into one or two large swathes. Maps (at least urban maps) tend to have a lot of pixels dedicated to small labels or street lines, and the colours used here rarely declare themselves when you look at the whole map.

Designing the map to a specific palette from the start

Swapping in the colours from an image with good colour harmony, in order of prevalence, sort of works, but now let’s look at our other alternative. This is simplify the styling of the map within QGIS to, say, 12 colours, and then changing these colours before exporting the map.

One advantage of this was that I would be able to apply a new palette to my map carefully. I wanted, for instance, to be able to put dark colours where I needed dark colours (labels, high contrast borders), to pair related colours together (parks and their sports fields, or neighbourhoods and their buildings) and to deploy complementary colours to occupy features that would allow them to give just the right amount of zing to the overall impression.

I also wanted to avoid this sort of thing, where the label and its halo in the original OSM map were of contrasting values, but after GIMP swapped in a different colour palette they had similar values.

Text was indexed colour “8” in the original map (left), which was a dark, warm grey. Its halo was colour “5”, a light tan, and this resulted in good contrast. In the new palette (right), colour “8” was a dark blue, and colour “5” was an only slightly lighter teal; as a result, there was very little contrast.

QGIS has an excellent feature for building a map from a specific colour palette: Project Colors. This is under Project>Properties>Default Styles. Essentially, you define a colour palette, each colour having a text string to identify it.

These palettes can be saved (and imported) as .GPL text files. Then, when styling a feature, you don’t assign it just any old colour: you use a data-defined override to assign it one of the colours from the Project Colors list:

With this in place, when the Project Color with that name changes, all features tied to it change colour as well.

Unfortunately, QGIS does not automatically pick up the colours of a project and list them for you. I still went through the process of indexing the image in GIMP and entering those colours as project colours. I named them by the typical features that used them. Note that once you begin assigning project colours to features, you can’t change the names of the colours (e.g., “building, street casings”) because in the QGIS project file the colours used in feature styling are actually referred to by these names.

It took quite a while to go through all of the styling rules and style everything using only the 12 project colours. Pleasingly, it did not change the overall look of the map very much.

Then I made an alternate project palette from the 12 Georgia O’Keefe colours, not worrying too much this time about what percentage of the original image they occupied, but assigning them where I thought they would work best. It seemed like it might be time to let go of applying colours in the order of their prevalence, because it was not necessarily producing good effects.

And that produced this map:

Georgia O’Keefe “Light Iris” colour palette, version 2

It has to be said that this method gives much better control. Given the colours you want in your palette, you can move them around to different classes of features until they represent the right proportion of the image. You can alter them slightly. You can save palettes you like and then read them back in later. (Note: when you read in a .GPL colour file, QGIS appends them to the project colours. Then you have to select and then delete the previous project colours .)

Furthermore, once the initial work of styling the map to a certain number of colours is done (and I soon expanded to 18 project palette colours, since 12 is pretty tight!), the process of going from a new model image…

“Early winter”

…to its palette…

…to a new version of the map…

“Early Winter” colour palette

… is quite fast. What made this quick was that, after GIMP indexed the photograph to 18 colours, I took a screenshot of its colormap, rotated that 90°, and then used the QGIS Colour Picker tool to assign these colours to the project palette. (And then saved that palette.)

A final piece of colour palette inspiration was this map from the Second Military Survey of the Habsburg empire, Galicia and Bucovina (1861–1864), which I obtained by screenshot from arcanum.com.

I love its muted but warm bronze tones, the pale green for flatlands, and the thin blue stream running through it. (This is the area today around Verkhovyna/Верховина, Ukraine). It also has a bit of complementary colour: some barely noticeable red elevation notations.

GIMP simplified this map to these 18 colours, which do capture the overall palette, but miss the red and the blue:

It seemed reasonable to just add them in:

And that produced this:

“Second Military Survey of the Habsburg Empire” colour palette

In summary

I’m sold on the second method: designing a map around a specific palette, and then swapping in alternate palettes with good colour harmonies stolen from other images. This was the approach that gave more satisfying results. You do more work up front, but in the process you also come to understand the underlying colour structure of your map: which colours are adjacent to which, which need to be related, which are fragmented and which occur in large continuous fields.

An 18-colour palette was much more reasonable than 12-colour palette. For this kind of urban map I calculated that I needed:

  • one colour with four different values. In the original map these were the greys: residential zones, buildings, neighbourhood labels.
  • other light-dark pairs, maybe six of these. For example you need a light blue and dark blue (bike routes and water), a light red and a dark red (motorway fill, motorway casing), a light yellow and a dark yellow (secondary road fill, secondary road casing), etc.
  • black
  • white

The GIMP is an invaluable tool in assessing the colours of an image. If nothing else, it can help you quantify the “colour signature” of an image.

If the overall colour harmony of your map matters (and it should!), palette-cloning and palette-switching are powerful tools. You could tweak your map so its colours reflect the natural colours of the area depicted. You could harmonize a map’s colours with those of the larger display in which it is being placed. You could deliberately echo the colours of an earlier map.

Advertisement

Making shaded relief directly from DEMs projected in degrees?

In a presentation I gave at GeoIgnite 2021, I was explaining the process of using QGIS to make shaded relief from DEMs, more or less as detailed in this post. There wasn’t time to explain the process of re-projection, so I rather boldly asserted that you don’t need to re-project your DEMs before making shaded relief. Instead, I said, you can just use the magic scaling number 111120 to convert from degrees to metres.

Really? Is this true?

Well, yes, in many cases it works pretty well. In QGIS you pull in an unreprojected DEM (its native CRS using degrees as units) but set the Project CRS to some reasonable metric projection, such as UTM. Then, opening the GDAL hillshade tool from the Toolbox, you use the DEM as your Input layer, set scale to 111120, and enter the parameters of your choice for Z factor, Azimuth of the light and Altitude of the light.

And, look: shaded relief:

This seems so much simpler than having to first re-project your DEM to square cells measured in metres. But there are some caveats. This technique fools with the terrain a bit, and could potentially distort how it looks. If you understand how it distorts it, you’ll be able to evaluate if you want to use it in your map.

The magic number

The magic number 111120 comes from the GDAL documentation, specifically the documentation for the gdaldem tool:

The gdaldem generally assumes that x, y and z units are identical. If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, the scale (-s) option can be used to set the ratio of vertical units to horizontal. For LatLong projections near the equator, where units of latitude and units of longitude are similar, elevation (z) units can be converted to be compatible by using scale=370400 (if elevation is in feet) or scale=111120 (if elevation is in meters). For locations not near the equator, it would be best to reproject your grid using gdalwarp before using gdaldem.

https://gdal.org/programs/gdaldem.html

Gdaldem is the tool that I recommend using for creating hillshades in QGIS, and there in the docs they say it: unless your locale is near the equator, it’s best to re-project your grid before using gdaldem.

The obvious problem with entering a scale factor of 111,120 metres per degree is that as you move away from the equator, lines of longitude become closer to one another. For example, in Canada at 54° north, one degree of latitude is close to 111,0000 metres, but a degree of longitude is only about 65,000 metres. In other words, a DEM cell is not square on the ground. So when you enter 111,120 as the scaling factor to convert degrees to metres, you will get a hillshade calculated on the basis of a square cell that has been created by stretching a rectangular cell in an east–west direction.

For example, with a 1″ DEM, each cell is 0.0002777777° (that’s 1 ÷ 3600) on a side. If we scale both axes at 111,120, the hillshading algorithm thinks it has a square cell 30 metres on a side. However the reality at 54° north is that on the ground the cell is about 30 m north–south and 18 metres east–west.

How bad is it?

You’d expect that because the landscape is stretched east–west, the slope would be under-calculated for cells that slope east–west. Basically, these cells would appear to be on a slope less steep than it really is.

And slope is important in the hillshading algorithm: it plays a role in determining how bright a hillshade pixel is. Pixels whose aspect and slope indicate that they are orthogonal to the sun direction are the brightest white, while those facing away from the sun direction are darkest black.

So if, for example, you had a 45° slope facing west, and the sun in the west at an elevation of 45 °, you’d expect pixels on the slope to have the highest possible hillshade value, 255. But if you calculated the hillshade from a degree-DEM using a scale of 111,120, that slope would be under-estimated; maybe it would appear to be 35°. As a result it would not get assigned a 255, but something less, maybe 240. There might be a nearby area, also facing west, whose true slope was 55°, but it now gets assigned the brightest value because it appears to be 45°.

In general then, for east or west-facing surfaces, we’d expect the hillshade built from the un-re-projected degree-DEM to distribute its brightest and darkest patches differently, particularly when the sun is due east or west.

That’s just what you see.

However, when the sun is not due east or west it’s hard to see this difference in slope calculation.

You can also experiment with other, lower scale factors, such as 90,000, 80,000 or 70,000. These, in effect, cause slopes to be less under-estimated on east–west aspects, but over-estimated on north–south aspects. In other words you are distributing the error around.

This suggests that in many cases, making a hillshade directly from a DEM projected in degrees is fine.

Other questions

A second distortion arising from creating a hillshade using a degree-DEM and a scale factor of 111,120 concerns the sun’s azimuth. If you are mapping in a local UTM projection, you know that UTM grid north is usually different from true north. This was visible in the figure above with the contour lines: the cell that was “square” in degrees turns out to be rectangular in a UTM projection, and also a bit tilted to the side.

With UTM projections, this difference is not more than 3°. At worst, then, you might ask for a sun azimuth of 315° and instead get one of 312°. I doubt whether you can notice this difference.

The last problem we might wonder about is whether shadows calculated by the hillshading algorithm are distorted. Are they foreshortened in an east–west direction when the hillshade is re-projected into a metric projection? The answer to this may be surprising: the GDAL hillshading algorithm does not calculate shadows. It calculates shadowed surfaces, but this is purely on the basis of the slope and aspect of the cell it is looking at (and its eight neighbours).

In other words, a piece of terrain at the bottom of a deep gorge will be rendered purely on the basis of its slope and aspect, as if the sun could shine on it right through the surrounding mountains.

Hence, no, we do not need to worry that shadows will be shortened.

So…

All of these things considered, I come to the surprising (for me) conclusion that yes, you actually can often just hillshade a DEM that’s still projected in degrees, using the scale factor 111,120. It seems heretical, but when I can’t tell the difference in the results, I have to admit that for practical purposes it’s a good solution.

QGIS 3 and Vector map tiles

If you’ve ever tried printing a map from the Open Street Map tile service (available in QGIS through QuickMapServices plugin) you’ve experienced the annoyance of having all your labels come out tiny. In other words, you see this on the screen…

But then you get this out of your 300 dpi print layout….

The reason those labels and icons come out so small with raster tiles is that your screen display is something like 96 pixels/inch, but the print layout has to pack them down to 300 pixels/inch. The print layout also winds up asking for tiles at a higher zoom level. Text labels and icons, like everything else, are “baked” into the raster tiles at a reasonable size for display at 96 dpi. At 300 dpi, everything is much smaller.

In 2018 I described, at the end of a post, a way to work around this, but now there’s a better option. If you work instead with vector tiles, the text labels are not baked into the tile: they are rendered by QGIS. This means that if the label is supposed to be 2.5 mm high, it’ll be 2.5 mm high no matter what dpi you export at out of the print layout.

QGIS 3.14 and above can handle vector tile servers as data sources. This means the server sends tiles of vector data, which are then styled at the user’s end by QGIS.

How they work

Vector map tiles (VMTs) were not developed for GIS use. Their main use is for people viewing maps through browsers. You commonly have a map that’s been embedded in a web page that includes the javascript libraries for mapbox rendering, such as here.

But vector tiles offer the cartographer a lot. As raster tiles, OpenStreetMap data was ideal for areas where you just could not get decent data—if, say, you were mapping the streets of a city in Turkmenistan. But now, as vector tiles, this data allows you to control colours, line weight, font, and, of course the visibility of elements. Imagine your own version of OpenStreetMap, with pink highways, no railroads and emphasized parks, and you begin to get the idea.

There are multiple formats for tiles, but here we’ll just look at the Mapbox vector tiles format, since that is what QGIS can support.

Because vector map tiles are designed to be viewed over the internet, the most basic principle of VMTs is that you typically to have a tile server, which is firing off vector tiles of data, and a GL style file that understands what data layers are delivered in this server’s tiles, and tells the browser how to render that data. Often these are at different URLs.

The URLs of vector tile servers look like this:

The z, y and x in curly brackets get replaced by the zoom level, row and column, and the data comes as a PBF (Protobuf) file or a MVT (Mapbox vector tile) file.

The GL style file, on the other hand, is a text file in JavaScript Object Notation (JSON). It might begin like this, with information on the sources for tiles, sprites and glyphs (VMT lingo for markers and fonts):

{
  "version": 8,
  "name": "Positron",
  "metadata": {...  },
  "sources": {
    "openmaptiles": {
      "type": "vector",
      "url": "https://api.maptiler.com/tiles/v3/tiles.json?key={key}"
    }
  },
  "sprite": "https://openmaptiles.github.io/positron-gl-style/sprite",
  "glyphs": "https://api.maptiler.com/fonts/{fontstack}/{range}.pbf?key={key}",

and then it will continue with a list of the layers in the tiles, and how they should be styled.

  "layers": [
    {
      "id": "background",
      "type": "background",
      "paint": {"background-color": "rgb(242,243,240)"}
    },
    {
      "id": "park",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "park",
      "filter": ["==", "$type", "Polygon"],
      "layout": {"visibility": "visible"},
      "paint": {"fill-color": "rgb(230, 233, 229)"}
    },
    {
      "id": "water",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "water",
      "filter": [
        "all",
        ["==", "$type", "Polygon"],
        ["!=", "brunnel", "tunnel"]
      ],
      "layout": {"visibility": "visible"},
      "paint": {"fill-antialias": true, "fill-color": "rgb(194, 200, 202)"}
    },     ...

The styling information here is written in something called Mapbox GL style, and here are some GL Style file URLs:

Vector tiles are relatively new, and almost no one gives them away for free, but we do have a few heroic companies offering this. Most require API-keys or access-tokens, but these are frequently free.

The first question you might have is what you get if you ask for vector tiles from the tile server, but don’t combine them with a GL Style file. VMTs without a style file aren’t useless—but what you get is data where all points share a single style, as do all lines and polygons.

Data from the Nextgen vector tile server without styling

A simple VMT service brought in using QGIS3’s native abilities

Our best candidate at the moment for free OpenStreetMap (OSM) data delivered through vector tiles is the ESRI OSM vector tile layer. (Thank you ESRI!). The tile server is here:

https://basemaps.arcgis.com/arcgis/rest/services/OpenStreetMap_v2/VectorTileServer/tile/z/y/x.pbf

and the style file is here:

https://www.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/root.json?f=pjson

We can add this tile server to QGIS by

  1. opening the Browser panel
  2. Right-clicking the new Vector Tiles category and choosing New ARCGIS vector tile service connection
  3. Filling it out like this:

Notice that the Service URL ends with “VectorTileServer.” There’s no “tile/{z}/{y}/{x}.pbf.” QGIS will append that automatically before making a tile request.

When we add this layer to our map canvas, and zoom to Baku, Azerbaijan at 1:8,028 scale, we get this pleasing map that looks remarkably like OpenStreetMap:

However, we also get a pile of errors, under the banner “Vector tiles: Style could not be completely converted.” If you click on “Details,” some will be like this:

Referenced font Arial Unicode MS Regular is not available on system

You can repair these font errors by installing the missing font – or simply accepting the local font that has been substituted.

But there are also a pile of errors along the lines of

amenity area/bowling alley: Could not retrieve sprite 'amenity area/bowling alley'

I will explore how to solve these sprite problems below. First let’s talk about what you’ve got so far.

Styling the vector tile layers

If I open the styling panel on the right of the map I get something truly new.

Wow, that’s the longest list of rule-based styling I’ve ever seen!

Can I literally uncheck or change the style for every type of layer? Are these layers within a layer?

Yes and yes. Vector map tiles contain many layers. (The set of layers is called the schema of a tile service.) Out of each layer (e.g., “landcover”), features are selected using a filter, which is quite similar to a rule in rule-based styling. For example, in the image above there is a specific fill style for features in the landcover layer for which _symbol = 48,and the zoom is 9 or higher. It’s called landcover/park.

I can uncheck landcover/park and get a map without parks…

Or I can double-click on the sample fill swatch and change the parks to a pink/white gradient fill.

In fact I can apply any QGIS style to these polygons: shapeburst fill, line fill, marker fill… I don’t have full control—I can’t change the min and max zooms that define this style, or the filter rule—but within the predefined filters, the possibilities are mind-boggling.

Once I like the way in which I have tweaked the styling, I can save the style for this layer to a .qml file. (Open the layer’s Properties, find the Style button, and choose Save Style…) This will allow me to apply this custom styling to the same vector tile layer, but in another QGIS project.

My powers over the styling of the VMTs are clumsy—a GL Style file typically defines scores of styles, and if I want to change, say, all my water styles from blue to grey, I will have to do that to each one individually. But it’s still very powerful.

There are a few other things to point out. Note at the bottom of the styling panel you have a checkbox for Visible rules only. Visible rules only cuts out styles that are not displayed at the current zoom. (It does not remove styling for features not visible in the current map canvas.)

The current zoom is also displayed there. This is always shown as an integer, but with vector map tiles you are not restricted to integer zoom levels. There’s no “Zoom to Native Resolution (100%)” like there is when you right-click a raster tile layer. The zoom level shows as 14, but the map canvas is at 1:10,000, which is somewhere between 1:17,061 (the break point for vector tiles changing from zoom 13 to 14) and 1:8,530 (the break point for vector tiles changing from zoom 14 to 15). If you do the math, you can conclude we’re at something like zoom level 14.17. This might make you queasy! And, welcome to the world of vector tiles.

[Note that zoom levels for vector tiles are one more than the zoom levels for typical raster map tiles (e.g., Google Maps, OSM ). This is a consequence of MapBox vector tiles being 512 pixels on a side. As well, I have a suspicion that, in QGIS, the detail shown on ESRI OSM vector tiles at zoom level 14 is the same as shown on OSM raster tiles at zoom level 13. But I’m not quite sure why this is.]

Styles are implemented in order, from the top of the list down, and you can drag them to reorder.

Finally, it’s important to understand that the styling rules you see here in QGIS are not the GL Style document itself. They are a translation of that document. You cannot (at, least, as far as I know) translate back from QGIS styling to a GL Style file.

Labels

If you move to the labelling tab in the style panel, you see this:

So, you also have extensive control over the labelling. Double click on a label’s style swatch, and change the font, size, colour… all of the usual label parameters are available.

For example, if I would like street names to be blue (bold-italic) and subway stations to be labelled in red, I can just do it by altering three styles.

Finally, note that while you can’t open an attribute table on your vector map tiles layer, you can use the Information tool. You will get results for every single layer under the mouse click.

I can now print at 300 dpi from this, and I get just what I see on the screen!

Sprites

All that might have been enough for you, and now you’re off to try customizing OpenStreetMap yourself. However, if you’re interested in fixing those sprite errors we encountered above, continue reading.

When you add a vector tile service to QGIS you frequently get this:

When you check out the Details, many of the errors refer to “could not retrieve sprite.”

At the same time, in the QGIS’s Network Log window, you’ll see something like this:

2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite@2x.json - server replied: Not Found
2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite.json - server replied: Not Found

So what’s going on?

If I paste the URL for the GL Style file into my browser (I’m using Firefox), the result is that whole code for the styling gets delivered to me in an unformatted blob.

It’s a bit dense reading. But near the top I can pick out this piece:

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/../sprites/sprite"

Sprites, in vector map tile parlance, are markers. QGIS needs to use this URL look up the sprite image (a PNG file with all the sprites packed into it) and its partner sprite JSON file (which explains which piece of the image a given sprite is). But, for QGIS, the “styles/..” part of this particular URL is confusing. We might understand it as meaning descend into the styles sub-directory and then pop back up out of it, but QGIS doesn’t know this.

I can make a better address for the sprites by saving the GL Style file locally, editing it so the URL of the sprites is simpler, and then telling QGIS to now use my local version of the GL Style file.

I save the GL Style file locally by selecting all the text out of my browser, and copying and pasting into my text editor. Then I change the sprite URL by removing the “styles/..” part, so it now reads

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/sprites/sprite"

Once I’ve saved the file as ESRI_OSM_Style.json, I return to QGIS, where I edit the Vector Tiles connection I made for this server, and specify this new file as the “Style URL.” I could post it on a server somewhere and give a http:// URL to it, but, with Ubuntu, I can also just use the file:// protocol:

Sprites make quite a difference at high zoom levels.

No sprites, zoom level 17
With sprites, zoom level 17

A few final points

There are three different types of support for vector tiles in QGIS. The original one was the Vector Tiles Reader plugin, which received the vector tiles, converted them to geoJSON data, and made a different layer for each kind of features in the tiles. There is the new native support, which I have used here: it receives the tiles and applies filters, or rules, to style the many different kinds of features in them. Finally, there is the MapTiler plugin, which has a great deal of potential but requires a free account with MapTiler/OpenMapTiles, and this has an unfortunate data cap of 100,000 tile requests/month. (It seems like a lot, but while learning to use and style vector map tiles, I went through my 100,000 tiles in eleven days.)

The QGIS implementation of vector tiles permits you to add “New generic connection” if you have URLs for the server and the GL Style file. Yet outside of MapTiler and ESRI I have found precious few servers that also provide GL Style files. Here’s a quick survey:

  • The city of Wien (Vienna)—the data is limited to Austria. [Server URL] [GL Style URL]
  • Ordnance Survey UK—the data is limited to the U. K. and requires a free API key. [Server URL to which you have to append your free API key] [GL Style URL to which you have to append your free API key] Because of the API-key, QGIS does not yet seem to be able to make the request for the sprites correctly.
  • MapTiler/OpenMapTiles (OSM data)—beware the 100,000 tile limit. (Literally every time you zoom, pan or resize the screen you can trigger a tile request.) [Server URL to which you have to append your free access key] [GL Style URL for the “Basic” style]
  • Nextzen (OSM data)—lovely data, but I haven’t found a GL Style file. [Server URL to which you append your API key.]
  • ESRI generously offers another vector tile server based on more data sources than just OpenStreetMap. [Server URL] It comes with many different styles which you can find here. To obtain the GL Style file for one of these, you typically Open the style in Map Viewer and then in the upper right is a button called View style. Copy the link out if it and use it, as we did above, in the Style URL field. (Generally speaking, you’ll have to also do the same hack with the sprites as I showed above.)
  • Mapbox, if you can decode their mysterious URLs (mapbox://mapbox.mapbox-streets-v11), serves up vector tiles, and the public access key is free (without data limits), but I have not been able to find a GL Style file for its schema that QGIS can fully translate. [Server URL for Mapbox-streets-v8] [GL Style URL for MapBox-streets-v11]
  • Thunderforest has nice vector tile servers, such as their Transport layer, and their Outdoors layer, that you can have access to with a free API key—but I have not been able to find GL Style files for the schemas of these servers. [Server URL to which you have to append your free API key]

I understand that converting data to vector tiles, and keeping it up to date, costs money, so I don’t begrudge these companies charging fees or imposing limits. On the other hand, I really appreciate those which don’t!

As time passes there will no doubt be more vector map tile servers, and the way schemas and styling work will change. Nonetheless working with vector map tiles is addictive, and it will be exciting to watch this evolve.

Shaded Relief using Skymodels, courtesy of Raster Chunk Processing

A couple of weeks ago I watched some excellent presentations from the How To Do Map Stuff day organized by Daniel Huffman. One that I particularly enjoyed was Jake Adams’s talk on building shaded relief out of hillshades. Toward the end of his talk he brought in something called a skymodel.

In this post, I’ll explain what Skymodels are, and how to get Jake’s Raster Chunk Processing software running in a linux environment so you can use the Skymodel feature to make your own versions of this unique kind of shaded relief.

Skymodels were introduced in a paper by Pat Kennelly and James Stewart in 2014. The basic idea behind the skymodel is that, in the real world, we see terrain illuminated by light coming from all parts of the sky in various amounts, not just from the sun. The usual hillshading algorithm, on the other hand, calculates what’s illuminated and what’s in shadow based entirely on light coming from a single point.

An overcast sky is one example of a skymodel: in the overcast sky, light comes mostly from above, but not from any particular direction.  In James Stewart’s SkyLum software, this is represented by the redder, “hotter” dots in the dome of the model, while we have yellow and then cooler blue dots down toward the horizon representing less light coming in from those directions.

skymodel_overcast

Contrast that with this model for what SkyLum calls a Type 6 sky, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.”

Type 6 partly cloudy, no gradation towards zenith, slight brightening towards the sun

Or this, the Type 13, “CIE standard clear sky, polluted atmosphere.”

type 13 CIE cstandard lear sky, polluted atmosphere

Each of these will produce a different kind of shaded relief.

Baltoro_assortment
The Baltoro glacier, Pakistan. Left: overcast. Middle: type 6, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.” Right: type 13, “CIE standard clear sky, polluted atmosphere.”

In Blender you can position multiple light sources, but SkyLum writes a CSV file (which I will call the illumination file) that defines 250 different light sources and assigns a relative weight to each. This is new level of complexity in lighting.

180,45,0.000962426
175,37,0.000929834
193,47,0.00118523
187,38,0.000949383
167,47,0.00143157
181,54,0.00168069
169,41,0.00141631
...

Each line in the illumination file above gives the azimuth, elevation, and weight for one illumination source. The weights for all 250 points will add up to 1.

So how does one use this to generate shaded relief?

Well, Jake Adams (thank you, Jake!) has written a clever piece of code called Raster Chunk Processing (RCP, hereafter), which he presents in a three-part blog post. RCP divides up large DEMs into smaller tiles (“chunks”), each of which is processed separately. All of the results are then merged back together for a final result. The point of RCP is that it allows you to work with DEMs that ordinarily would max out your RAM and cause processing to grind to a halt.

This is similar to the way I have divided large DEMS into tiles for processing by Blender, but Jake’s RCP code allows you to use this “divide-and-conquer” strategy to do a whole host of things. One of them is to build a skymodel hillshade, given a DEM and an illumination file from SkyLum.

The RCP code is written in python, which is platform independent, so although Jake gives us instructions for installing it under Windows, we can get it running under linux with just a few changes.  In this case I did this on a system running Ubuntu 18 Bionic, and since I was using this machine for mapping I already had QGIS and GDAL installed.

To begin, head over to Jake’s Git repository, and download the RCP software using the Clone or Download button. Once unzipped, this will produce a directory called rcp-master. Within it you will see, among other things, the main program (raster_chunk_processing.py), a subdirectory called SkyLum (which contains the SkyLum program) and a couple of other python files that will be important to us, like settings.py and methods.py.

RCP is written to be run under python 3, so from a Ubuntu point of view these are the packages you will need to install:

  • python3-numba
  • python3-astropy
  • python3-gdal
  • python3-skimage
  • python3-numpy

(Note that python3-scimage replaces the module Jake calls scikit-image.)

Another thing that RCP needs is a working copy of mdenoise, the software that implements Sun’s algorithm for denoising topographic data. Or rather, RCP needs to at least think it has a copy. So you have a choice: if you want to be able use RCP to denoise DEMs, you should compile yourself a copy of the mdenoise binary; there are instructions here, and it’s not too painful. If not, just use sudo to place an empty file called mdenoise in /bin.

Then use a text editor on the file settings.py to alter its single line about where the mdenoise executable is.

MDENOISE_PATH = '/bin/mdenoise'

The last piece you have to put in place is a way to run SkyLum. SkyLum only comes as a Windows binary, SkyLum.exe, so to run this you need to have wine installed. The good news is that SkyLum runs quite well under wine.

Right-click SkyLum.exe, and choose Open With… >Other Application. In the list of application choose, “winebrowser.”

winebrowser

SkyLum should open right up and show you a piece of terrain with a sky dome over it.

fresh skylum

The complete instructions for SkyLum are in the README file included with it, but I will give a summary here of what I find useful:

  • By default when you open SkyLum the sun (sun position is shown at the bottom) is at 45° elevation and azimuth 180° (due South). 0° is north, and azimuths increase clockwise, as is standard
  • Hit ? for a help screen
  • Move the sun with the arrow keys.
  • Rt-click and choose an illumination model. See Kennely and Stewart’s paper for more on these.
  • Hit p to have SkyLum.exe calculate points after you have positioned the sun and chosen a skymodel
  • Hit o to have SkyLum write a sky model file. It’s conventional to give it a .csv extension.
  • Use a text editor to delete the header lines. All you want left are the comma-separated lines with the azimuth, elevation and weight, as shown above.

Now, how to run RCP?

Let’s assume you have a DEM  called myDEM.tif, and an illumination file you made with SkyLum called illum.csv. You’ve already deleted the header lines from illum.csv. You want to divide your DEM into 1000x1000px chunks, with an overlap of 200 px. We’ll also assume you have 4 processors and you want to use 3 of them for this operation.

Drop both myDEM.tif  and illum.csv in the rcp_master directory. Open a terminal there.

The general form of the RCP command line is:

python3 raster_chunk_processing.py -m [method] {general options} {method-specific options} input_file output_file

Notice that the first element in the command line is python3, not just python.

For my skymodel in this case the command line will be…

 python3 raster_chunk_processing.py -m skymodel -s 1000 -o 200 -p 3 -l illum.csv --verbose myDEM.tif myDEM_skymodel.tif

What you’ll see in the terminal is a bunch of information about the job to be done, and then you’ll see RCP submitting the sub-jobs (the chunks) to be skymodelled. You can walk away, or watch, fascinated, as the chucks get worked on…

Preparing output file myDEM_skymodel.tif...
Output dimensions: 4046 rows by 4515 columns.
Output data type: <class 'numpy.float32'>
Output size: 69.7 MiB
Output NoData Value: -32767.0

Processing chunks...
Tile 0-0: 1 of 25 (4.000%) started at 0:00:00.701178 Indices: [0:1200, 0:1200] PID: 8777
Tile 0-1: 2 of 25 (8.000%) started at 0:00:00.701253 Indices: [0:1200, 800:2200] PID: 8778
Tile 0-2: 3 of 25 (12.000%) started at 0:00:00.701284 Indices: [0:1200, 1800:3200] PID: 8779

All of the switches and parameters are explained quite well in Jake’s post. My additional notes are:

  • the input file is generally a geotiff, but it can also be a TIF/TFW pair. It has to have a NoData value defined or you will get an error No NoData value set in input DEM.
  • RCP will exit with an error if the output file already exists
  • By default RCP applies a vertical exaggeration of 5X to the DEM, because this was what Kennely and Stewart did in their paper. However, you can change this if you would prefer a different vertical exaggeration. Open methods.py in your text editor and go to line 272. This line originally says  in_array *= 5. You can change that 5 to whatever number you would prefer. (No recompile necessary.)
  • the overlap (-o) parameter is based on how far you think shadows may stretch. The shadowing algorithm checks outward for 600 pixels to see if a given point is shadowed by anything. For this reason, there is no point in making overlap larger then 600.
  • the chunk size parameter (-s) is based on how much RAM each process requires while running. You can experiment with this and watch on the system monitor to see how close you are to spilling over into swap.
  • When shadows are being calculated in the skymodel,

The output file can be dragged into QGIS, where I often find I want to increase brightness.

adjusdting_brightness
Left: hillshade image fresh out of RCP. Right: with brightness +100

Another idea that seems to produce some nice results is to combine two or more skymodels with different transparencies and styling. You might also want to check out a gallery of the results of applying all of the different skymodels in SkyLum to the same piece of terrain.

Now that you have RCP running, of course you’ll need to try all of the different skymodels on your favourite DEM to see which you like best. But it’s also worth checking out the other kinds of processing RCP can do on DEMs. It is a very fast conventional hillshader (method hillshade). It runs a nice, quick gaussian blur (method blur_gauss), much faster than the SAGA Gaussian filter module in QGIS. And I haven’t even tried the Clahe and TPI processing yet!

 

 

 

How do we understand the size of Syria?

Syrian_comparison_inset_2004_CIA

This inset appears on a CIA map of Syria from 2004. We can assume it’s meant to give the map reader a sense of the size of Syria by comparing it to a region that he or she is familiar with.

I’m fascinated, first, by the assumption that the decision-maker reading the map is located in the mid-Atlantic states, (Washington, presumably). This kind of regionality runs throughout American politics, a geography of where important people live, and where they don’t, that one carries in the memory and consults without even knowing it. 

This is an attempt to make the map *personal* but I almost think it should be captioned, “This might be helpful to you if you happen to live in the northeast.”

The second thing I’m fascinated by is, assuming we have to stick to the northeast USA, would it have been a smarter decision  to have aligned Damascus with Washington, DC? These are the two national capitals. Such a juxtaposition would, in the context of the Syrian Civil War, allow the President to imagine New York no longer doing his bidding, and Providence functioning like an independent state.

CIA inset Washington and Damascus aligned

Perhaps this juxtaposition puts too much of Syria out to sea, but it does put New York City more or less where Raqqa, the ISIS capital, was.

Of course one must always be careful when constructing these comparison maps. In GIS software a polygon (Syria) whose coordinates are in degrees, dragged to another latitude, will be the wrong size. The safer method, which I have used here, is to make separate maps at the same scale of both Syria and the northeast US, and then juxtapose them in a photo editing program (the GIMP, in my case).  Syria is about 780 km from end to end on its long axis, about the distance from Boston to Richmond, Virginia.

This brings up the question of what it means to understand distance. In the original map, what does the distance from, say, Philadelphia to the eastern tip of Tennessee, mean to people living in Washington? I suspect they rarely go very far to the southwest from the city, and when they do they experience slow, twisting roads that have difficulty passing through the Appalachian mountains. So if the Washington-dweller says “It would take me six hours to drive to the tip of Tennessee!” is there anything meaningful at all there for his understanding of Syria?

I kind of prefer this juxtaposition, both for the similarity of desert landscape, and the sense of distance.

CIA inset LA and Damascus aligned

Los Angeles stands in for Damascus here, and Las Vegas finds itself somewhere up on the Euphrates. (“Las Vegas on the Euphrates” is not yet, but may someday be, the tourism slogan of Deir Ez-Zur or Raqqa.) If you are familiar with the distances and deserts of the American southwest, it is remarkable how much smaller Syria looks when shown like this.

Of course, as Appalachian Trail hikers know, there happens to be a town in Virginia called Damascus. It’s actually right there, just across the border from that eastern tip of Tennessee. So perhaps the best juxtaposition of all aligns Damascus, Virginia, with Damascus, Syria.

CIA inset Damascus and Damascus aligned