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:
+proj=lcc +lat_1=36 +lat_2=52.8 +lat_0=40 +lon_0=31.5 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs
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.
There is no new raster image with this linear transformation, just a small text file with the same name as the map and a .WLD extension.
To clip off the collar of the georeferenced map, in case I want to mosaic it with adjacent maps, I can create a polygon in a temporary layer that covers the part of the map I want to keep…
and then do Raster>Extraction>Clip raster by mask layer. I like to check Create an output alpha band and Keep resolution of output raster.
And the result is a clipped raster with transparency in band 4.