Bug in psconvert to GeoTiff?

Q&A
1

Hi,

I think I encountered a bug when converting to GeoTiff with psconvert. This is is with version 6.1.1 A minimum example demonstrating the issue is below.

First, this works and produces a tiff with the right coordinates

gmt grdimage -JM-81/46/15 -B0.1 -R-81.3/-80.9/46.6/46.8 -BneSW -P -Ctopo @earth_relief_01s > out.ps
gmt psconvert -W+g out.ps
gdalinfo out.tiff

Driver: GTiff/GeoTIFF
Files: out.tiff
Size is 2036, 1388
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[“Mercator_1SP”],
PARAMETER[“central_meridian”,-81],
PARAMETER[“scale_factor”,1],
PARAMETER[“false_easting”,0],
PARAMETER[“false_northing”,0],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,“9001”]]]
Origin = (-26421.455146726864768,4090887.571554347872734)
Pixel Size = (17.486072234763000,-17.483454968943999)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( -26421.455, 4090887.572) ( 81d14’14.45"W, 34d38’26.93"N)
Lower Left ( -26421.455, 4066620.536) ( 81d14’14.45"W, 34d27’37.61"N)
Upper Right ( 9180.188, 4090887.572) ( 80d55’ 3.12"W, 34d38’26.93"N)
Lower Right ( 9180.188, 4066620.536) ( 80d55’ 3.12"W, 34d27’37.61"N)
Center ( -8620.634, 4078754.054) ( 81d 4’38.79"W, 34d33’ 2.45"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue

Now omitting -P from grdimage

gmt grdimage -JM-81/46/15 -B0.1 -R-81.3/-80.9/46.6/46.8 -BneSW -Ctopo @earth_relief_01s > out.ps
gmt psconvert -W+g out.ps
gdalinfo out.tiff
Driver: GTiff/GeoTIFF
Files: out.tiff
Size is 2036, 1388
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[“Mercator_1SP”],
PARAMETER[“central_meridian”,-81],
PARAMETER[“scale_factor”,1],
PARAMETER[“false_easting”,0],
PARAMETER[“false_northing”,0],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,“9001”]]]
Origin = (-75102079391.892868041992188,75094956416.964126586914062)
Pixel Size = (17.486072234763000,-17.483454968943999)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-75102079391.893,75094956416.964) ( 94d27’28.33"W, 90d 0’ 0.00"N)
Lower Left (-75102079391.893,75094932149.929) ( 94d27’28.33"W, 90d 0’ 0.00"N)
Upper Right (-75102043790.250,75094956416.964) ( 94d 8’16.99"W, 90d 0’ 0.00"N)
Lower Right (-75102043790.250,75094932149.929) ( 94d 8’16.99"W, 90d 0’ 0.00"N)
Center (-75102061591.071,75094944283.446) ( 94d17’52.66"W, 90d 0’ 0.00"N)

So the coordinates of the corners are clearly wrong. The documentation for psconvert mentions that -W+g automatically sets -P but it does not explain that you will get wrong results if you omit it when creating the figure. So if this is not a bug, then I think it would be helpful to explain this in the documentation.

Regards

Max

2

Thanks, it’s a bug (this reminds me some old shit with the -P).
Even specifying -P directly doesn’t make a difference. But it works though when converting to EPS.

3

And that -P thing is a minor bug. The worst is that even with the -P the geographical coordinates are very wrong. It says lat is around 34 but -R says 46.6. This is so weird because the patagonia example in the psconvert man page works fine.

EDIT: The issue is related with a true scale != equator in the Mercator projection. See open issue #4193