Convert geodetic coordinates to geocentric (cartesian) I'd like to retain as much float precision as possible without having to resort to specialized libraries so I need this code to work with only standard modules. 

I'm doing this in QGIS where specialized libraries aren't available by default. numpy would be OK to use.  

Sources for the computations can be found on Wikipedia.

Here is my code:

import math

## Ellipsoid Parameters as tuples (semi major axis, inverse flattening)
grs80 = (6378137, 298.257222100882711)
wgs84 = (6378137, 298.257223563)


def geodetic_to_geocentric(ellps, lat, lon, h):
    """
    Compute the Geocentric (Cartesian) Coordinates X, Y, Z
    given the Geodetic Coordinates lat, lon + Ellipsoid Height h
    """
    a, rf = ellps
    lat_rad = math.radians(lat)
    lon_rad = math.radians(lon)
    N = a / math.sqrt(1 - (1 - (1 - 1 / rf) ** 2) * (math.sin(lat_rad)) ** 2)
    X = (N + h) * math.cos(lat_rad) * math.cos(lon_rad)
    Y = (N + h) * math.cos(lat_rad) * math.sin(lon_rad)
    Z = ((1 - 1 / rf) ** 2 * N + h) * math.sin(lat_rad)

    return X, Y, Z


Input:

lat = 43.21009
lon = -78.120123
h = 124
print(geodetic_to_geocentric(wgs84, lat, lon, h))


Output:

(958506.0102730404, -4556367.372558064, 4344627.16166323)


This online converter gives the following output:  

X: 958506.01027304
Y: -4556367.37255806
Z: 4344627.16160147


So, quite close for X and Y but slightly different for Z. Is this a rounding error ? If so, on which side ? Is there a better strategy to do the computations that may retain better precision ?
I can't explain the discrepancy with the online conversion tool without seeing the source code for the latter.

1. Review


The comment for the ellipsoid parameters should give the units for the semi-major axis. Also, I think "reciprocal flattening" would be a better description of the second parameter, given that you use the name rf later on in the code.
The docstring for geodetic_to_geocentric should give the units (degrees) for the latitude and longitude parameters.
"Height above ellipsoid" would be a clearer description for the h parameter.
The parameters to geodetic_to_geocentric could be spelled out in full, which would help understand their meaning. Parameter names are documentation too.
Comments would help understand the meaning of values like N.
The expressions (N + h) * cos(lat), sin(lat), and (1 - 1 / rf) ** 2 are evaluated twice. These could be stored in local variables to avoid duplicated work.


2. Revised code

from math import cos, radians, sin, sqrt

# Ellipsoid parameters: semi major axis in metres, reciprocal flattening.
GRS80 = 6378137, 298.257222100882711
WGS84 = 6378137, 298.257223563

def geodetic_to_geocentric(ellipsoid, latitude, longitude, height):
    """Return geocentric (Cartesian) Coordinates x, y, z corresponding to
    the geodetic coordinates given by latitude and longitude (in
    degrees) and height above ellipsoid. The ellipsoid must be
    specified by a pair (semi-major axis, reciprocal flattening).

    """
    φ = radians(latitude)
    λ = radians(longitude)
    sin_φ = sin(φ)
    a, rf = ellipsoid           # semi-major axis, reciprocal flattening
    e2 = 1 - (1 - 1 / rf) ** 2  # eccentricity squared
    n = a / sqrt(1 - e2 * sin_φ ** 2) # prime vertical radius
    r = (n + height) * cos(φ)   # perpendicular distance from z axis
    x = r * cos(λ)
    y = r * sin(λ)
    z = (n * (1 - e2) + height) * sin_φ
    return x, y, z