The code you cited as reference assumes h is ellipsoidal height, which is pretty standard when you specify WGS84 as your reference frame (I would avoid using the word Datum if not talking about a legacy local reference frame).
When you say:
x = (h + N) * cos_lambda * cos_phi;
y = (h + N) * cos_lambda * sin_phi;
z = (h + (1 — e_sq) * N) * sin_lambda;
You are using standard Lat+long+ellipsoid height to ECEF formulae.
Saying "above sea level" is often misleading, because sea-level is not exactly a equipotential surface of the earth. In case you have "above sea level" heights, what you really have are orthometric heights, your data provider must specify which and where the reference surface (often called W0) is. This is almost surely NOT your case.
On the other hand, your goal is ECEF to local tangent plane coordinate conversion. Once you know that the heights are ellipsoidal, the transformation is simple roto-translation, since both systems are cartesian coordinates.
Please Refer to This wikipedia page: https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_ECEF_to_ENU.
If heights were not ellipsoidal, read below for discussion, but you could:
- Make the conversion of heights using a geoid model prior to transformation. You will get coordinates on a plane tangent to the ellipsoid.
- If the work area is small, fairly flat and near sea level, you can do the math for horizontal East and North coordinates using
h = 0, and map the original heights to the Up axis. On small areas the biases introduced by the 0 height imposition would be unnoticeable. Here "small" depends on the scale of your final representation.
Think before going on that your math doesn't need to be more precise than your final product needs. Keep in mind that differences between ellipsoidal and orthometric heights (called geoid undulation) are up to 10-50 meters, but variations are smooth with geography, and may be ok for you if you are looking for relative between-points precision instead of absolute values.
More discussion on Heights
TL;DR
- Sea surface is not a surface of equal level. There is some evil thing called "dynamic sea topography" preventing it
- Any "above sea level" measurement means really "above a surface level that on this specific tide gauge is equal to the average or all measures we have of sea level"
- We should probably stop calling those "above sea level" heights and use the more accurate "orthometric heights".
- If the data provider didn't mention, heights are probably above ellipsoid or useless, may be both
- There is no single "radius of the earth".
More on "above sea level"
On old cartography, heights were measured as "above sea level" heights. This is a logical choice because intuition says that the sea is big enough to allow its surface to adopt the shape of an equipotential surface. It would be certainly true if the only forces driving the water were static forces, i.e. gravity, friction, etc. But there are certain dynamics of currents, heat transport, and ocean atmosphere interaction, that prevents the surface of the seas from conforming to a level surface. (No mention to the tides, which can theoretically be filtered out with time averaging)
Prior to the proliferation of artificial satellites, there where no other observable "vertical" line than the plumb line. The plumb line is driven by gravity, and is always perpendicular to equipotential surfaces of the earth. The origin of heights were usually defined by the mean sea level at some tide gauge near to the study region. From the sea to the mainland, the reference surface was propagated with a combination of methods including gravity measurements, triangulation networks, geodetic calculations, and geometric leveling. There are many interesting works from the 1950s on international geoids and leveling conventions.
With the beginning of the satellite era, the path to 3D measurements was open. We only needed accurate orbit and signal propagation models; and both things were needed anyway for maintenance tasks. 3D measurements leads to real geocentric reference frames, ECEF (Earth-Centered Earth-Fixed). Then the measurement of heights experienced a shift, or a bifurcation. Any aircraft-generated data will necessarily measure positions from the ECEF reference frame, and any topographic measure will result in heights relative to an equipotential surface.
Final note: Connecting both heights
For connecting both types of heights, geophysics and geodesy works on geoid models, both local and global. I refer the reader to any "physical geodesy" book, i.e. Hoffman and Moritz. Or you may want to visit this link and explore the models and its associated papers: http://icgem.gfz-potsdam.de/tom_longtime
It is useful to know this, because some global height models uses global geoid models to provide orthometric heights instead of the measured ellipsoidal heights (case: SRTM).
I admit I may have diverted from the original question. Time will tell.
Edit 2: on the radius of the (simplified model of the) earth.
Also, there's a mention to the earth radius on your first paragraph. It is worth noting that the radius of the ellipsoid is not unique nor constant. In fact, for each point on the surface of the ellipsoid there are two main radius often called M, the meridional radius, and N, the prime-vertical radius. N is closely related to the radius of a parallel, and is the same N in the code above.
It's not always noted, but the far above cited formulae work in part due to the fact that N, the prime vertical radius, is the length of the segment, in the direction perpendicular to the ellipsoid, from the surface to the intersection with the Z axis, the axis of rotation.
You may have also read about major and minor radius which refers to the parameters of the ellipsoid:
a, the equatorial radius, often called semi-major axisb, the polar radius, also called semi-minor axis
