A friend asked me to make them some contour maps of Kullberg, given I clearly had all the data necessary, and given it was not something I'd tried before, it seemed like a fun thing to learn. In the end it was quite easy, but with two slight gotchyas I had to work around.

The first thing you need is your elevation map, in this case a tile on the area of interest from Lantmäteriet 1m Markhöjdmodell:

In theory all you need to do to generate a contour map from that is use a single GDAL command, gdal_contour:

$ gdal_contour mhm-70_6_707_60_2500.tif contours.geojson -i 1

The -i parameter is a distance in map units between contours. In this map each pixel is a floating point value of height over sea level in metres, and so I'm going to generate a contour at 1m steps. This will give you this map:

Success, right? Well, it's fine, but it has issues. To help discuss the first one, let's quickly look at the local land cover map, taken from Nationella Marktäckedata:

In the land cover data you can see the peninsular and the island, and if you look closely at the contour map you can see that the river effectively floods over the island, and the peninsular is also disconnected.

The reason for this is that the water level in the map is 235.08 metres, and so what you're seeing as the first contour for the island is the line at 236 metres, which is considerably higher than the water level, so we have lost the outline of the island.

You can give gdal_contour a list of steps, so I tried to pass it 235.08 or slightly lower as the starting point, but it seems it always rounds to the nearest integer unit of height. So in the end I fixed this by editing the elevation map, and made all pixels that were 235.1 and lower down to 234.99 metres. Basically I lowered the water by 9cm :) That said, that's well within the bounds for this area, as this river is actually a reservoir and the water levels do change quite a bit, and I suspect the elevation map is towards the higher end of things. You can see in this screen grab I took today from Den Stora Älgvandringen the difference in water level here as marked on this rock which is in that narrow channel between the island and the peninsular: the water level the last few days has been really quite low (and a tip of the hat to Snuttenliten in the moose discord channel for spotting the handy water level marker!).

It's something you have to watch out for in working with tiled elevation data: the tile to the left of this one has the water level slightly lower than the one we're looking at today, at 234.81, as most likely it was recorded in a different pass. That's not a big difference, but there's potential for a reservoir like this to change level quite quickly, and so throw off attempts to match tiles together. I also know that the river is the lowest point in this particular tile, so I'm only adjusting the water level and nothing else. If that wasn't the case then I'd need to use the water level mask I created before.

To do the height adjustment I just used my own Yirgacheffe Python library, as it makes jobs like this quite succinct:

import yirgacheffe as yg
r = yg.read_raster("mhm-70_6_707_60_2500.tif")
l = yg.where(r <= 235.08, 234.98, r)
l.to_geotiff("lowered.tif")

That second to last line effectively says "in each pixel that is lower than or equal to 235.08, make it 234.98, and leave everything else as it was". I'm tooting my horn I realise, but I really love that Yirgacheffe exists at times like this.

Now if I re-run the command with the adjusted water height, we get a much clearer definition of the peninsular and the island:

A related note on the -i option: if you pick 2m say as your distance between lines, then it will always pick elevations that are multiples of 2, and so I'll get contours at 234, 236, and so forth. Given that this map has its effective water level at 235 metres, that means we've back to having a bad outline, as the first step is back up at 236. To work around this I had to ditch the -i option and instead use the -fl (fixed list) option, and then specify 235 237 239... all the way up to 311 for this area. Perhaps this is the point I should have scripted something!

---

The second issue with the contour map is that it's quite noisy, as every little undulation is captured, whereas contour maps normally have some level of smoothing to them. Thankfully we can use mapshaper for this. I'm going to use the command line tools here via npx, but there is a browser based GUI version too. I'm actually already using mapshaper to convert convert from GeoJSON to SVG file also, and it turned out there is just an extra option I can pass it to reduce the noise in the map:

$ npx mapshaper contours.geojson -simplify 2% contours_simplified.geojson

And from that we get a map this is closer to what one thinks of when thinking of contours:

I found I had to be very aggressive with the simplification - even at 10% detail, there was still quite a bit of noise in there. You still have some quite steep bits on the island at 1m, so in the end I also made a 2m step version:

You can now all print that out and colour it in :)

---

There is one other note in all this that is worth noting if you're doing this for actual GIS purposes, rather than making an SVG to display. Initially I was using QGIS to review the GeoJSONs I was making, rather than pushing them through to SVG files; QGIS is the standard tool used in geospatial work for this kind of thing. What I found was that although the output of gdal_contour displayed fine in QGIS, if I ran mapshaper's simplify and generated another GeoJSON rather than an SVG file, QGIS wouldn't display it. But the file was 12MB of data, so I knew there should be something visible.

After a bit of hunting I realised it was a problem with a difference in how GDAL and mapshaper interpret the GeoJSON standard. In earlier GeoJSON you could specify a Coordinate Reference System (CRS), and the data in the GeoJSON could use that CRS. However, when the GeoJSON standard was updated in 2016 the CRS field was dropped and all data stored in a GeoJSON needs now to be in WGS 84.

Some quick context to try make what happens next make more sense. WGS 84 means you use latitude and longitude in the GeoJSON file, and when that gets projected onto a flat rectangle like a map, or an SVG file, you treat the latitude values as your y value and the longitude values as you x value. But the problem with this in general is that as your latitude changes, so does the circumference of the earth, so just converting longitude to an x value when you plot the map gets increasingly distorted as you move to the poles. I've discussed problems with WGS 84 before, so I'll spare you the deep dive, other than to note that Sweden is quite near the poles, and so for its officially published geospatial data it doesn't use WGS 84 (nor does the UK for that matter). If you want to see why, you can compare the above maps with this one where the data is in WGS 84:

It looks squished, but that's because my blog software has tried to keep all the images on this page the same width, it should be the same height as the previous one, but much much wider, due to the aforementioned distortions. The point is, what is a relatively square area if you are looking down over it locally, becomes a very different shape in the WGS 84 Mercator projection.

Anyway, back to why QGIS didn't display anything. Because GDAL knew that the GeoTIFF data was in EPSG:3006, Sweden's local CRS, it made an old style GeoJSON with that CRS embedded, and all the point values in that CRS. EPSG:3006 actually uses metres from a fixed origin, so the points in that are like 500000 across and 7000000 up. Mapshaper however assumed there were very large latitude and longitude values, and made me a GeoJSON with simplified versions of those numbers in, but claimed the output was WGS 84. So QGIS presumably just threw it's hands up in the air as these values don't make sense in WGS 84, and didn't display anything (annoyingly it didn't seem to complain either).

It took me a while to figure this out, and in the end I had to use ogr2ogr, another tool in the GDAL suite to let me translate my GeoJSON files between the local CRS and WGS 84 and then back again afterwards.

$ ogr2ogr -t_srs EPSG:4326 coverted_to_wgs84.geojson source.geojson

I think the more I do work in GIS the more I think it's fine to have latitude and longitude as a storage format, but never as a rendering or rasterising format.