# Lunar Reconnaissance Orbiter (LRO) stereographic projection transformation equations issues

I have been trying to figure out the coordinate frames the Lunar Orbiter Laser Altimeter (LOLA) team uses for data.

In their dsmap_polar.cat file they state: \begin{align} X&=R \sin{\lambda}\,,\\ Y&=R \cos{\lambda}\,, \end{align} where $$\lambda$$ is the longitude and

$$R=2*1737400*\tan{(90-\theta_{\rm L})}\,,$$ $$\theta_{\rm L}$$ being the latitude.

The issue here is the meaning of their "longitude". The standard spherical system is defined:

\begin{align} x&= r \sin{\theta}\cos{\phi}\,,\\ y&= r\sin{\theta}\sin{\phi}\,,\\ z&= r\cos{\theta}\,, \end{align} where the angles are defined in the figure below. The standard chosen by LOLA is unusual as it seems to be rotated (they use the work by Snyder I think). A significant issue arises when e.g. one wants to convert to gnomonic map. When I use the above standard spherical system, i do not get back to the same point from the gnomonic map.

I am looking for the motivation and explanation why would someone choose sine of longitude for x coordinate and cosine for y coordinate so that I can better understand what i have to change in my transformation equations.

Also, if one transforms from stereographic projection latitude longitude system as per their equations, does one really end up with x-axis pointing to mean-Earth direction?

## Detailed background

What i really need is to find the error in my conversion formulas. I need to be able to convert between LOLA stereographic projection and gnomonic projection.

Gnomonic to cartesian planetocentric system (C1): \begin{align} x &= \frac{\hat{x}}{\sqrt{ 1 + \left(\frac{\hat{x}}{a}\right)^2 + \left(\frac{\hat{y}}{a}\right)^2 }}\,,\\ y &= \frac{\hat{y}}{\sqrt{ 1 + \left(\frac{\hat{x}}{a}\right)^2 + \left(\frac{\hat{y}}{a}\right)^2 }}\,,\\ z &= \frac{a}{\sqrt{ 1 + \left(\frac{\hat{x}}{a}\right)^2 + \left(\frac{\hat{y}}{a}\right)^2 }}\,, \end{align} where hat values are the coordinates in gnomonic projection.

Cartesian planetocentric to cartesian stereographic projection (C2): \begin{align} \tilde{x} &= \frac{2\,a\,y}{a + z}\,,\\ \tilde{y} &= \frac{2\,a\,x}{a + z}\,, \end{align} where the tilde values are cartesian coordinates in stereographic projection and $$a$$ is the reference lunar radius. The switched $$x$$ and $$y$$ in the numerator should reflect the LOLA team strange choice.

Now the inverse transformation are given by:

Stereographic cartesian system to latitude longitude system (provided in the LOLA catalogue file for the norther hemisphere) (C3): \begin{align} \lambda &= \arctan{\frac{\tilde{y}}{\tilde{x}}}\,,\\ \theta_{\rm L} &= 90 - 2\arctan{\frac{r}{2a}}\,, \end{align} where $$r = \sqrt{\tilde{x}^2 + \tilde{y}^2}\,.$$

From latitude and logitude to gnomonic cartesian system (C4): \begin{align} \hat{x} &= a \tan(90-\theta_{\rm L})\cos{\lambda}\,,\\ \hat{y} &= a \tan(90-\theta_{\rm L})\sin{\lambda}\,, \end{align}

From here, latitude and longitude to cartesian planetocentric $$(x,y,z)$$ is converted using the standard above.

Now, if i compare the results of the transformation (C1) and the (C3) converted to cartesian system by the standard above, the x and y axes are switched. The problem is, i am not sure which transformation is correct. What i really want is to be able to get to the planetocentric coordinates with x axis as mean earth direction.

## Specific numerical example

The code below uses Matlab syntax. $$a$$ stands for the reference radius of the lunar sphere. Let's start with some specific point in gnomonic projection:

a   =  1737400;
x_g = -411141.107140;
y_g = -411325.894877;


Now, convert this point from gnomonic coordinates to cartesian planetocentric which is derived using the standard spherical coordinates above:

x = x_g/sqrt( 1 + (x_g/a)^2 + (y_g/a)^2 );
y = y_g/sqrt( 1 + (x_g/a)^2 + (y_g/a)^2 );
z =   a/sqrt( 1 + (x_g/a)^2 + (y_g/a)^2 );


This gives:

(x,y,z) = (-389878.216487,-390053.447599,1647547.281362)


I continue to convert to stereographic projection using (note the interchanges $$x$$ and $$y$$ in the numerator on the right hand side which should bring one to the system of the LOLA team - hopefully):

x_s = 2*a*y/(a + z);
y_s = 2*a*x/(a + z);


which gives:

(x_s,y_s) = (-400407.334903,-400227.452318)


To check for consistency, I do the inverse transformation to make sure that i get the same point in gnomonic system. Taking this (x_s,y_s) point as a starting point, I convert it back to latitude,logitude system using official formulas of the LOLA team from the catalogue file above:

