1

I'm trying to use gdalwarp to project an FAA source chart in EPSG:9802 (Lambert Conic Conformal) to EPSG:3395 (World Mercator). The issue is that part of the data is east of longitude -180, and part is west. So when I try this:

/usr/bin/gdalwarp -t_srs EPSG:3395 -dstalpha -overwrite Lambert.tif Mercator.tif

I end up with a large, mostly empty image spanning 360 degrees of longitude, with the data small and distorted at the extreme left and right ends of the resulting TIFF.

I'm new to GDAL/projection but my research indicates that this might be because the central meridian of World Mercator is defined as 0 degrees.

I'm trying to find a way to use (for example) -180 degrees as the central meridian so that the resulting image in World Mercator is usable. If it makes it easier, I do not need the resulting World Mercator image to be georeferenced at all.

I've spent a couple of days trying various options to gdalwarp and gdal_translate but haven't made any progress. I've also spent a few hours trying to split the source TIFF into two (one west of -180 degrees and one east) and projecting them separately to World Mercator, so I could then manually combine the resulting images (which don't need to be georeferenced) with something like GIMP. Nothing has worked. Among various approaches I've tried --config CENTER_LONG 180, the -te option to gdalwarp, intermediate projection to Polar Stereographic, and others. Any ideas?

gdalinfo for Lambert.tif:


Driver: GTiff/GeoTIFF
Files: Lambert.tif
Size is 19587, 7468
Coordinate System is:
PROJCRS["AK_IFR",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257221999999,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["Lambert Conic Conformal (2SP)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",51,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-154,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",65,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",37,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
GeoTransform =
  -3327476.67088, 184.470564, 16.370996
  1268433.444758, 16.370185, -184.468772
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2022:12:07 10:59:53
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=GPL Ghostscript 8.71
  TIFFTAG_XRESOLUTION=400
  TIFFTAG_YRESOLUTION=400
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-3327476.671, 1268433.445) (152d33'50.09"E, 50d55'42.07"N)
Lower Left  (-3205218.073, -109179.345) (164d55'11.69"E, 41d26'35.27"N)
Upper Right (  285748.266, 1589076.258) (147d50'44.26"W, 65d27' 4.19"N)
Lower Right (  408006.864,  211463.469) (147d45'27.98"W, 52d47'50.30"N)
Center      (-1459734.903,  739948.457) (178d15' 4.16"W, 55d34'58.24"N)
Band 1 Block=19587x1 Type=Byte, ColorInterp=Red
Band 2 Block=19587x1 Type=Byte, ColorInterp=Green
Band 3 Block=19587x1 Type=Byte, ColorInterp=Blue
0

1 Answer 1

1

You could use a custom CRS which sets the central meridian to 180 instead of 0.

The proj4 definition of EPSG:3395 (WGS 84 / World Mercator) is:

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

If I change the central meridian to 180 and use that in the gdalwarp command instead of EPSG:3395:

gdalwarp -t_srs "+proj=merc +lon_0=180 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs" input.tif output.tif
gdalinfo output.tif

I get an output like so:

Driver: GTiff/GeoTIFF
Files: output.tif
       output.tif.aux.xml
Size is 20936, 14621
Coordinate System is:
PROJCRS["unknown",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Mercator (variant A)",
        METHOD["Mercator (variant A)",
            ID["EPSG",9804]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",180,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-3054171.234459184575826,9689316.348941301926970)
Pixel Size = (317.323442552925542,-317.323442552925542)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-3054171.234, 9689316.349) (152d33'50.09"E, 65d27' 4.19"N)
Lower Left  (-3054171.234, 5049730.295) (152d33'50.09"E, 41d26'37.55"N)
Upper Right ( 3589312.359, 9689316.349) (147d45'23.97"W, 65d27' 4.19"N)
Lower Right ( 3589312.359, 5049730.295) (147d45'23.97"W, 41d26'37.55"N)
Center      (  267570.562, 7369523.322) (177d35'46.94"W, 55d13'11.35"N)

enter image description here

Note: input.tif is just some dummy data with the same dimensions and CRS as Lambert.tif derived from the gdalinfo output included in the question.

1
  • Lifesaver! Thank you.
    – Roger
    Jan 11, 2023 at 2:57

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.