Producing offset after georeferencing Soviet topographical military maps from printed control point coordinates in QGIS

Question

I have 'acquired' a collection of scans of (previously top-secret) Soviet-era military topographical maps at a set of 1:25.000. I want to georeference a set of these scans.

I know the coordinate system and projection (Pulkovo 1942, Gauss-Kruger Zone 8, EPSG:28408) and am selecting this when pulling the scans into QGIS's Georeferencer Plugin.

To do the georeferencing, I'm using the coordinates printed at each corner of the map, which has been standardised at 7.5' divisions of latitude and 5' divisions of longitude.

I'm using WGS84 as the output coordinate system. However, when I then load the resulting GeoTIFF, it produces a very small latitude offset (approx 2") when compared to all other data I can find, including my own GPS traces and GCPs.

Within the area I'm interested in, the offset seems to be the same for every map, but I'm wondering if the issue is with my georeferencing process or if the maps themselves were all produced with this error?

Answer 1

As a matter of fact I georeferenced a lot of such topomaps and completely agree with the comments. In order to speed up the process of georeferencing you should develop a tool in Qgis or in another GIS that generates border of the scanned map by its nomenclature number in a layer created in the mentioned projection Pulkovo 1942, Gauss-Kruger Zone 8, EPSG:28408. If you would like to have a program code for generating coordinates of frames by nomenclature number I can provide you with it in Python for example. In fact, there are three functions have been implemented in plpgsql. I tested them using PostgreSQL 9.6 and Postgis 2.5 and invoked them in Qgis 3.2 Examples of usage and the code.

For making boundary of 1:100 000 sheet select * from __bvv_sk42_bound_100('M-38-134')

CREATE OR REPLACE FUNCTION public.__bvv_sk42_bound_100(
    nom_number text)
    RETURNS TABLE(nom text, x1 double precision, y1 double precision, zone integer, epsg integer, trapezium geometry) 
    LANGUAGE 'plpgsql'

    COST 100
    VOLATILE 
    ROWS 1000
AS $BODY$

DECLARE
    /*Bounds of the sheet in coordinate system 1942*/
    Low_mil_lat double precision;
    Lef_mil_lon double precision;
    nl integer;
    rw integer;
    x1 double precision;
    y1 double precision;
    x2 double precision;
    y2 double precision;
    epsg integer;
    zone integer;
    col_number integer;
BEGIN
    col_number = cast(substring(nom_number from 3 for 2) as integer);
    Low_mil_lat = (ascii(upper(nom_number)) - 65) * 4.0;
    Lef_mil_lon = (col_number - 31.0)*6;

    nl = cast(substring(nom_number from 6 for 3) as integer);

    rw = (nl - 1) / 12;

    x1 =  Lef_mil_lon + (nl - rw * 12 - 1) * 0.5;
    y1 =  (Low_mil_lat) + 4 - (rw + 1) * (1.0/3);
    x2 = x1 + 0.5;
    y2 = y1 + 1/3.0;

    zone = col_number - 30;
    epsg = 28400 + zone;

    return query select nom_number, x1, y1, zone, epsg, ST_Transform(ST_MakeEnvelope(x1, y1, x2, y2, 4284), epsg);
END
$BODY$;

For making boundary of 1:50 000 sheet select * from __bvv_sk42_bound_50('M-38-134-A'). Only acceptable A, B, V, G

CREATE OR REPLACE FUNCTION public.__bvv_sk42_bound_50(
    nom_number text)
    RETURNS TABLE(nom text, longitude double precision, latitude double precision, zone integer, epsg integer, trapezium geometry) 
    LANGUAGE 'plpgsql'

    COST 100
    VOLATILE 
    ROWS 1000
AS $BODY$

DECLARE
    r record;
    lon double precision;
    lat double precision;
    d_lon double precision;
    d_lat double precision;
    letter text;
    col_number integer;
BEGIN
    /* 
    Returns bounding box of a sheet in scale 1:50000 as a
    quarter of a sheet in scale 1:100000
    */

    for r in select x1, y1 from __bvv_sk42_bound_100(nom_number)
    loop
        lon = r.x1;
        lat = r.y1;
    end loop;

    col_number = cast(substring(nom_number from 3 for 2) as integer);
    zone = col_number - 30;
    epsg = 28400 + zone;

    d_lat = 1.0/6;
    d_lon = 0.25;
    letter = substring(nom_number from 10 for 1);
    if letter = 'A'
    then
        lat = lat + d_lat;
    elsif letter = 'B'
    then
        lat = lat + d_lat;
        lon = lon + d_lon;
    elsif letter = 'G'
    then
        lon = lon + d_lon;
    end if;

    return query select nom_number, lon, lat, zone, epsg, ST_Transform(ST_MakeEnvelope(lon, lat, lon + d_lon, lat + d_lat, 4284), epsg);