r   = sqrt(x_s^2 + y_s^2);
lon = atan2d(y_s,x_s);
lat = 90 - 2*atand(0.5*r/a);


and this results:

(lat,lon) = (71.492702,-135.012873)


Now, converting into cartesian system using:

theta = 90-lat;
phi   = lon;
x = a*sind(theta)*cosd(phi);
y = a*sind(theta)*sind(phi);
z = a*cosd(theta);


results in:

(x,y,z) = (-390053.447599,-389878.216487,1647547.281362)


You can now see that compared to the previous result that i got when converting from the gnomonic, here the x and y axes are interchanged here already. So i do not continue with the final transformations to gnomonic system. Also note that this is one specific point, so i am not sure, if generally, i only get interchanged $$x$$ and $$y$$ axes or for other points, other issues arise such as different signs.

• What is a" LOLA team"? It's unhelpful to use specific acronyms that many readers will not recognize, and it's always good practice to write out the phrase in full the first time before using its acronym later. Can you update your post to spell it out and maybe add a link? Thanks!
– uhoh
Jun 5, 2021 at 3:49
• I think we need a numeric example to see and understand the kind of error you get.
– Uwe
Jun 5, 2021 at 15:20
• @Uwe I did include a numerical example including Matlab code to reproduce the issue. Jun 5, 2021 at 16:14

I tried to do the same using python:

import numpy as np

a   =  1737400
x_g = -411141.107140
y_g = -411325.894877

x = x_g/np.sqrt( 1 + np.square(x_g/a) + np.square(y_g/a) )
y = y_g/np.sqrt( 1 + np.square(x_g/a) + np.square(y_g/a) )
z =   a/np.sqrt( 1 + np.square(x_g/a) + np.square(y_g/a) )

print('(x, y, z) = ', (x, y, z))

x_s = 2*a*y/(a + z)
y_s = 2*a*x/(a + z)

print('(x_s, y_s) = ', (x_s, y_s))

r   = np.sqrt(np.square(x_s) + np.square(y_s));
lon = np.arctan2(y_s,x_s);

print('(lat, lon) = ', (lat, lon))

x_d = a * np.tan(np.deg2rad(90.0)-lat) * np.sin(lon)
y_d = a * np.tan(np.deg2rad(90.0)-lat) * np.cos(lon)

print('(x_d, y_d) = ', (x_d, y_d))


result :

(x, y, z) =  (-389878.2164872003, -390053.4475990476, 1647547.281362467)
(x_s, y_s) =  (-400407.3349029025, -400227.4523177882)
(lat, lon) =  (1.247783033129123, -2.356419165143969)
(x_d, y_d) =  (-411141.10714000027, -411325.8948770001)


Compare your equations above to :

X = RSIN(LONPI/180)

Y = RCOS(LONPI/180)

You used cos for x and sin for y, but in LOLA sin is used for x and cos for y.

x and y values are not switched now and

x_g = -411141.107140

y_g = -411325.894877

(x_d, y_d) = (-411141.10714000027, -411325.8948770001)

are equal now.

So I think I finally found the error. I made some additional errors when I searched for your error. Hopefully there are no more errors in my answer.

• In matlab i am using sind et al (whose argument is in degrees, not in radians) and I believe my code is consistenly using functions that accept degrees. If you look at you (x,y,z) results, they also have switched x and y coordinates. I am quite sure, this is not the issue. I believe the issue is really that unusual choice of from where to measure $\phi$ or azimuthal angle. Jun 5, 2021 at 19:53
• Ok, i will look into it but this way i get the same "messed up" coordinate system in gnomonic projection, that is, azimuth taken from y axis. I guess the same solution would be achieved if one calculated the lat,lon system to spherical coordinates with sine and cosine of azimuth interchanged. Also there is still the issue what points to the mean-Earth direction in the planetocentric frame, is it x or y? Actually, this is the crucial issue because this gives the orientation in space. Jun 6, 2021 at 1:31
• This is what coordinate specification for LRO say(lunar.gsfc.nasa.gov/library/451-SCI-000958.pdf): "These are right-handed spherical coordinates where the z-axis is the mean axis of rotation and the x-axis is the intersection of the Equator and the Prime Meridian. ... The y-axis is orthogonal to the x- and z-axes. ... Longitude is the angle between this vector and the plane of the Prime Meridian measured in an eastern direction." and in part 6.1 "The Prime Meridian (0 ̊ Longitude) is defined by the mean Earth direction" Jun 6, 2021 at 1:40
• Based on the above, x axis in planetocentric, cartesian system points to the mean earth direction, but based on their transformation equations, it is the y axis. Or am i wrong? (because the azimuth is with respect to the y axis) Jun 6, 2021 at 1:41
• @leosenko "Based on the above, x axis in planetocentric, cartesian system points to the mean earth direction, but based on their transformation equations, it is the y axis." May be this problem is a result of your switching between x and y axis by switching sin and cos?
– Uwe
Jun 7, 2021 at 8:17