Setting proj parameters of old map

Question

I would like to georeference these fantastic old public domains maps https://legacy.lib.utexas.edu/maps/ams/italy_50k/

I started from this one, and it has these projection parameters.

Then I have tried to georeference it, using the commands below and this proj parameters +proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +ellps=bessel +units=m +no_defs +k_0=0.998992911

# delete existent images
rm ./output_gcp.tif
rm ./output_warped.tif
# download the image
curl -L "http://legacy.lib.utexas.edu/maps/ams/italy_50k/txu-pclmaps-oclc-6540719-milano-45-iii.jpg" >./output.jpg
# add control points
gdal_translate -of GTiff -gcp 532.292 2433.53 405000 551000 -gcp 2618.57 2316.82 424000 551000 -gcp 2388.9 467.338 423000 567000 -gcp 537.666 569.428 407000 567000 -gcp 933.404 1364.4 410000 560000 -gcp 2090.5 1299.92 420000 560000 -gcp 1462.72 517.946 415000 567000 -gcp 1576.18 2376.68 415000 551000 ./output.jpg ./output_gcp.tif
# apply the proj parameters
gdalwarp -r near -order 1 -co COMPRESS=PACKBITS  -t_srs "+proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +ellps=bessel +units=m +no_defs +k_0=0.998992911" ./output_gcp.tif ./output_warped.tif

But I have a not so good result, a displacement of about 150 meters.

I think that the problem here is the definition of the standard parallel/s.

With this kind of projection (as I have read here), "two parallels and a scale of 1 are equivalent to one parallel somewhere between them and a scale less than 1. Since there is a scale, which is slightly less than 1, I think that there is one parallel"; it should be the Origin I can read in the projection info (45°54'N), and I have used it.

But it seems not work. Probably I must set two standard parallels, but which should i use?

EDIT

I have tried to use geographical coordinates of the 4 corners (using EPSG:1660):

# download the image
curl -L "http://legacy.lib.utexas.edu/maps/ams/italy_50k/txu-pclmaps-oclc-6540719-milano-45-iii.jpg" >./output.jpg
# add control points
gdal_translate -of GTiff -gcp 418.832 2486.39 -3.5 45.3333 -gcp 2689.8 2501.38 -3.25 45.3333 -gcp 2686.8 349.879 -3.25 45.5 -gcp 425.704 337.633 -3.5 45.5 ./output.jpg ./output_mm.tif
# apply the proj parameters
gdalwarp -r near -order 1 -co COMPRESS=PACKBITS  -t_srs "+proj=longlat +ellps=intl +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +pm=rome +no_defs" ./output_mm.tif ./output_mm_warped.tif

It's a better result, but I would like to find a right way using the Lambert Projction.

Moreover it's wrong assumption because "The Monte Mario datum based on the Hayford (International 1909) ellipsoid was adopted in 1940, but this map explicitly says that it is projected from the Bessel ellispsoid" (see this reply).

Answer 1

The main problem is in the image warping parameters, they are not giving you the right gridded image.

About the projection itself, you are using the Lambert Conical projection in the right way. You have just not information about the datum transformation.

A datum transformation is going to be necessary if we want to compare this map with one based on the WGS84 datum.

About warping the image from the geographic coordinates of the corners with a linear transformation seems to me a wrong process, beyond that the result looks better than the previous one.

---

I have loaded the original .jpg file and digitized your Ground Control Points:

---

They seems good. But when you warp the image, you are forcing computation to first degree polynomials (-order 1), which is something like performing a linear transformation in the image (the overabundance of points is used only in the compensation). This is the result with a rectangular grid superimposed:

---

Since eight GCPs are enough to a second order polynomial, this is the result with -order 2 parameter (or without it, because it will be the default for 8 GCPs):

---

Seems better, but the scanner has its deformations. This is the result forcing the Thin Plate Spline transformation (-tps):

---

That's the best that we can do with those GCPs. Let's see how it is georeferenced:

---

The Monte Mario datum based on the Hayford (International 1909) ellipsoid was adopted in 1940, but this map explicitly says that it is projected from the Bessel ellispsoid.

From Geodetic datums of the Italian cadastral systems:

In the first decades of the twentieth century, the Italian Institute of Military Geography (Istituto Geografico Militare; I.G.M) developped four geodetic network, all on the Bessel 1841 ellipsoid; the Genova 1902, the old Monte Mario, the Castanea delle Furie 1910 and the Guardia Vecchia datums, for northern, central and southern Italy and Sardinia, respectively ( Mori, 1922 ).

In the section 9.1.3.4 of Practical Geodesy Using Computers, from Marteen Hooijberg there are the Helmert parameters derived from WGS84 to Genoa 02:

(A conversion from radians to seconds and the scale factor to p.p.m. is needed.)

The authors of the previous paper derived different parameters, to WGS84:

dX= +656.5 m; dY= +138.2 m; dZ= +506.5 m; rotX= –5.187 arc sec; rotY= +2.540 arc sec; rotZ= –5.256 arc sec; scale factor= –12.61 ppm (according to ‘coordinate frame rotation’; average horizontal error, except Sicilia and Southern Italy: 2.5 m; maximum horizontal error: 5.6 m.)

(We need to invert the rotation parameters because the convention for the -towgs84 parameter is "position vector" instead of "coordinate frame".)

gdalwarp -r near -tps -co COMPRESS=PACKBITS -t_srs "+proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +k_0=0.998992911 +ellps=bessel +towgs84=656.5,138.2,506.5,5.187,-2.540,5.256,-12.61 +units=m +no_defs" ./output_gcp.tif ./output_warped.tif

You will check the accuracy for the whole image, it is really very good.

Answer 2

Looking at this page which proposes this proj string for this Italy South grid:

+proj=lcc +lat_0=39.5 +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 
+y_0=600000 +ellps=bessel +units=m +no_defs

Seems to indicate you are on the right track, that may be as accurate as the projection gets?