2

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?

  • 4
    The degrees are not WGS84, but Pulkovo 1942. – AndreJ Dec 4 '18 at 7:18
  • @AndreJ I mentioned this in the original question, and that I have been using it as the input CRS. – djpeanut Dec 4 '18 at 7:34
  • 1
    The map is projected so you should be georeferencing to the grid not the graticule. Have a close look with a ruler and you'll see the graticule bends slightly but the grid is straight. The Soviet era maps aren't top secret, far from it, they are a bit dated now but if that's all you have covering an area they are fairly good (better than nothing). – Michael Stimson Dec 4 '18 at 7:38
  • 1
    @djpeanut this "input CRS" has no impact on the georeferencing. Your target CRS should be EPSG:4284 instead of 4326, because you are using russian latlon degrees, not WGS84 ones. Have a look at epsg.io/4284 to see what transformations are available. – AndreJ Dec 4 '18 at 12:48
  • 1
    @djpeanut please self-answer your own question in the answer section below, for the benefit of future readers. For the JOSM problem, you might open a new topic. – AndreJ Dec 4 '18 at 17:08
2

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 Sheets of different scales

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. Boundaries of 1:25 000 trapeziums

  • Thank you very much for this. You might have just sped up my workflow by several orders of magnitude :) I'll get back to you if any questions... – djpeanut Dec 5 '18 at 4:52
  • You are welcome. Ready to answer any question. – Vadym Dec 5 '18 at 5:22
  • select * from __bvv_sk42_bound_100 ('M-38-144') - why does the algorithm not work? You tested it? – Cyril Dec 6 '18 at 14:34
  • On a sheet frame of scale 1: 100 000 only 6 points (4 corners and 2 so-called deflection points)... – Cyril Dec 6 '18 at 16:24
  • Did you add all the three functions in your PostgreSQL database? – Vadym Dec 6 '18 at 17:02
1

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!

  • That means the coordinates at the four corner of the paper map were in lats/longs (???). Confusion arose (for me) because EPSG 28408 is in meter. It may also mean the seven parameters is not the cause of the problem. – Ralph Tee Dec 5 '18 at 5:35
  • @RalphTee Yes, the printed co-ordinates are in dd mm ss format – djpeanut Dec 5 '18 at 5:46
  • 1
    @RalphTee The map bounds are in degrees, but the GK meter grid and coordinates are also printed along the collar. You can reference to both systems. While P42 GK zones always have the 7-parms datum shift, selecting EPSG:4326 for degrees has not, while EPSG:4282 has. – AndreJ Dec 5 '18 at 7:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.