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
Improve this question
5
  • 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, 2014 at 21:00
  • Why don't you give your update as an answer?
    – AndreJ
    Dec 31, 2014 at 17:04
  • Would it help? mkennedy answered it, not me.
    – scruss
    Jan 1, 2015 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, 2015 at 5:03
  • Although I've now added an answer, I think the issue you raise above is more of a SE design problem than something users have to work around.
    – scruss
    Jan 3, 2015 at 2:23

1 Answer 1

Reset to default
3

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
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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