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.
