I boiled the problem down to that python example:
import sys
try:
from osgeo import gdal
from osgeo import osr
except:
import gdal
import osr
if __name__ == '__main__':
version_num = int(gdal.VersionInfo('VERSION_NUM'))
print( "GDAL Version: " + str(version_num))
wkt = """
PROJCS["Popular Visualisation CRS / Mercator",
GEOGCS["Popular Visualisation CRS",
DATUM["Popular_Visualisation_Datum",
SPHEROID["Popular Visualisation Sphere",6378137,0,
AUTHORITY["EPSG","7059"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6055"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4055"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3785"]]"""
X = 2662981.191207001451403
Y = -1935374.062384711578488
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs, srsLatLong)
(int, lat, height) = ct.TransformPoint(X, Y)
# Report results
print('X: %s\t\tY: %s' % (X, Y))
print('Lat: %s\t\tLon: %s' % (gdal.DecToDMS(lat, 'Lat', 2), gdal.DecToDMS(int, 'Long', 2)))
my output is:
GDAL Version: 1900
X: 2662981.19121 Y: -1935374.06238
Lat: 17d 1' 1.01"S Lon: 23d55'19.08"E
However, when I run this from proj4, I get a different answer:
invproj +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
2662981.191207001451403 -1935374.062384711578488
23d55'19.081"E 17d7'29.797"S
The latitude is off by 6'28" or so... ?????? What could be going wrong? Many thanks, Pat.