[GIS] Get the geotransform in Mercator given a pixel size away from the origin

coordinate systemgdalproj

I have a radar image with the following information:

  1. The projection is Mercator (Mercator_1SP in GDAL)
  2. The projection origin is 0º lat 0º lon
  3. At a given geographic coordinates (e.g. lon=10º lat=40º), the scale is 1kmx1km each pixel. I know at wich pixel in the image this coordinates are (they are usualy at the center of the image)
  4. The datum can be WGS84 or others, so it's not always the spherical earth

How can I determine the geotransform? I don't know how to transform this pixel scale distance into the projection's scale at the given point.

–EDIT–

So I need the upper left pixel coordinates and the pixel width and height in the Mercator coordinates from the data explained above.

Thank you.

Best Answer

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.

Related Question