# Get the geotransform in Mercator given a pixel size away from the origin

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.

• The information given is not sufficient to solve the problem: you need also to specify the dimensions of each pixel in degrees. Nov 8, 2012 at 13:49
• Well, their program is able to do it, with the pixel size in meters. Actually, with a spherical earth I am able to calculate it, by converting the meters do degrees, but this is "quite" more difficult on a ellipsoid Nov 8, 2012 at 14:00
• Whose program? What additional information does it have? The problem with your question is that "1km x 1 km each pixel" is not a scale. A scale, by definition, is a ratio of distances. To know that ratio, we must know how large a pixel is. BTW, the solution for an ellipsoid is surprisingly simple: I was prepared to post it when I realized you did not provide enough information to carry it out. Nov 8, 2012 at 14:20
• The program that generates the data (Vaisala's IRIS). No extra info I can get :( What they tell is that at the specified coordinates, the pixel is 1km width and 1km height. In a spherical earth, one can get the upper left geo coordinates by simply deltaY = numPixelspixelHeight/earthRadius and deltaX = numPixelspixelHeight/(earthRadius * cos (referenceLat)) (all in radians), and then transform it into mercator. The pixel size can be then calculated by knowing the deltas in one pixel in latlon. Maybe I didn't explain it peoperly, sorry. Nov 8, 2012 at 14:29
• Did you notice that the formula requires the reference latitude? What is it? Could you provide a link to the material you are citing in "What they tell"? Also, what are you really asking concerning "transform[ing] this pixel scale distance into the projection's scale"? What exactly is the "pixel scale distance" and do you really just want to calculate the scale of the projection at each location, or perhaps do you need to compute the parameters of the projection so that the data can be unprojected? Nov 8, 2012 at 14:37

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.

• I have to think better about how to use it, but the information seems excellent, thank you very much! Nov 8, 2012 at 20:07