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'$ and the 2nd half of the answer will explain that that 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.
So if we chose a geodetic latitude $\phi$, calculate choose a longitude of zero so that $Y=0$ and calculate $X, Z$ and plug them into this:
$$\phi' = \arcsin\left(\frac{Z}{\sqrt{X^2+Z^2}} \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.
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?
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 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.arcsin(Z/np.sqrt(X**2 + Z**2))
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()