1

I'm trying to transform this shaded map of Venus, a 2000x1000 TIFF in what I believe is equirectangular projection, without GIS metadata, to Mercator.

original image

I first tried

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395 \
  venus_shaded.tif venus_shaded-mercator.tif

but this fails with

ERROR 1: The transformation is already "north up" or a transformation
between pixel/line and georeferenced coordinates cannot be computed for
venus_shaded_cyl.tif. There is no affine transformation and no GCPs.
Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this
check.

Simply adding SRC_METHOD=NO_GEOTRANSFORM is clearly not the answer, as it produces garbage:

failed warp

I then tried using gdal_translate to add GCPs:

gdal_translate \
  -gcp 0 0 0 90 \
  -gcp 0 333 0 30 \
  -gcp 0 666 0 -30 \
  -gcp 0 999 0 -90 \
  -gcp 666 0 120 90 \
  -gcp 666 333 120 30 \
  -gcp 666 666 120 -30 \
  -gcp 666 999 120 -90 \
  -gcp 1333 0 -120 90 \
  -gcp 1333 333 -120 30 \
  -gcp 1333 666 -120 -30 \
  -gcp 1333 999 -120 -90 \
  venus_shaded.tif venus_shaded-gcp.tif 

(and then running the same gdalwarp command on the new file), but the results make it clear I don't know what I'm doing.

enter image description here

I don't know whether the problem is that I'm misunderstanding the format of the -gdp option to gdal_translate, misunderstanding how to apply it to pixel coordinates, misunderstanding the whole concept of GCPs, or all of the above. It's also possible EPSG:4326 is the wrong source coordinate system.

Regardless, I feel like this should be a pretty straighforward operation, and there's probably a better way to do it. What is that way?

2
  • what does gdal_translate give for venus_shaded.tif now? You probably don't need s_srs
    – mdsumner
    Commented Feb 9, 2020 at 22:59
  • Ah no, you'll need to set -ullr argument with gdal_translate to give the extents in the source - you won't need GCPs if it's a simple georef - of course the earth's datum won't be suitable, but you'll get close enough for casual use
    – mdsumner
    Commented Feb 9, 2020 at 23:02

1 Answer 1

1

This case of referencig is linear, because you are referencing a rectangular grid (row, pixel) to a rectangular grid (latitude, longitude).

The easiest way to do it, is with a World File.

Create a new text document and paste the following lines in it:

0.18
0.0
0.0
-0.18
-179.91
89.91

Save it as venus_shaded.wld, in the same folder as venus_shaded.tif.

0.18 is the pixel cell width and -0.18 is the pixel cell height.
(89.91, -179.91) are the referenced coordinates of the first row first pixel cell center.

Let's check the referencing:

C:\GA\GIS\Pruebas\test>gdalinfo venus_shaded.tif
Driver: GTiff/GeoTIFF
Files: venus_shaded.tif
       venus_shaded.wld
Size is 2000, 1000
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.180000000000000,-0.180000000000000)
Metadata:
  TIFFTAG_DOCUMENTNAME=test.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000)
Lower Left  (-180.0000000, -90.0000000)
Upper Right ( 180.0000000,  90.0000000)
Lower Right ( 180.0000000, -90.0000000)
Center      (   0.0000000,   0.0000000)
Band 1 Block=2000x1 Type=Byte, ColorInterp=Red
Band 2 Block=2000x1 Type=Byte, ColorInterp=Green
Band 3 Block=2000x1 Type=Byte, ColorInterp=Blue

Seems good to me.

Leaving aside that this map is not of the Earth, just note that Mercator projection tends to infinite at the poles.

GDAL knows it, and lets us warp the whole spheroid to Mercator, clipping the poles:

C:\GA\GIS\Pruebas\test>gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395 venus_shaded.tif venus_mercator.tif
Creating output file that is 1738P x 1407L.
Processing venus_shaded.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

C:\GA\GIS\Pruebas\test>gdalinfo venus_mercator.tif
Driver: GTiff/GeoTIFF
Files: venus_mercator.tif
Size is 1738, 1407
Coordinate System is:
PROJCRS["WGS 84 / World Mercator",
    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["World Mercator",
        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",0,
            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["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - between 80┬░S and 84┬░N"],
        BBOX[-80,-180,84,180]],
    ID["EPSG",3395]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,16165130.295636426657438)
Pixel Size = (23050.260090030031279,-23063.530245522277255)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DOCUMENTNAME=test.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-20037508.343,16165130.296) (180d 0' 0.00"W, 80d59'27.57"N)
Lower Left  (-20037508.343,-16285256.760) (180d 0' 0.00"W, 81d 9'30.36"S)
Upper Right (20023843.694,16165130.296) (179d52'38.09"E, 80d59'27.57"N)
Lower Right (20023843.694,-16285256.760) (179d52'38.09"E, 81d 9'30.36"S)
Center      (   -6832.325,  -60063.232) (  0d 3'40.95"W,  0d32'35.47"S)
Band 1 Block=1738x1 Type=Byte, ColorInterp=Red
Band 2 Block=1738x1 Type=Byte, ColorInterp=Green
Band 3 Block=1738x1 Type=Byte, ColorInterp=Blue

Note this:

Corner Coordinates:
Upper Left  ... (180d 0' 0.00"W, 80d59'27.57"N)
Lower Right ... (179d52'38.09"E, 81d 9'30.36"S)

You can open a question about that seemingly arbitrary clipping.

At this moment, the referencing issue was solved.

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.