I found the problem....
The problem is in defining the geo_transform. I had the following:
x_rotation = 0
y_rotation = 0
cell_width_meters = 50
cell_height_meters = 50
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]
dataset.SetGeoTransform(geo_transform)
The Gdal documentation isn't real clear about what these values are. (See SetGeoTransform) Searching around the internets I derived that the passed values should be (in order):
- top_left_x
- cell_width_meters
- x_rotation
- top_left_y
- y_rotation
- cell_height_meters
Which seems right, BUT re-reviewing the GDAL API Tutorial I noticed that the last value, cell_height_meters was shown being given in a negative value. It seems that this was all that was needed to properly output the data in the expected orientation.
So now I've changed the geo_transform definition line to:
(Notice the added "-")
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, -cell_height_meters ]
It is not clear precisely what form the data are in, but ultimately any solution will have to use the equations for the Mercator projection of an ellipsoid. It sounds like it's possible to identify the latitude and longitude for at least two pixels, such as a given one and one at an origin. Along with the datum, this will suffice.
The datum, among other things, describes the size and shape of the ellipsoid. We will need its eccentricity e, which is related to the flattening via
e^2 = f * (2 - f).
For WGS 84, for instance,
1 / f = 298.257 223 563
giving e^2 = 0.00669438 and e is approximately 0.08181919084. (For the sphere, f and e are both zero, considerably simplifying the following equations.)
We also need the equatorial radius (semi-major axis) of the ellipsoid, a; for WGS 84, a = 6,378,137.0 meters.
Every image or map has a (global) scale: it is the amount by which the projection coordinates are uniformly multiplied to place projected points (which are in meters) on the images (which is in pixels). Let's call this scale S; its unit of measurement is therefore pixels per meter. Then, relative to an origin on the image, the x (easting) and y (northing) coordinates of the Mercator projection of the point with longitude lambda and latitude phi are
x = S * a * lambda
y = S * a * [log(tan(pi/4 + phi/2)) + (e/2) * (log(1 - e * sin(phi)) - log(1 + e * sin(phi)))]
Although this may look nasty, notice that e and a are known constants, S is an unknown constant, and phi is known for each pixel (presumably, it is the latitude of the pixel's center). Let's encapsulate that information by abbreviating the formula for y as
y = S * a * M(phi; e, a)
(Incidentally, M(phi; 0, a) = log(tan(pi/4 + phi/2)) has a particularly simple form for the sphere.)
Suppose we have a pixel at row i0 and column j0 corresponding to (lambda0, phi0) on the earth, and another pixel at row i1 and column j1 corresponds to (lambda1, phi1). Then we can deduce S in two ways from the two equations. The x equation will work provided the pixels are in different columns and the y equation will work provided the pixels are in different rows:
S = (j1 - j0) / (lambda1 - lambda0) / a
S = (i1 - i0) / (M(phi1; e, a) - M(phi0; e, a)) / a
At least one of these equations will work. With S now known, we can compute the Mercator coordinates for any point (lambda, phi) on the globe directly. We can also invert the projection numerically or by using a (rapidly) iterative algorithm or with a trigonometric series, as described by Snyder on p. 44.
From the information in the question we might guess that S equals 1 pixel / 1 Km = 0.001 pixels per meter. However, that presumes the projection has true scale at the Equator. The allusion to a reference latitude in a comment to the question suggests that the projection has already been rescaled by some amount. That is why it appears necessary to compute S. Otherwise, we would not need information about two pixels: the latitude, longitude, row, and column of a single pixel would suffice.
Reference
Snyder, John P. Map Projections--A Working Manual. USGS professional paper 1395, 1987.
Best Answer
By reading this GDAL ticket http://trac.osgeo.org/gdal/ticket/4977 it should be possible to flip a GeoTIFF by using a negative pixel size in ModelPixelScaleTag. On Y-axis this would mean positive pixel size for GDAL. However, as you can read from the ticket, this trick does not work for you and even if it works you should still place the origin to lower left.
Instead or editing the world file you can probably use .VRT file http://www.gdal.org/gdal_vrttut.html and edit the GeoTransform element e.g.
Even in this case you must also move the origin to lower left in the VRT file. I believe that finally you will decide to write a few rows of Python code for overriding the default GeoTransform and turn the Y pixel size into positive and slide the Y coordinate by (pixel size Y) * (image height).