12

Background This is my second question related to georeferencing naked raster maps in order to re-visualize them on different coordinate systems and in conjunction with other data layers. The previous question is at Convert an arbitrary meta-data-free map image into QGIS project

Problem My goal is to georeference this map:

Eurasian steppe, Encyc. Brit.?

This does not appear to be Plate-Carrée. So in QGIS, I created several reasonable control points, which for completeness I have attached at the bottom [ref:1]. I provide QGIS Georeferencer the same target SRS as my project file, EPSG:4326. I get exceptionally poor results with Helmert and the polynomial transforms but get a reasonable image with thin plate spline (which makes the resulting geoestimate go through my control points). However, even this result is poor, e.g., at higher latitudes (see the Russian coast north of Japan). This is a screenshot of my QGIS screen using a Natural Earth background.

QGIS georeference result, thin plate spline

Alternative path I tried a similar exercise with the much easier-to-use tool at MapWarper: see the result and control points at http://mapwarper.net/maps/758#Preview_Map_tab where I get poorer results (probably due to the fact that I added fewer control points).

Questions in nutshell

  1. Are there are any tricks I'm missing to getting a good georeference?
  2. Is this projection instantly recognizable?
  3. At Unknown Coordinate System on old drawing, gdaltransform is suggested to transform several coordinate points into a several target SRS, with the goal of actually uncovering the projection parameters used to generate the original map. I tried something like this: after saving my QGIS list of points, I did some string processing to get a list of space-separated long/lats via cat eurasian-steppe-gcp.points | tail -n+2 | cut -d, -f1-2 | sed 's/,/ /'> tmp.txt and inputting the resulting file into gdaltransform: gdaltransform -s_srs EPSG:3785 -t_srs EPSG:4326 < tmp.txt and switching the s_srs and t_srs flags (the project uses EPSG:4326). I know I'm shooting in the dark, hoping to get lucky, so I wasn't surprised when I couldn't make sense of the outputs.

Can someone expand on how I would use this method to find the best estimate of the source map's projection and projection parameters?

My thinking behind this is that rather than messing with placing myriad control points for a good georeference, might it be easier to get a near-perfect georeference with fewer control points, just looping through all the common coordinate systems?

Does it involve cross-validation of each point against all the others, for each CRS under test?

I'd like to get an understanding of either this algorithm or of georeferencing so I can automate the process---I run into this issue all the time, and until content creators stop treating their maps as one-off creations never to be integrated with other content, I don't expect to stop.

References

[ref:1] QGIS GCP file:

mapX,mapY,pixelX,pixelY,enable
142.632649100000009,54.453595900000003,505.941176470588232,-95.220588235293974,1
154.934252200000003,59.559921699999997,536.411764705882206,-52.779411764705742,1
80.080158100000006,9.657192300000000,291.558823529411711,-322.661764705882206,1
10.448442600000000,57.819128900000003,21.676470588235190,-103.926470588235134,1
34.007173000000002,27.761438299999998,101.117647058823422,-244.852941176470466,1
50.950890399999999,11.862196600000001,171.852941176470495,-313.955882352941046,1
29.713217199999999,60.024133200000001,90.779411764705799,-92.499999999999829,1
60.000000000000000,0.000000000000000,208.308823529411683,-362.382352941176350,1
69.867506500000005,66.639146199999999,224.088235294117567,-33.191176470588061,1
27.276107100000001,71.049154799999997,89.147058823529306,-21.764705882352814,1
140.000000000000000,0.000000000000000,536.955882352941217,-362.926470588235190,1
20.000000000000000,0.000000000000000,43.441176470588132,-362.926470588235190,1
20.196882700000000,31.243024100000000,47.249999999999901,-231.794117647058698,1
9.171861099999999,42.848309999999998,8.073529411764603,-175.205882352941046,1
131.955786100000012,43.196468600000003,481.999999999999943,-162.691176470588090,1
73.813303700000006,45.169367200000003,256.735294117646959,-161.602941176470438,1
50.602731800000001,44.589102900000000,168.044117647058727,-167.588235294117510,1
121.394975900000006,18.941421099999999,455.882352941176407,-284.029411764705742,1
103.987047000000004,1.417439300000000,389.499999999999943,-357.485294117646959,1
109.325478599999997,55.962283100000001,380.249999999999943,-98.485294117646902,1
31.454010100000001,46.562001500000001,95.132352941176379,-158.882352941176322,1
43.639560299999999,68.844150499999998,137.573529411764611,-40.264705882352814,1

Non-update

