..

How to calibrate maps

Ņemot vērā, ka šī tēma interesēs ļoti mazu skaitu manu lasītāju, rakstīšu rakstu angļu valodā. Ja tu esi viens no tiem, kuru tas tomēr interesē, uzdod jautājumus komentāros.


I love maps. I love looking at them, creating them, and doing all kinds of strange things to them. For a long time I had wondered, how it was possible to put old scanned maps online, on top of a modern map, to see how places look now compared to then. Maybe even identify lost places and roads, or finding objects that are now somewhere in the middle of woods. 

A few years ago, I found that the scanned map had to be georeferenced to GeoTIFF and then somehow cut into tiles, and hosted on the web using OpenLayers. I opened the website of Openlayers, and after a short while closed it. It was too much for me. 

Fast forward to a few months ago. I discovered Leaflet.js and subsequently also Mapbox.js, which is more or less the same thing, but prettier. 

Mapbox.js is a simple JavaScript file, which put into your page header can make a map. All you need is a few lines of HTML code. Like this:

<html>
<head>
<script src='https://api.tiles.mapbox.com/mapbox.js/v1.6.2/mapbox.js'></script>
<link href='https://api.tiles.mapbox.com/mapbox.js/v1.6.2/mapbox.css' rel='stylesheet' />
<style>
body { margin:; padding:; }
#map { position:absolute; top:; bottom:; width:100%; }
</style>
</head>
<body>
<div id='map'></div>
<script>
var map = L.mapbox.map('map', 'examples.map-9ijuk24y')
.setView([40, -74.50], 9);
</script>
</body>
</html>

And you have your map. To add markers, one more line. To add tracks, one more line. To load CSV markers from file — yes — just one more line. Mapbox.js has wonderful and simple documentation, with examples for any situation. 

So I created my first Leaflet/Mapbox webpage. I am not a programmer, I have no idea how to write JavaScript, I made it all with copy and paste, and a little help from some of my twitter followers. And so Dodies.lv was born. It is a not for profit website for hiking trails in Latvia. 

If you want to show a map, with anything on it, this is the way to go. 

So some time passed, and I wanted something new, something more. And then I remembered my old idea about hosting scanned historical maps. Hey, I rememered, Mapbox can do it! 

So all I needed was to scan, calibrate and tile. This will explain some of the technicalities of this process. It is not easy, but I will try to make it possible to repeat what I did. 

Let’s get down to business

Once you have obtained the scanned maps, you need to reference them. This means telling some program which point on the map is where in real life. In real life means WGS84 coordinate system, the standard that most mapping programs and GPS units use. There are a ton of coordinate systems, but that’s a horror story for another night. 

Georeferencing can be done in several ways, and most people recommend using QGIS. This is a very powerful tool (I have written about it before), but in my opinion, the referencing part is not that powerful or easy to use. There is another powerful tool, although not for free – OziExplorer. This tool is more meant for viewing maps, or driving around with them, but it also has a nice calibrator built in. Also, it’s map format is very popular in various applications, and the approach there is in my opinion easier to use. 

In OziExploer you clik “Choose and calibrate an image”, then find the edge of the map, where map creators have incidated the latitude and longitude of the map edge. Choose “point 1”, click on the corner of the map border, and enter the coordinates you can read from the map. Repeat this for all the four corners. It looks like this: 

You can also enable the cropping tool, which allows you to place markers where the map should be cut. This is useful for maps that have a white border around them, and I suggest you to use it. So for each map page you need to place 4 markers with coordinates, and 4 crop markers, where the map will be cut later. 

Of course this is not always as easy. Notice that in my image above, the coordinates are a little weird. The text next to them explains it, it says “East of Ferro“. This means that the map was made during a time when Greenwitch was not the Prime Meridian as it is today (the zero longitude).  And so for the example above, I need to substract the difference to know the longitude in todays format (39.30’00”-17.39’46”) and enter the result as the actual longitude. There might be a way to teach OziExplorer to understand I am using the Ferro longitude, but I haven’t been successful in doing that. 

So once you are done, save the calibration in the OziExplorer “.map” format, which fortunately is understood by many programs, such as GDAL. I can’t even begin to fathom the power of GDAL so I will not comment on what it is.  

Even after learning to do some things with it, GDAL to me still remains something akin to a drink of the gods. Or black magic.

GDAL itself can do image tiling, but I found a script that uses GDAL, but understands OZI explorer and can output tiles in Mapbox friendly folder structure (ZYX, as opposed to GDAL native TMS). 

Now to install GDAL, you have to many things. I suggest to follow guides made for your own OS. In short, for MacOS you need to install “Xcode5”, “HomeBrew” and then “brew install gdal”. You also need Python in your system, fortunately MacOS has it by default. 

So once you have calibrated map images with .map files alongside them, and you have a running gdal installation, you need the tiler_tools script. 

This is what you do next (assuming everything is running as it should, but I assure you, this will never be the case at first try. There are many dependencies and conflicts to solve to get GDAL running):

./gdal_tiler.py --cut --release *.map

The above command will read your calibration files, and then create a folder structure of tile images, for each zoom level. You will end up with thousands of small files, each in a separate folder for each map file called “somename.zyx”.  

If you have many such map files, you will need to combine the resulting folders into one, like this: 

./tiles_merge.py *.zyx karte

and if the resulting folder is huge, you can install pngquant, and compress the PNG files like this: 

find . -name '*.png' -print0 | xargs -0 -P4 pngquant --ext .png --force 256

Now all you need is to open up the resulting “karte” folder like in my example, and open the automagically created .html file to view your brand new creation, which should look similar to this:

Of course, the script included only OpenLayers and Google map examples, which is fine. But if you love Mapbox, like me, you might create a much smaller and more elegant HTML file with this content below the initial Mapbox stuff I was showing above: 

L.control.layers({
'Name of layer': L.tileLayer('/tiles1/z{z}/{y}/{x}.png').addTo(map),
'Name of other layer': L.tileLayer('/tiles2/z{z}/{y}/{x}.png')
}).addTo(map);

Isn’t it beautiful? 

As a bonus, I will show you another great tip 

Lets say you have obtained a list of points in CSV format. Just a table of points. You open it up, and – it turns out to be in some weird map projection, in my case, our local LKS92 projection. 

You could put these on a Mapbox map using a plugin, but why not convert them to WGS84 regular coordinates? 

If you have a csv file with fields “x,y” called “visivisi.csv”, make a new file called visivisi.vrt, and this inside: 

<OGRVRTDataSource> 
<OGRVRTLayer name="visivisi">
<SrcDataSource>visivisi.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs</LayerSRS>
<GeometryField encoding="PointFromColumns" x="x" y="y"/>
<Field name="latitude" src="x" />
<Field name="longitude" src="y" />
</OGRVRTLayer>
</OGRVRTDataSource>

To explain, this is a standard GDAL VRT file, it includes a line about the source CSV file name, and a line about the CSV file projection. In my case it is the LKS92 projection in “Proj4” format. Procections are easily found with your favorite search engine. 

Then, run this magic command:

ogr2ogr -overwrite -f CSV -lco GEOMETRY=AS_XY -t_srs EPSG:4326 test.csv visivisi.vrt

You will end up with a new file “test.csv” which will include additional columns, this time in normal latitude/longitude format (WGS84).