2

I have a map file in a tiff format that is not geo-referenced, but I know the following projection information: +proj=lcc +lon_0=10.856 +lat_0=51.322 +lat_1=42 +lat_2=57 +datum=WGS84. I would like to transform it into EPSG:3857.

In my first attempt I tried to fix four arbitrary points and perform warping

gdal_translate -a_srs "+proj=lcc +lon_0=10.856 +lat_0=51.322 +lat_1=42 +lat_2=57 +datum=WGS84" -of GTiff -gcp 19725 11865 1 56 -gcp 103755 12990 23 56 -gcp 112755 80400 23 46 -gcp 8925 78990 1 46 original.tif translated.tif

But the coordinates of the corners of the warped image where not as expected and the map was distorted. 

Corner Coordinates:
Upper Left  ( 1208480.629, 6678543.890) ( 10d51'21.48"E, 51d19'21.08"N)
Lower Left  ( 1208480.629, 6678541.123) ( 10d51'21.48"E, 51d19'21.03"N)
Upper Right ( 1208486.180, 6678543.890) ( 10d51'21.66"E, 51d19'21.08"N)
Lower Right ( 1208486.180, 6678541.123) ( 10d51'21.66"E, 51d19'21.03"N)
Center      ( 1208483.404, 6678542.507) ( 10d51'21.57"E, 51d19'21.06"N)

After reading many posts, here and elsewhere, I came to this:

  1. Converting the UL and LR corners to x and y
cs2cs +proj=latlong +lon_0=10.856 +lat_0=51.322 +lat_1=42 +lat_2=57 +datum=WGS84 +to +init=epsg:3857

-3.63 57.28 >   -404089.75 7817564.89
23.76 44.39 >   2644951.10 5525995.28
  1. Geo-referencing the tiff
gdal_translate -a_srs "+proj=lcc +lon_0=10.856 +lat_0=51.322 +lat_1=42 +lat_2=57 +datum=WGS84" -a_ullr -404089.75 7817564.89 2644951.10 5525995.28 original.tif translated.tif

The result of gdalinfo translated.tif is as follows:

Size is 14735, 11333
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",42],
    PARAMETER["standard_parallel_2",57],
    PARAMETER["latitude_of_origin",51.322],
    PARAMETER["central_meridian",10.856],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-404089.750000000000000,7817564.889999999664724)
Pixel Size = (206.925066168985410,-202.203265684284787)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=144
  TIFFTAG_YRESOLUTION=144
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -404089.750, 7817564.890) (146d18'53.31"E, 73d27'53.56"N)
Lower Left  ( -404089.750, 5525995.280) (158d45'33.92"W, 88d 1'24.64"N)
Upper Right ( 2644951.100, 7817564.890) (172d26'19.88"W, 64d26'59.68"N)
Lower Right ( 2644951.100, 5525995.280) (138d13'57.20"E, 73d22'12.76"N)
Center      ( 1120430.675, 6671780.085) (161d52'17.63"W, 79d37'35.75"N)

However, the coordinates of the corners above are completely off. When I attempted to warp it, the image could become so huge that I had to stop the process after a couple of minutes. 

gdalwarp -overwrite -t_srs EPSG:3857 -r near -co COMPRESS=LZW translated.tif warped.tif

The result of gdalinfo warped.tif is as follows:

Corner Coordinates:
Upper Left  (-20037508.339,29704755.766) (180d 0' 0.00"W, 88d54'44.28"N)
Lower Left  (-20037508.339, 9464937.842) (180d 0' 0.00"W, 64d26'59.58"N)
Upper Right (20037503.600,29704755.766) (179d59'59.85"E, 88d54'44.28"N)
Lower Right (20037503.600, 9464937.842) (179d59'59.85"E, 64d26'59.58"N)
Center      (      -2.370,19584846.804) (  0d 0' 0.08"W, 84d41'15.52"N)

What did I miss?

|improve this question
  • Try projecting the lat-lon points to the Lambert system, then use those for the input gcp. I'm also not sure that gdal_translate is the right function. Maybe gdalwarp which allows you to specify both the input and output coordinate systems? – mkennedy Jun 8 '18 at 20:04
  • I have tried your first suggestion with the same results. As for the second, you cannot set gcps with gdalwarp and I use gdal_translate to try geo-referencing the image. In any case gdalwarp alone doesn't solve it. – yaronh hnoray Jun 11 '18 at 13:09
1

The first thing to do is to re-read the documentation of gdal_translate https://www.gdal.org/gdal_translate.html and then save your ground control points correctly

-gcp pixel line easting northing

So start counting from the top-left pixel and feed first the location of the pixel (horizontal - vertical) and then the coordinates of that pixel. You did

-gcp 19725 11865 1 56

and that means "the pixel at location 19725 to the right and 11865 downwards from the top-left corner has coordinates 1 E, 56 N.

|improve this answer

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.