Analysis of van der Grinten I wrote a Python tool to fit GCPs to any projection that Proj4 supports (via Pyproj) and applied it to the couple of projections suggested in the answers. The source code as well as updated GCPs are available at https://github.com/fasiha/steppe-map

The van der Grinten has only 1 parameter to tune, and here's the resulting image (using the latest image from Britannica, which is a high-res and updated map (though it still lacks projection data)).

Van der Grinten fit

Van der Grinten has a relative error of 0.035 between the GCPs and best-fit points, which is the worst of the bunch I tried, and the coastline overlay bears that out qualitatively.

(It may help if you open this image in its own tab, it's quite high-res. You'll also see green arrows indicating the georeferenced points (they should match significant landmarks on the image) as well as red arrows indicating where those points are fitted to (they should match the same landmarks on the coastline overlay)---the deviation between the two can help the eye see the differences between the image and the fit.)

Analysis of Albers equal-area Trying the same thing with the Albers equal-area projection (which I think is the same as "Albers conformal Conic"). This fit, involving a 4-dimensional parameter fit, is better, with a relative error of 0.025, but looks pretty poor nonetheless.

Albers equal-area fit

Analysis of Robinson and Eckert V projections I fit a number of pseudocylindrical projections supported by Pyproj (all that I could find that had one free parameter) and found that the Robinson and the Eckert V projections did the "best" in terms of relative error between the GCPs and the fitted points, both with relative errors of 0.015.

Here's the Robinson:

Robinson fit

And here's the Eckert V.

Eckert V fit

Note the deviations of the fitted coastline from the image's coastline. I think with this I can conclude that the map is none of the above.

Winkel tripel: possible solution?

After sequentially trying every projection in this Proj manual from 1990 (updated 2003) ftp://ftp.remotesensing.org/proj/OF90-284.pdf I finally came to the Winkel tripel projection. This produces the lowest quantitative errors (0.011) and the coastline is uniformly quite good (or equivalently, uniformly slightly bad). I read that this is the projection of the National Geographic Society, which means it's famous, and this adds weight to the candidacy of this projection for Britannica's map. The fitted SRS: +units=m +lon_0=47.0257707403 +proj=wintri.

Winkel tripel fit

I will try to tweak my GCPs to try and drive the error down lower.

0

2 Answers 2

3

From a first analysis, we can see that only equator and the central meridian (passing approximatively through Moscow, about 37.36 E of longitude) are straight lines, while all the other parallels and meridians are curve lines. Among the global projections, these features suggest to try the Van der Grinten one (as bhell did it in a comment of the question), i.e. the standard world map projection of the National Geographic Society.

So I've set the project SR in QGIS as custom SRS World_Van_der_Grinten_I (Moscow):

+proj=vandg +lon_0=37.36 +x_0=0 +y_0=0 +R_A +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

derived directly from EPSG:54029 (World_Van_der_Grinten_I), changing only the longitude of the central meridian. Then I've tried to georeference the image using the projective transformation with the following GCPs chosen on an imaginary graticule passing through the available tick marks:

mapX,mapY,pixelX,pixelY,enable
-1930600.82502013398334384,-91.97837987072125543,43.52016129032256231,-362.54032258064518146,1
11412132.24015056900680065,-91.97837987233651802,537.54173867127497033,-362.73421498975926625,1
11804723.15208163484930992,8693410.00194966979324818,551.40707210100163138,-43.85864138853239069,1
-1651013.74878141912631691,7637895.72632359713315964,54.16126727366525984,-80.84230997717837397,1

This is the result in the custom World_Van_der_Grinten_I (Moscow) with GCPs (world country admin boundary in red):

enter image description here

This georeference could be further improved choosing better accurate GCPs (the great difficult is the absence of the tick marks or a graticule), however it seems that map contains some artefacts especially in Central and Northern Europe. For instance, consider Germany which seems collapsed in reference to Italy or the relative position of the Kolguyev Island (that fits quite well) and the Russian coast (that doesn't fit at all).

The same result reprojected in WGS84 (EPSG:4326):

enter image description here

0
1

there is a strong distortion on the east ("going up"): my guess is a Albers conformal Conic. Then the vertical meridian seems to be at 40° (-> your central meridian) and the horizontal parallel probably around 40°N. You should then measure your XY coordinate from the axis at the intersection of meridian 40 and parallel 40, then have a try...

Note : The parallels are not parallel, so you can exclude Cylindrical (Mercator) and pseudo-cylindrical projections (sinusoidal, Eckert)

0

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