Before I became a cartographer, I spent ten years teaching wilderness travel at the National Outdoor Leadership School. I picked up GIS skills from a diamond exploration outfit in Canada’s Northwest Territories As a cartographer I split his time between topographic trail maps and doing historical mapping.
This is a procedure that came up in a discussion with a friend, and I think it is tricky enough to be worth recording here.
Specifically we are using QGIS 3 to georeference a 1941 map of the Odessa, Ukraine, area, one of the 1:1,000,000 International Map of the World Series.
This map is bounded by latitudes 46° and 51° north, and longitudes 27° and 36° east.
We want as little loss of image quality as possible, therefore we want to avoid warping (re-projection). If warping the map were not a concern, we could georeference it in a geographic projection (e.g., EPSG 4326 or 4267) with a few ground control points (GCPs) in degrees, using the intersections of latitude and longitude lines. But in this case we need to georeference it in the original Lambert projection, as it was printed. The transformation will be “linear” and in fact only a world file will be written. The world file will enable QGIS to read in the map image without warping it.
This will not be possible unless we can figure out what the original projection was. Fortunately at the bottom of the map the person who did the scanning has left the statement of projection.
A Lambert conic projection usually relies on four parameters. There are the two parallels of latitude at which the cone touches the globe: these are the two numbers listed here: 36° north and 52° 48′ (52.8°) north. They will be called lat_1 and lat_2 in the projection definition.
Then there are the coordinates of the origin point of the projection: lon_0 and lat_0. The meridian that runs straight down through the centre of the map is clearly the central meridian of the projection, because it runs perfectly vertically, the only longitude line that does so on the map. The centre of the map falls halfway between longitudes 31 and 32, so lon_0 is 31.5°.
Lat_0 is a little harder to figure out. However, in my experience it doesn’t really matter what you choose for lat_0. I chose 40° north.
The first step then is to create a new CRS (coordinate reference system) in QGIS for this Lambert Conic projection. We go Settings>Custom Projections, and click the “+” button to add a new CRS.
Plugging in our parameters we determined above into the normal Lambert Conic definition, we get this:
QGIS will assign an EPSG number for your new CRS. In this case I got 100030, but they are always different (and greater than 100000).
The next step is for me to put my main map window in QGIS into the new projection, and open an array of latitude and longitude lines, such as the Natural Earth 1:10 million scale one-degree graticule layer. And I turn on labels for these lines so I can see what degrees I am looking at.
The reason I do this bears directly on the central technique we are going to use here. Because the original map is in Lambert, and I want to register it in Lambert, I will have to enter Lambert coordinates for each of the GCPs. But looking at the map I only see degrees: I don’t see Lambert coordinates. Fortunately QGIS can tell me what the Lambert coordinates are for each grid intersection, as long as I am displaying my grid in the same CRS as I want to use to register the map. You can witness this by zooming in on and hovering your mouse over a grid intersection and looking at the Coordinate text box in the bottom margin.
There are the Lambert coordinates—at least for this specific Lambert Conic projection we are using—for 46° north, 26° east: -421460 metres east, 674061 metres north.
We don’t need to write these down though. We will get them assigned to the map image in a more elegant way.
Now it’s time to open the georeferencer (Raster>Georeferencer) and bring in the map image. Immediately you will be asked which CRS you want to georeference this map in. Choose the custom CRS you just made (in my case, 100030).
Before you place GCPs, go to the settings of the georeferencer and set up your transformation parameters to ensure no warp will result.
You want transformation type to be Linear. The resampling method is not important because no resampling should occur. The target SRS is your custom CRS. And you have checked the box called Create world file only (linear transform). (I also like to check the Load in QGIS when done box.)
Note that in the georeferencer the bottom line now looks like this:
I’m going to place the first GCP in the lower left corner (the southwest corner), which is 46° north, 27° east. Before I do this, I go into the main map window and zoom in on just that intersection, to a fairly large scale, say 1:10,000.
Now in the georeferencer I place my GCP point on 46° north, 27° east, and QGIS asks me for the X and Y of that point. Or, I can click From map canvas.
Once I click From map canvas, the georeferencer is hidden, the cursor becomes a cross, and I am invited to pick that same point on the main map canvas. I carefully click right on the intersection of 46° north and 27° east, and immediately the Lambert coordinates of that point as filled in for me in the dialogue box:
I click OK, and I’m ready to do my next GCP. Remember, the first thing I will do is pan the main map to the coordinate intersection where I’m about to add this second GCP.
A linear transformation never requires more than three GCPs, but I like to do four so I get an estimate of my error. I do the four corners of the map.
At this point I can see, from the GCP table at the bottom of the georeferencer, how my error is.
The residual looks like 3 to 5 pixels, and the mean error is 6.5 pixels, quite acceptable in this image which is 4700 x 5700 pixels.
Now I hit the Start Georeferencing button (green triangle “play” button) and a world file is written. Because I’ve checked Load in QGIS when done, I immediately get asked for the CRS of the new georeferenced map, and again I select the new custom CRS I created for this map.
The map appears in the main map window, and I drag it under the graticule layer, to see that it is properly georeferenced.
For example, consider this pattern in a window in the external courtyard wall at Baku’s Taza Pir mosque.
This beautiful pattern, with its eight-pointed stars set within octagons, turns up on plate 67 in Jules Bourgoin’s 1879 Les Éléments de L’Art Arabe (which you can download from archive.org).
It’s wallpaper group is the fairly common *442 (p4m) and it is generated by tessellating a square cell.
Construction of this pattern is straightforward. The eight-pointed star in the centre is inscribed in a circle whose radius is one quarter the side of the square. The vertices of the octagon are found by extending the sides of the star. The rest of the construction lines are extensions of the octagon sides, and lines connecting star dimples that are three apart.
But one enters the Taza Pir compound via a stairway from the street. The panels in the stairwell are related, but subtly different from the window design!
What did they do here? There is the same eight-pointed star in the centre, and the same enclosing octagon, but in this case they’ve trimmed back, to the borders of the octagon, the square that one repeats.
As a result, the square tile borders stand out strongly as lines, and around each point where four tiles meet we get a big diamond holding four small diamonds.
(This also belongs to the *442 wallpaper group.)
Now, in the Old City, I came across a piece of octagon-based decoration that illustrates what happens if one doesn’t follow best practices, as explained by Eric Broug. This pattern involves starting with the same pattern as in the Taza Pir windows (above; a.k.a. Bourgoin’s Plate 67), but then repeating a somewhat random subset of it. In other words, incorrect tessellation.
It would appear that the manufacturer of these pre-cast concrete blocks selected a piece out of the overall pattern that was not the all-important basic square, but rather a rectangle.
Hence each of the concrete blocks looks like this:
When you put them together, lines match up, but the effect of the original design is lost.
The wallpaper group of this pattern would be*2222 (pmm).
Elsewhere in the Old City, there were pre-cast patterns that did tessellate pleasingly, again with octagons.
Eric Broug, from the School of Islamic Geometric Design, writes that sometimes while travelling he sees a piece of contemporary Islamic geometric design and recognizes it as, well, let’s say less than best practice. (He sometimes posts images of this sort of thing on Instagram under the hashtag #cpigd, which stands for Common Problems in Islamic Geometric Design.)
What sort of mistakes are they? He explains on his page on best practices, but to me the most common one is incorrect tessellation, where a block of pattern is repeated in ways that cause lines to abruptly stop instead of continuing on.
I thought it would be interesting to take a look at the designs I found in Azerbaijan, in the cases where they were identifiably Islamic, and ask the same question: are they examples of good, traditional design.
So, bearing that in mind, let’s look at a grill window I found in the Divankhana of the Palace of the Shirvanshahs in Baku, Azerbaijan.
The Palace of the Shirvanshahs is the premier piece of historical architecture in the Old City of Baku, or, as it’s called in Azerbaijani, İçərişəhər or Icheri Sheher. The original buildings have been deduced to date from the 15th century, but most of the palace was heavily renovated/restored in the 20th century so it’s not immediately clear whether the details one sees are original or the work of a restorer.
The Divankhana (which is also variously called the Divan-Khane or Divanhane) is a structure in its own courtyard just off the outer courtyard of the palace. It holds a pleasing octagonal pavilion. whose original function is unknown (there are many theories). The pavilion is domed and consequently two stories in height, so it stands above the courtyard wall and can be seen from the palace courtyard. In one corner of the Divankhana, there is a staircase leading up to a locked door on the upper storey of the pavilion.
This window is at the top of that staircase, and looks out into the palace’s outer courtyard. Here it is seen from inside.
What I thought when I saw this is There’s no way this can be a best practices design. There were so many wacky elements that I had never seen in an Islamic geometric design before. For one, I couldn’t find a single axis of reflection in it, anywhere. For another, it contained a number of strange, 3-way intersections.
Is it a bad design, perhaps a modern artist not working within traditional lines, or could this be authentic traditional design?
Let’s look at how the pattern works.
The whole pattern begins with a pair of adjacent large octagons.
These fill the window, as shown.
Centred on the vertices of each octagon are eight smaller octagons, sized so that when they overlap they bisect each other’s sides.
As you might design it on paper
Clipped to the window opening
These circles of smaller octagons define an empty space in the centre of each of the larger octagons, a space which is an eight-pointed star.
So far, so good, and very symmetrical and, in fact, infinitely tile-able. The cleverness comes with what they did inside each 8-pointed star.
They divided this space with a four-armed pattern which has rotational symmetry but no mirror symmetry.
They ran it in opposite directions in the top and bottom halves of the window. So, looking from inside the window, there is a counter-clockwise star in the top half and a clockwise star in the bottom half.
This is fairly mind-boggling for traditional Islamic design — I think. (I’m no expert.) The little pattern of four “hammerhead” shapes that circles within each 8-pointed star looks more M. C. Escher than standard geometric design.
Mathematically, we might ask: does this pattern at least pass the test of being able to be continued in all directions?
Well, yes. One would simply add more big octagons above, below and to the sides, and then add the smaller octagons, etc. This could go on forever.
Each big octagon hosts an eight-pointed star in its centre. If you alternate big octagons that host clockwise stars with big octagons that host counter-clockwise stars, following a checkerboard-like pattern, the centre of each star would be a four-fold centre of rotation. The corners where four octagon tiles come together are two-fold centres of rotation. And there are axes of reflection along the lines where adjacent big octagons touch.
Axes of reflection in blue; four-fold centres of rotation as green squares; two-fold centres of rotation as red diamonds.
The fundamental tile, from which the entire pattern can be generated through reflection and rotation, is triangular.
And, it gets even a little more complicated.
At the Palace of the Shirvanshahs I never shot a picture of the outside of the window, so when I got home I found that, while there’s no Google Street view in Azerbaijan, someone had conveniently taken a photosphere in the main courtyard of the palace four months before I was there.
Here, in the photosphere, is the facade of the Divankhana as seen from the outer courtyard. The Divankhana is the building with the whitish dome.
Last month I went to Tbilisi, the capital of the Republic of Georgia. Before I went, the only map I had of the city was this inset on the 2007 (third) edition of ITMB’s Georgia. It covers the downtown area that most visitors are interested in: on the right bank of the Mtkvari river from the Metekhis Khidi (Metekhi Bridge) at the downstream end to Respublikis Moedani (Republic Square) at the upstream end.
There are many things I like about this map. It has a basic visual hierarchy: Rustaveli avenue is given extra width to show its importance, and the high-traffic roads are in yellow. Also, they’ve made an effort to use Georgian words (moedani, khidi, gamziri, kuča) for common features (square, bridge, avenue, street) — even though the transcription scheme is not what is recommended today (kuča would be written kucha). This is surprisingly rare: most maps, in fact every other map I encountered on this trip, translated those terms. But if you think about it, leaving them in Georgian, with a handy little table of geographical equivalents in the legend, has a distinct advantage. It immediately gets you started with a few words of vocabulary.
Now, this map is not detailed enough for finding your way around the old town, which is the tangle of curved streets near the Metekhis Khidi. For example, if you look near “Blue Bath,” you see a lot of generalized, unnamed streets.
Before I went then, I made myself a paper map from Open Street Map. I cruised around in OSM at zoom level 17, stamping out PDFs. I printed these out and then cut-and-taped them together. Besides being fun, it was an excellent way to familiarize myself with the city’s geography, and I highly recommend this exercise.
(Ignore the yellow highlighter on the map. That was just my way of keeping track of where I’d been.)
One thing that’s good about OSM is that it is generally in the local language. In the case of Tbilisi, the streets are mostly labelled in Georgian, and this gives you a reason to learn the Georgian alphabet – which in turn makes reading street signs easier. (If you would like the Georgian alphabet at hand while reading this post, the Wikipedia page on Georgian scripts is quite helpful. This is the mkhedruli script.)
But if you look closely there are also a surprising number of features labelled in other languages: there is some English, like “SIXT rent-a-car” or “The Statue of King Vakhtang Gorgasali.” There’s also some Russian (Храм Метехи). I guess this is the nature of volunteered geographic information. On the positive side, from this kind of thing you can deduce some things about who visits Tbilisi, and who’s interested in mapping it. Actually, in my two-week experience, Russian-speaking tourists in Georgia outnumbered those speaking English.
According to the OSM wiki, zoom level 17 is roughly equivalent to 1:4,000 scale, making this the most detailed paper map you are likely to find. I found it very accurate, even the pedestrian streets (blue tinged) and the incredibly valuable pedestrian underpasses. It’s also quite up-to-date. For example, where ITMB had Leselidzis Kuča in 2007, the street has now been renamed to K’ot’e Apkhazis Kucha.
However, OSM can be quirky. The major Tavisuplebis Moedani (თავისუფლების მოედანი), or “Freedom Square,” is not labelled at all, although the monument in the middle of it (Tavisuphlebis Monumenti/თავისუფლების მონუმენტი) was labelled on my screen as I was looking at OSM. And that label disappeared in the export to PDF. A bit weird.
The area around Tavisuphlebis Moedani/Freedom Square in OSM, zoom level 17, as seen in my browser (left, a screenshot) and in the PDF generated with Share>Image (right)
Once I had arrived in Tbilisi. I began hoovering up maps. The first one I got, since I was staying at the Envoy Hostel, was the Envoy hostel’s branded map. (“Envoy Hostel Tbilisi City Map”) I imagine this map turns up in many Tbilisi locations under different branding. It’s dated 2017 and a tiny note in the corner says it is made by “Saniday.”
There are lots of good things to say about this free map. The streets are clean and mostly labelled. Points of interest are flagged, although many of these are just the establishments for eating, drinking and shopping that have sponsored the map. Most notably, they’ve picked up a couple of widenings in the streets which catch the eye as landmarks, and have named them, , like “Jerusalem Square” and “Meidan Square.”
I proceeded to Prospero’s Books on Rustaveli Avenue (a delightful bookstore) and bought three maps. I bought the Tbilisi City Miniguide and Map (undated with no publisher mentioned)…
the GeoLand Tbilisi City Map 2015…
and Tbilisi: A Walk Around Old City, by Idea Design Group (undated)….
Let’s start with the Tbilisi City Miniguide and Map. If you look at the area between the Metekhi Bridge and the Peace (Mshvidobis) Bridge it’s apparent that this piece of cartography is a disaster.
On a recent trip to Thessaloniki I acquired quite a few maps of the city. It turns out that finding a really good one is not trivial. Here’s the one I kept going back to:
This is the Fraport map which is available for free at the airport. (I found a display of them in the baggage claim area.) I’ve written in the bus routes on the major avenues myself. The cover looks like this:
But I was not able to get a hold of this map before arriving. The map that I printed out at home was this one:
Some of the other maps you can acquire in Thessaloniki itself are worth looking at. The next three I picked up at the tourist information point at the bottom of Aristolelous. Here is their Thessaloniki City Map:
As you can see, not all streets are labelled. Its cover looks like this:
And there is the even more spare Thessaloniki Museums’ Map:
and its cover…
And the Thessaloniki Monuments Map:
with its cover (this is the Greek edition):
I have to say the Museums’ and Monuments Maps are so skeletal as to be useless for navigation in the city. Also, there’s some inconsistency. The Museums’ map and the Monuments map both identify the street that leads south from Antigonidon Square to Egnatia as “Παπαζωλη/Papazoli,” whereas the other maps agree that it is called “Antigonidon.”
A map that I bought in a bookstore was this “Best of” map:
It’s pretty good, but the lamination actually caused me not to use it. Too hard to draw on.
Last but not least, that all-important bus route map. You can get this at the tourist info point as well.
This is a tiny map, and not quite as detailed as you might like, but oh-so-valuable. You don’t use it for finding your way around on foot, only for figuring out where bus routes go and which one you want.
The 315 black-and-white hillshade from openmaps.gov.bc.ca. The black-and-white hillshade with the sun at an azimuth of 315° is typically the most useful of the layers on this WMS server.
It’s got some limitations:
there are only two sun azimuths to select from: 315° (northwest) and 225° (southwest). More flexibility would be good, because each landscape seems to have a different ideal azimuth to bring out the landforms that you want to bring out. More about this down below.
there’s no height exaggeration possible
occasionally there are odd artifacts, like straight lines running across the hillshade
On the other hand, you can adjust its brightness and contrast, and use the Multiply blend mode, which means you can do some nice things with it.
When you set the blend mode of a hillshade to MULTIPLY, it shows through upper layers.
But, if you want to adjust the azimuth or exaggerate height, you’ll need to find a DEM and make your own hillshade. It’s well-established (though not well-explained) that the human eye needs to see light from above, preferably from above and to the left. If your map is, say, south-up, you will need a hillshade where the light comes from the southeast (which will be in the upper left).
Now south is at the top (map is rotated 180°) and the shaded relief appears flat or inverted. For this map you would need a hillshade where the sun’s azimuth is southeast, or 135°.
Here the map is still rotated 180°, but the hillshade has been manufactured with an azimuth of 135°. The mountains look like mountains again.
A 3x vertical exaggeration of the terrain emphasizes the detail in the valley bottoms.
The overall process
Let’s talk about the three principles.
The hillshading algorithms require a DEM in a metric projection. That means that DEMs projected in degrees won’t work: you have to re-project them first. Unfortunately, just about all DEMS come projected in degrees.
The scale of your final map determines what sort of cell size you want in your re-projected DEM. A DEM with 10m cells is far too detailed for a map at 1:500,000, and the file would be enormous. On the other hand, a cell size of 500m would make a very coarse hillshade at 1:500,000. As a general guideline, divide the denominator of your scale by 5000 to get roughly what cell size you want. So if your map is 1:100,000, you’d be looking for a DEM with (roughly) 20m cells.
It is an option to style any DEM as “Hillshade” (other options are singleband grey, multiband colour, paletted, etc.) but the hillshading algorithm in the toolbox (Processing>Toolbox) produces a better result.
Acquiring DEM data
The first thing you need to do is pick your resolution. Because DEMs usually come projected in 4326, DEM resolution is typically expressed in arc-seconds, or seconds of latitude. Because degrees of latitude in BC are bigger than degrees of longitude, these cells are not square. They are upright rectangles.
What you want to know is how these non-square cells measured in arc-seconds will convert to square cells measured in metres. Here’s a handy chart.
You’ll notice that I favour the DEMs created from the Shuttle Radar Topography Mission (SRTM). These have a few advantages over the DEMs produced by Natural Resources Canada or the provincial mapping agencies:
they are usually more recent (dating to about 2000)
they represent actual measurements, as opposed to a grid generated from contour lines
they go right across provincial and international boundaries
From here, I’ll demonstrate how this works with an actual example. In this case I want shaded relief for a series of maps that I’m making. All the maps are in the same area, and they range in scale from 1:45,000 to 1:130,000. Looking at the handy chart above, this scale suggests I’m going to want to use SRTM 1 DEM.
Displaying an area in a print composer with grid of lines in CRS 4326 is one way to figure out which tiles you’ll need.
So I go to Recordlist, enter the bounding box and select SRTM1.
Downloading DEM tiles from Recordlist
Once they are downloaded and unzipped, I’ll read them into QGIS to confirm that they cover the right area.
Merge and Clip
You want to merge all of the individual DEMs into one, using Raster>Miscellaneous>Merge. (Incidentally, they do not have to be read into QGIS to do this.)
OpenStreetMap data is often the best large-scale data you can find in regions of the world where governments are not yet distributing free, open geodata.
There are several ways to download OSM data (QGIS plugin, Geofabrik, direct from OSM), but if the format you prefer is shapefile, and you have a specific area you’re interested in, the old Weogeo service may be the easiest and the best way.
Weogeo seems to have been bought or otherwise taken over by the for-profit company Trimble. I appreciate Trimble hosting this service and keeping it alive, but since Trimble integrated it into their pre-existing sales system, you have to go through a number of odd steps to get your free data. Since it’s free, open data, you hope for one of those good-feeling interactions that suggest the sharing-without-strings that OSM represents. Instead be prepared here for a less comfortable, more corporate, interaction where they’ll want to you to go to their “data marketplace,” make an account, “add to cart,” and so on. But it works really well.
Two important notes up front:
If you use the Firefox plug-in “Privacy Badger,” it breaks this site. You have to disable it for trimbledata.com.
You won’t get your data immediately. Depending on the order size you may have to wait up to 24 hours to receive the email saying that you can download it.
Click Shop OpenStreetMap data now. You are now at the Trimble Data Marketplace.
Sign In (upper right corner). (Registration is free, and necessary.) Your name should now appear in the upper right hand corner, and below it, over the northwest portion of the map, you should see a summary of your current order, with headings for Region, Layers, Datum-Projection and File Format.
Initially the map shows the entire world (in a ghastly pseudo-mercator projection, but that does go hand-in-hand with OSM). Region is “Entire map,” Layers are “26/26” (26 of 26 available layers), Datum-Projection is “Lat/Long-WGS84 (Native),” File Format is “ESRI Shape,” Estimated size is “1.36 TB” and cost is “Free.
Don’t be fooled by the 1.36 TB. That’s the size of the entire world database for OSM. Trimble will restrict you to a 5 GB limit in a single order.
On the map you can now pan and zoom to the region you’re interested in. Let’s say you go to the area around Budapest, Hungary.
Notice that even though you are zoomed in, the Region is still “Entire Map” (meaning the whole world) and Estimated Size is still 1.36 TB. To narrow down your data area you have to either upload a KML file with a polygon of your area of interest, or draw a polygon on the screen.
To draw a polygon, click the pencil icon next to Region, click Draw, and begin placing points around your polygon. Clicking on the initial point closes the polygon. Now Region and Estimated Size change to something more reasonable.
You may not want all 26 layers, and you can click the pencil icon next to Layers to select the subset you want. The 26 layers are explained in depth at the OSM wiki page on Map Features, but for simplicity I’ll just list the 26 categories here, with links to the OSM wiki page:
Note that these categories are often subdivided into separate shapefiles for point, line and polygon features.
In this case let’s say I want only Waterway, Highway, Railway. I deselect all, select these three, and click the “X” to close the Layers list. Layers has now changed to “3/26” and Estimated Size has now dropped to 381 MB.
Then click Order. Accept the Content License and click Add To Cart. Click View Cart.
Everything should still read “Free,” so click Checkout. Go through the [annoying] Address Information and Select Payment steps (“No Payment Required”) and then finally Place Order.
The next step is that you receive a series of emails acknowledging your order, telling you that your order is being prepared, and finally that your order is ready for download. It can be almost immediate for small orders, or up to 24 hours for large ones.
When you do download your data it will be in a ZIP file called “weogeo_<order number>.zip.”