4

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
  • 3
    Try setting lat_0 to 39.5 also (and maybe lat_2). The fact that you're getting very low latitude values makes me think that the latitude of origin is defaulting to zero. – mkennedy Dec 7 '14 at 21:00
  • Yes, that's it, thanks. lat_0 was all it needed! – scruss Dec 8 '14 at 5:22
  • Why don't you give your update as an answer? – AndreJ Dec 31 '14 at 17:04
  • Would it help? mkennedy answered it, not me. – scruss Jan 1 '15 at 2:46
  • Without an answer, the question will pop up regularily by the Community user, and we can not mark a new question as duplicate of this because it does not have an answer. – AndreJ Jan 2 '15 at 5:03
2

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

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.