Radius of the earth for a particular location

[izzatzubir]

Skyfield provide a neat positioning, using WGS84.

However, WGS84 considers the earth as an oblate spheroid, thus the radius of the earth at different latitude/longitude are different. There are formula to calculate the radius, using geodetic latitude, polar radius and equatorial radius, such as here: https://planetcalc.com/7721/

I would like to request to call this method, given the position of the location. This will be useful to calculate arc of vision (ARCV) which is historically used to calculate lunar crescent sighting.

Thanks.

[brandon-rhodes]

Do you want the radius from the Earth's center to the latitude/longitude point, or do you want the length of the vector in the "down" / "nadir" direction to where it intersects the Earth's equatorial plane? Maybe a diagram will help. From the Wikipedia:

https://en.wikipedia.org/wiki/Geodetic_coordinates

Are you wanting the actual radius — the length of a line from the center O to the observer position N? Or are you wanting the (shorter) length of the line CN?

[izzatzubir]

It is from center of the earth, O. The parameter that I would like to find is the 'geocentric altitude' of the moon. The topocentric altitude can be found using (earth + location).at(t).observe(moon).apparent().altaz(). The geocentric altitude would be the altitude with respect to a 'translated horizon', where the new frame of reference shares the same z-axis (along zenith-nadir line), while the x,y frame has been translated along the z-axis by the radius of the earth.

Here are some codes that I have experimented:

from mpmath import  sqrt
from skyfield import api
from skyfield.api import wgs84, N, E, load

eph = api.load('de440s.bsp')
earth, moon = eph['earth'], eph['moon']
location = wgs84.latlon(5.411, 100.2) #set the lat,lon 

for i in range (5):
    ts = load.timescale().utc(2023,6,10,12,i) #test a date
    vector = location.at(ts).xyz.km #find the vector of the location. Assuming the origin of the vector is the center of the earth
    #print(vector) #1
    x,y,z = vector
    radius = sqrt(x**2 + y**2 + z**2) #find the magnitude of the vector, by right, should be the radius of the earth at the location, in km.
    revert_to_center = wgs84.latlon(5.411, 100.2, -radius*1000) #choose the same location, but elevation set to the negative radius above. By right, this shoul be at the center of the earth. The frame of reference shares the same zenith.
    observe_from_earth = earth.at(ts).observe(moon).apparent().distance().km #observe from earth vector (399)
    observe_from_reverted_radius = (earth+revert_to_center).at(ts).observe(moon).apparent().distance().km #observe from the 'de-elevated' location.

    print(abs(observe_from_earth-observe_from_reverted_radius)) #find the difference. Should be zero if assumptions are correct. 

[izzatzubir]

New information that I just found: The wgs84 origin point is at Earth's center of mass. I assume it is not necessarily at the center of the 'oblate spheroid' earth, but could be shifted. Am I wrong?

If this is the case, does the ICRS earth vector, actually points to this origin or no?

[brandon-rhodes]

I have always seen the center of the Earth's mass and the center of the oblate spheroid treated as the same exact point.

Your calculation sqrt(x**2 + y**2 + z**2) is, I believe, exactly correct, for producing the value you originally asked for: the distance from the center of the Earth to a given geographic location. So I think you have now solved your original question!

The reason your "revert to center" position doesn't exactly align with the center of the Earth can be seen if you look back at the diagram above: a negative elevation moves the position down along the C-N line, not the O-N line, so you miss the center of the Earth and instead wind up a little bit into the spheroid's southern hemisphere.

I'm going to go ahead and close this, since you worked out the formula you needed, but please feel free to comment further!