Setting proj parameters of old map
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.
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?