<p>The title of the question reads:</p>
<blockquote>
<p>Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)?</p>
</blockquote>
<p>and the body expands:</p>
<blockquote>
<p><strong>Question:</strong> What exactly are the cos and sin parallax coefficients, and how closely can <span class="math-container">$\arctan2(\sin, \cos)$</span> be expected to serve as a proxy for latitude? Within a degree or so?</p>
</blockquote>
<h2>First the &quot;How closely can (they) serve as a proxy for latitude?&quot; part</h2>
<p><strong>tl;dr:</strong> to within +/-0.2 degrees; the range in the differences between <em>geocentric</em> and <em>geodetic</em> latitude.</p>
<p>The table provides the sin and cos of the same angle <span class="math-container">$\phi'$</span> multiplied by <span class="math-container">$r/R_E$</span> and the 2nd half of the answer will explain that <span class="math-container">$\phi'$</span> is the angle that a radius drawn from the observatory to the Geocenter makes with respect to the Earth's equatorial plane. Everyday latitude or GPS latitude which is a kind of <em>geodetic latitude</em> (what I want to use for my map) is a little bit different because it's drawn on an ellipsoid.</p>
<p>How different?</p>
<p>Well it turns out that <a href="https://space.stackexchange.com/a/27194/12102">I wrote an answer</a> about <em>geodetic latitude</em> about 18 months ago. (@DavidHammen wrote a <a href="https://space.stackexchange.com/a/35654/12102">nicer answer</a> elsewhere as well) Here is an excerpt from my answer:</p>
<blockquote>
<p>The 3D cartesian coordinates <span class="math-container">$X, Y, Z$</span> in Earth-centered, Earth-fixed coordinates assuming an ellipsoidal shape is given by:</p>
<p><span class="math-container">$$X =  \left(N(\phi) + h \right) \cos\phi \cos\lambda $$</span></p>
<p><span class="math-container">$$Y =  \left(N(\phi) + h \right) \cos\phi \sin\lambda $$</span></p>
<p><span class="math-container">$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi $$</span></p>
<p>where <span class="math-container">$\phi, \lambda, h$</span> are latitude, longitude, and altitude, and <span class="math-container">$a, b$</span> are the equatorial and polar radii of the ellipsoid used, and</p>
<p><span class="math-container">$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}. $$</span></p>
<p>Lat, lon, alt in &quot;GPS coordinates} is based on <a href="https://en.wikipedia.org/wiki/World_Geodetic_System#Longitudes_on_WGS_84" rel="nofollow noreferrer">WGS 84</a> with <span class="math-container">$a, b$</span> of 6378.1370 and 6356.7523 kilometers, respectively.</p>
</blockquote>
<p>The parallax coefficients are <span class="math-container">$\rho \cos(\phi')$</span> and <span class="math-container">$\rho \sin(\phi')$</span> and <span class="math-container">$\rho$</span> is the distance from the observatory to the Geocenter divided by &quot;the Earth radius&quot;. That makes them:</p>
<p><span class="math-container">$$\rho \cos(\phi') = \frac{\sqrt{X^2+Y^2}}{R_E}$$</span></p>
<p><span class="math-container">$$\rho \sin(\phi') = \frac{Z}{R_E}$$</span></p>
<p>So if we chose a geodetic latitude <span class="math-container">$\phi$</span>, choose a longitude of zero so that <span class="math-container">$Y=0$</span> and calculate <span class="math-container">$X, Z$</span> and plug them into this:</p>
<p><span class="math-container">$$\phi' = \arctan2\left(Z, X \right) = \arctan2\left(\frac{Z}{R_E}, \frac{X}{R_E} \right) = \arctan2\left( \rho \sin(\phi'), \rho \cos(\phi') \right)$$</span></p>
<p>we can get the geocentric latitude, and then calculate the error by subtracting the two.</p>
<p><strong>It turns out to be about +/- 0.2 degrees with extrema at about +/- 45 degrees.</strong> It is of course zero at the equator and poles.</p>
<p>So to an error of 0.2 degrees, <em>yes that can indeed be used as a proxy!</em></p>
<p><a href="https://i.sstatic.net/NtQeH.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/NtQeH.png" alt="geocentric minus geodetic latitude" /></a></p>
<h2>Now the &quot;What exactly are...&quot; part</h2>
<p>As pointed out <a href="https://astronomy.stackexchange.com/questions/35535/can-i-use-the-parallax-coefficients-for-observatories-as-a-proxy-for-latitude-us#comment69049_35535">in comments</a>, the link explains how the cos and sin parallax coefficients are defined:</p>
<blockquote>
<p>The following list gives the observatory code, longitude (in degrees east of Greenwich) and <strong>the parallax constants (rho cos phi' and rho sin phi', where phi' is the geocentric latitude and rho is the geocentric distance in earth radii</strong>) for each observatory. It is updated nightly.</p>
</blockquote>
<p>What is <a href="https://en.wikipedia.org/wiki/Latitude#Geodetic_and_geocentric_latitudes" rel="nofollow noreferrer">geocentric latitude</a>? It's the &quot;easy one&quot; to understand; it's just the angle measured at the Geocenter between the radius to the point and the closest point on Earth's equator. So while the definition *does not specify a specific Earth radius, it doesn't have to because the table only gives the unit-less parts cos phi' and sin phi'.</p>
<h2>But why both sin and cos? Can't you just use sin between -90 and +90 degrees?</h2>
<p><em>To just define an angle</em> (without the radius), mathematically you might think so. However sin is insensitive near the poles and in the world of real computing numerical errors can crop up. Near the poles sin will still provide the appropriate <em>sign</em> and cos will provide the sensitivity. This way you will have low numerical errors (for this particular step) everywhere.</p>
<p><em>But this is the latitude and distance together</em>, so you have to have two numbers. They could be a latitude angle and a radius explicitly, but this also works. (revised per <a href="https://astronomy.stackexchange.com/questions/35535/can-i-use-the-parallax-coefficients-for-observatories-as-a-proxy-for-latitude-us#comment69135_35537">@DH's comment</a>)</p>
<h2>But what are the parallax coefficients <em>for</em> exaclty? What are they <em>really?</em></h2>
<p>I don't know, but hopefully someone with a valid &quot;astronomer's ID card&quot; will come along and post an answer. I'm to scared to ask a new question on this topic given all the comments under the question! :-)</p>
<hr />
<p>Script for the plot</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
a, b = 6378.1370, 6356.7523 # only their ratio matters in this case
todegs = 180/pi

phi = np.linspace(-halfpi, halfpi, 101)
h = 0.
lam = 0.
N = a**2 / np.sqrt(a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2)
X = (N + h) * np.cos(phi) * np.cos(lam)
Y = (N + h) * np.cos(phi) * np.sin(lam)
Z = ((b/a)**2 * N + h) * np.sin(phi)

phi_prime = np.arctan(Z/X) # two quandrants so no need for arctan2

plt.figure()
plt.plot(todegs*phi, todegs*(phi_prime-phi))
plt.xlabel('geodetic (&quot;GPS&quot;) latitude (degs)', fontsize=12)
plt.ylabel('geocentric - geodetic latitude (degs)', fontsize=12)
plt.xlim(-90, 90)
plt.show()
</code></pre>

CC BY-SA 4.0
2020-03-20T16:25:30.380Z
The parallax coefficients are typically used when applying the topocentric corrections (for a given site on the Earth's surface as opposed to the geocenter). They let you form a 3d vector which can be combined with precession/nutation/frame bias matrix which is what describes the orientation between the celestial coordinate system and the CIO and the terrestrial reference and the TIO. Code examples include [find_orb](https://github.com/Bill-Gray/find_orb/blob/master/mpc_obs.cpp) (search for `rho_cos_phi`) and [SOFA cookbook](http://www.iausofa.org/sofa_pn_c.pdf)
astrosnapper
CC BY-SA 4.0
2020-03-22T03:08:20.010Z
@astrosnapper I've fixed this answer and will read those soon, thank you!
uhoh
CC BY-SA 4.0
2020-03-22T04:32:39.560Z
Re *But why both sin and cos?* That's the wrong question. The given values are $\rho\cos\phi'$ and $\rho\sin\phi'$, not $\cos\phi'$ and $\sin\phi'$ . The given values encode both radial distance and geocentric latitude. The right question is then *why not radial distance and geocentric latitude?* I don't know the answer to that question. An even better question is *why not altitude and latitude?* I don't know the answer to that question, either. (Note well: I used an unqualified "latitude" in my better question. Latitude without a qualifying prefix almost always means geodetic latitude.)
David Hammen
CC BY-SA 4.0
2020-03-23T23:03:55.470Z
@DavidHammen I see what you mean and have revised accordingly
uhoh
CC BY-SA 4.0
2020-03-20T12:37:38.590Z
2020-03-23T23:02:04.777Z
2020-06-17T09:47:40.653Z
uhoh
<p>In <a href="https://astronomy.stackexchange.com/q/35530/7982">How are the three character IAU Minor Planet Center's Observatory Codes assigned? Why are some letters so popular?</a> I looked at the distribution of the first alphanumeric character versus the value of the second two (36 x 100). </p>

<p>The table lists longitudes but no latitudes for these observatories, but using the cos and sin parallax coefficients in <span class="math-container">$\arctan2(\sin, \cos)$</span> and plotting as a function of tabular longitude produces a map that looks like it could be the distribution of telescopes on Earth.</p>

<p><strong>Question:</strong> What exactly are the cos and sin parallax coefficients, and how closely can <span class="math-container">$\arctan2(\sin, \cos)$</span> be expected to serve as a proxy for latitude? Within a degree or so?</p>

<p><a href="https://i.sstatic.net/3OvZi.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/3OvZi.png" alt="arctan2(sin, cos) as proxy for latitude"></a></p>

<pre><code>import numpy as np
import matplotlib.pyplot as plt
# blob is the unformatted lines from https://www.minorplanetcenter.net/iau/lists/ObsCodes.html
threes = [line[:3] for line in blob.splitlines()]
key = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
a = [[x] for x in key]
for t in threes:
    a[key.find(t[0])].append(t)
b = np.zeros((36, 100))
for i, thing in enumerate(a):
    c = [int(x[1:]) for x in thing[1:]]
    for d in c:
        b[i, d] = 1
plt.imshow(-b, interpolation='nearest', cmap='gray')
plt.gca().set_aspect(2)
plt.show()
pairs = []
todegs = 180/np.pi
for line in blob.splitlines():
    try:
        cos, sin = [float(x) for x in (line[13:21], line[21:30])]
        lat = todegs * np.arctan2(sin, cos)
        lon = float(line[4:13])
        pairs.append([lon, lat])
    except:
        pass
print(len(pairs), len(blob.splitlines()), float(len(pairs))/len(blob.splitlines()) )
lon, lat = np.array(list(zip(*pairs)))
lon = np.mod(lon+180, 360) - 180
plt.figure()
plt.plot(lon, lat, 'ok', markersize=1)
plt.ylim(-90, 90)
plt.xlim(-180, 180)
plt.gca().set_aspect('equal')
plt.show()
</code></pre>

CC BY-SA 4.0
2020-03-20T10:07:16.800Z
@StephenG -- Right here on this corner of the SE network, there's the question [What is ′ in orbital mechanics?](https://astronomy.stackexchange.com/questions/19930)
David Hammen
CC BY-SA 4.0
2020-03-22T03:09:13.277Z
@DavidHammen I've posted an answer. I think I've got it right, and the answer is about +/- 0.2 degrees. How does it look?
uhoh
<observational-astronomy><coordinate><observatory><parallax>
Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)?
uhoh
# Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)?

(Asked by: username_0 on 2020-03-20)

In How are the three character IAU Minor Planet Center's Observatory Codes assigned? Why are some letters so popular? I looked at the distribution of the first alphanumeric character versus the value of the second two (36 x 100). 
The table lists longitudes but no latitudes for these observatories, but using the cos and sin parallax coefficients in $\arctan2(\sin, \cos)$ and plotting as a function of tabular longitude produces a map that looks like it could be the distribution of telescopes on Earth.
Question: What exactly are the cos and sin parallax coefficients, and how closely can $\arctan2(\sin, \cos)$ be expected to serve as a proxy for latitude? Within a degree or so?

import numpy as np
import matplotlib.pyplot as plt
# blob is the unformatted lines from https://www.minorplanetcenter.net/iau/lists/ObsCodes.html
threes = [line[:3] for line in blob.splitlines()]
key = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
a = [[x] for x in key]
for t in threes:
    a[key.find(t[0])].append(t)
b = np.zeros((36, 100))
for i, thing in enumerate(a):
    c = [int(x[1:]) for x in thing[1:]]
    for d in c:
        b[i, d] = 1
plt.imshow(-b, interpolation='nearest', cmap='gray')
plt.gca().set_aspect(2)
plt.show()
pairs = []
todegs = 180/np.pi
for line in blob.splitlines():
    try:
        cos, sin = [float(x) for x in (line[13:21], line[21:30])]
        lat = todegs * np.arctan2(sin, cos)
        lon = float(line[4:13])
        pairs.append([lon, lat])
    except:
        pass
print(len(pairs), len(blob.splitlines()), float(len(pairs))/len(blob.splitlines()) )
lon, lat = np.array(list(zip(*pairs)))
lon = np.mod(lon+180, 360) - 180
plt.figure()
plt.plot(lon, lat, 'ok', markersize=1)
plt.ylim(-90, 90)
plt.xlim(-180, 180)
plt.gca().set_aspect('equal')
plt.show()

**Tags:** observational-astronomy, coordinate, observatory, parallax
> **Comments:**
> **username_1**: @StephenG -- Right here on this corner of the SE network, there's the question [What is ′ in orbital mechanics?](https://astronomy.stackexchange.com/questions/19930)  
> **username_0**: @DavidHammen I've posted an answer. I think I've got it right, and the answer is about +/- 0.2 degrees. How does it look?  

---

## 1 answers

### Answer by username_0 on 2020-03-20. Score: 0 | Accepted answer

The title of the question reads:

Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)?

and the body expands:

Question: What exactly are the cos and sin parallax coefficients, and how closely can $\arctan2(\sin, \cos)$ be expected to serve as a proxy for latitude? Within a degree or so?

First the "How closely can (they) serve as a proxy for latitude?" part
tl;dr: to within +/-0.2 degrees; the range in the differences between geocentric and geodetic latitude.
The table provides the sin and cos of the same angle $\phi'$ multiplied by $r/R_E$ and the 2nd half of the answer will explain that $\phi'$ is the angle that a radius drawn from the observatory to the Geocenter makes with respect to the Earth's equatorial plane. Everyday latitude or GPS latitude which is a kind of geodetic latitude (what I want to use for my map) is a little bit different because it's drawn on an ellipsoid.
How different?
Well it turns out that I wrote an answer about geodetic latitude about 18 months ago. (@DavidHammen wrote a nicer answer elsewhere as well) Here is an excerpt from my answer:

The 3D cartesian coordinates $X, Y, Z$ in Earth-centered, Earth-fixed coordinates assuming an ellipsoidal shape is given by:
$$X =  \left(N(\phi) + h \right) \cos\phi \cos\lambda $$
$$Y =  \left(N(\phi) + h \right) \cos\phi \sin\lambda $$
$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi $$
where $\phi, \lambda, h$ are latitude, longitude, and altitude, and $a, b$ are the equatorial and polar radii of the ellipsoid used, and
$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}. $$
Lat, lon, alt in "GPS coordinates} is based on WGS 84 with $a, b$ of 6378.1370 and 6356.7523 kilometers, respectively.

The parallax coefficients are $\rho \cos(\phi')$ and $\rho \sin(\phi')$ and $\rho$ is the distance from the observatory to the Geocenter divided by "the Earth radius". That makes them:
$$\rho \cos(\phi') = \frac{\sqrt{X^2+Y^2}}{R_E}$$
$$\rho \sin(\phi') = \frac{Z}{R_E}$$
So if we chose a geodetic latitude $\phi$, choose a longitude of zero so that $Y=0$ and calculate $X, Z$ and plug them into this:
$$\phi' = \arctan2\left(Z, X \right) = \arctan2\left(\frac{Z}{R_E}, \frac{X}{R_E} \right) = \arctan2\left( \rho \sin(\phi'), \rho \cos(\phi') \right)$$
we can get the geocentric latitude, and then calculate the error by subtracting the two.
It turns out to be about +/- 0.2 degrees with extrema at about +/- 45 degrees. It is of course zero at the equator and poles.
So to an error of 0.2 degrees, yes that can indeed be used as a proxy!

Now the "What exactly are..." part
As pointed out in comments, the link explains how the cos and sin parallax coefficients are defined:

The following list gives the observatory code, longitude (in degrees east of Greenwich) and the parallax constants (rho cos phi' and rho sin phi', where phi' is the geocentric latitude and rho is the geocentric distance in earth radii) for each observatory. It is updated nightly.

What is geocentric latitude? It's the "easy one" to understand; it's just the angle measured at the Geocenter between the radius to the point and the closest point on Earth's equator. So while the definition *does not specify a specific Earth radius, it doesn't have to because the table only gives the unit-less parts cos phi' and sin phi'.
But why both sin and cos? Can't you just use sin between -90 and +90 degrees?
To just define an angle (without the radius), mathematically you might think so. However sin is insensitive near the poles and in the world of real computing numerical errors can crop up. Near the poles sin will still provide the appropriate sign and cos will provide the sensitivity. This way you will have low numerical errors (for this particular step) everywhere.
But this is the latitude and distance together, so you have to have two numbers. They could be a latitude angle and a radius explicitly, but this also works. (revised per @DH's comment)
But what are the parallax coefficients for exaclty? What are they really?
I don't know, but hopefully someone with a valid "astronomer's ID card" will come along and post an answer. I'm to scared to ask a new question on this topic given all the comments under the question! :-)

Script for the plot
import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
a, b = 6378.1370, 6356.7523 # only their ratio matters in this case
todegs = 180/pi

phi = np.linspace(-halfpi, halfpi, 101)
h = 0.
lam = 0.
N = a**2 / np.sqrt(a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2)
X = (N + h) * np.cos(phi) * np.cos(lam)
Y = (N + h) * np.cos(phi) * np.sin(lam)
Z = ((b/a)**2 * N + h) * np.sin(phi)

phi_prime = np.arctan(Z/X) # two quandrants so no need for arctan2

plt.figure()
plt.plot(todegs*phi, todegs*(phi_prime-phi))
plt.xlabel('geodetic ("GPS") latitude (degs)', fontsize=12)
plt.ylabel('geocentric - geodetic latitude (degs)', fontsize=12)
plt.xlim(-90, 90)
plt.show()

> **Comments:**
> **username_2**: The parallax coefficients are typically used when applying the topocentric corrections (for a given site on the Earth's surface as opposed to the geocenter). They let you form a 3d vector which can be combined with precession/nutation/frame bias matrix which is what describes the orientation between the celestial coordinate system and the CIO and the terrestrial reference and the TIO. Code examples include [find_orb](https://github.com/Bill-Gray/find_orb/blob/master/mpc_obs.cpp) (search for `rho_cos_phi`) and [SOFA cookbook](http://www.iausofa.org/sofa_pn_c.pdf)  
> **username_0**: @astrosnapper I've fixed this answer and will read those soon, thank you!  
> **username_1**: Re *But why both sin and cos?* That's the wrong question. The given values are $\rho\cos\phi'$ and $\rho\sin\phi'$, not $\cos\phi'$ and $\sin\phi'$ . The given values encode both radial distance and geocentric latitude. The right question is then *why not radial distance and geocentric latitude?* I don't know the answer to that question. An even better question is *why not altitude and latitude?* I don't know the answer to that question, either. (Note well: I used an unqualified "latitude" in my better question. Latitude without a qualifying prefix almost always means geodetic latitude.)  
> **username_0**: @DavidHammen I see what you mean and have revised accordingly  

---

