I'm helping a friend convert his late father's war diary coordinates to WGS84. I have a proj.4 definition that appears to work in QGIS, but if I transform a scanned map with gdalwarp, it's not even close. What am I doing wrong with my CRS definition, please?
Coordinates were given in the Modified British System, which is a variant of the OSGB grid, but extended over several zones in Europe. Projection parameters (given to me by the linked site's author) are as follows:
Projection - Lambert Conical Orthomorphic
Ellipsoid: Bessel 1841
False Easting : 700000
False Northing : 600000
Central Meridian : 14.0°
Central Parallel : 39.5°
Scale Factor : 0.99906
Everything except the scale factor can be confirmed from 100k Index to WWII topo maps of Italy at McMaster University. From this, I developed a proj.4 definition:
+proj=lcc +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel +units=m +no_defs
After converting alphanumeric coordinates to grid coordinates and using the above definition, I get expected values for some sample locations:
- rN175845 → 717500, 784500 → 41.1629 °N, 14.2086 °E
- rN138862 → 713800, 786200 → 41.1783 °N, 14.1646 °E
When I try to use gdalwarp to translate a georeferenced map sheet to WGS84, however, I end up with latitudes around 3-4 °N. The command line I used is:
gdalwarp -s_srs "+proj=lcc +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000 +y_0=600000 +ellps=bessel +units=m +no_defs" -t_srs EPSG:4326 -r cubic -of GTiff in_modified.tif out-4326.tif
gdalinfo reports the following for the input file:
Size is 6641, 5162
Coordinate System is `'
Origin = (694095.553313978714868,954540.625156613648869)
Pixel Size = (8.495191956798934,-8.495191956798934)
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 694095.553, 954540.625)
Lower Left ( 694095.553, 910688.444)
Upper Right ( 750512.123, 954540.625)
Lower Right ( 750512.123, 910688.444)
Center ( 722303.838, 932614.535)
Band 1 Block=6641x1 Type=Byte, ColorInterp=Red
Band 2 Block=6641x1 Type=Byte, ColorInterp=Green
Band 3 Block=6641x1 Type=Byte, ColorInterp=Blue
(the sample coordinates are on another map sheet, btw.)
Update: a working proj.4 definition, based on mkennedy's comment below, is:
+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
Answer
Expanding on mkennedy's comment, and translating generic map terms to proj.4 parameters:
Projection - Lambert Conical Orthomorphic → +proj=lcc
Ellipsoid: Bessel 1841 → +ellps=bessel
False Easting : 700000 → +x_0=700000
False Northing : 600000 → +y_0=600000
Central Meridian : 14.0° → +lon_0=14
Central Parallel : 39.5° → +lat_0=39.5 +lat_1=39.5
Scale Factor : 0.99906 → +k_0=0.99906
(other proj.4 terms) +units=m +no_defs
which results in a working definition of
+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
No comments:
Post a Comment