Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)? 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()
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 &quot;How closely can (they) serve as a proxy for latitude?&quot; 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 &quot;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 &quot;the Earth radius&quot;. 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 &quot;What exactly are...&quot; 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 &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'.
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 &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! :-)

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 (&quot;GPS&quot;) latitude (degs)', fontsize=12)
plt.ylabel('geocentric - geodetic latitude (degs)', fontsize=12)
plt.xlim(-90, 90)
plt.show()