END

$BODY$;

For making boundary of 1:25 000 sheet select * from __bvv_sk42_bound_25('M-38-134-A-a'). For the last letter only acceptable a, b, v, g.

CREATE OR REPLACE FUNCTION public.__bvv_sk42_bound_25(
    nom_number text)
    RETURNS TABLE(nom text, longi double precision, lati double precision, zone integer, epsg integer, trapezium geometry) 
    LANGUAGE 'plpgsql'

    COST 100
    VOLATILE 
    ROWS 1000
AS $BODY$

DECLARE
    r record;
    lon double precision;
    lat double precision;
    d_lon double precision;
    d_lat double precision;
    letter text;
    col_number integer;
BEGIN
    /* 
    Returns bounding box of a sheet in scale 1:25000 as a
    quarter of a sheet in scale 1:50000
    */

    for r in select longitude, latitude from __bvv_sk42_bound_50(nom_number)
    loop
        lon = r.longitude;
        lat = r.latitude;
    end loop;

    col_number = cast(substring(nom_number from 3 for 2) as integer);
    zone = col_number - 30;
    epsg = 28400 + zone;


    d_lat = 1.0/12;
    d_lon = 0.125;
    letter = substring(nom_number from 12 for 1);
    if letter = 'a'
    then
        lat = lat + d_lat;
    elsif letter = 'b'
    then
        lat = lat + d_lat;
        lon = lon + d_lon;
    elsif letter = 'g'
    then
        lon = lon + d_lon;
    end if;

    return query select nom_number, lon, lat, zone, epsg, ST_Transform(ST_MakeEnvelope(lon, lat, lon + d_lon, lat + d_lat, 4284), epsg);
END

$BODY$;

Below is some results of testing

I'd like to add some additional function that can help accelerate georeferencing process in coordinate system 1942.

CREATE OR REPLACE FUNCTION public.__bvv_generate_topo_25000(
    nom_num text)
    RETURNS void
    LANGUAGE 'plpgsql'

    COST 100
    VOLATILE 
AS $BODY$

DECLARE
    i integer;
    j integer;
    k integer;
    nom_100 text;
    nom_50 text;
    nom_25 text;
    trap geometry;
    r record;
    tn text;
    nn text;
    expr text;
    epsg integer;
    col_number integer;
    zone integer;
    ar_50 text[];
    ar_25 text[];

BEGIN
    ar_50 = Array['A', 'B', 'V', 'G'];
    ar_25 = Array['a', 'b', 'v', 'g'];

    col_number = cast(substring(nom_num from 3 for 2) as integer);
    SELECT to_regclass('public.topo_100000') into tn;
    if tn = 'topo_100000'
    then
        DROP TABLE public.topo_100000;
    end if;

    zone = col_number - 30;
    epsg = 28400 + zone;

    expr = format('CREATE TABLE public.topo_100000
    (
        nom_number text COLLATE pg_catalog."default" NOT NULL,
        trapezium geometry(polygon, %s),
        CONSTRAINT topo_100000_pkey PRIMARY KEY (nom_number)
    )', epsg);

    execute expr;

    for i in 1..144
    loop
        nn = lpad(cast(i as text), 3, '0');
        nom_100 := nom_num || '-' || nn;
        for j in 1..4
        loop
            nom_50 = nom_100 || '-' || ar_50[j];
            for k in 1..4
            loop
                nom_25 = nom_50 || '-' || ar_25[k];
                select trapezium from __bvv_sk42_bound_25(nom_25) into trap;
                insert into topo_100000 values(nom_25, trap);
            end loop;
        end loop;
    end loop;
END

$BODY$;

The function takes the name of 1:1 000 000 sheet. For example. In order to generate all the 1:25000 trapeziums that fall into the trapezium M-38 the function should be invoked like this.

select * from __bvv_generate_topo_25000('M-38')

After its performing the table topo_100000 2304 records is created.

Answer 2

Changing the target CRS in the Transformation Settings dialog of QGIS's Georeferencer plugin to EPSG:4284 (Pulkovo 1942) solved the problem and created a correctly referenced GeoTIFF as an output, which I could then load into Google Earth. The various Gauss-Kruger zones didn't seem to make any noticeable difference when reprojecting at such a large scale. Also, with these specific topo maps, after colour-correcting and cropping in Photoshop, I got the best results by using a Thin Plate Spline transformation with the cubic resampling method. Hope that helps!