There are quite a few situations where you want to convert latitude / longitude to a position in a Cartesian coordinate system, or to get a straight line distance from latitude / longitude ** (not).

In fact, I think many other people have written this kind of article. However, I think that the coordinate system you want to use may differ depending on the application. So, I learned how to convert to ** a three-dimensional Cartesian coordinate system that does not move with the center of gravity of the earth as the origin **, so I would like to summarize it here with the meaning of a memorandum.

It's a little about, so if you have any suggestions along the way, please.

About the representation of latitude and longitude

First of all, I would like to touch on the representation of latitude and longitude as a premise. The following two types of expressions are mainly used to describe latitude and longitude. Be careful not to confuse it.

Notation in degrees, minutes and seconds

A format such as 135 ° 12'34.56" in sexagesimal. It's the one you often hear. It is read as "135 degrees 12 minutes 34 seconds 56".

It may be written as 135 ° 12'34" 56 or read as "135 ° 12 minutes 34.56 seconds". Then it seems that it may be written as 1351234.56 or 135.12.34.56. Personally, it's a word of ** please stop because it's really confusing **, but you can distinguish it from the following decimal notation by the way the period is attached.

Decimal notation

As the name implies, this is a decimal notation, such as 135.2096. The "minutes" and "seconds" are represented by a small number of "degrees". The 135 ° 12'34.56" in the above example is

135 + (12 / 60) + (34.56 / 3600)

Is calculated and converted to the decimal notation 135.2096. The reading is "135.2096 degrees".

Some degrees, minutes, and seconds notation looks very similar to decimal notation, but you can see that ** you can tell which one you use by looking at the way the decimal point is added **.

About the coordinates to be obtained here

This article explains how to convert latitude and longitude information measured by GPS into a position in an absolute three-dimensional coordinate system. The important point is that it explains ** just the procedure for getting 3D coordinates **.

In other words, what I mean is ** it doesn't matter if the earth is actually there **. Even if the coordinates of two arbitrary points are obtained by this method, the distance between those two points is ** a straight line distance through the earth **. That's why the method in this article is useful when thinking about airplanes and spacecraft (because it doesn't move along the surface of the earth).

When thinking about buildings and automobiles on the ground, for example, ** a method to obtain two-dimensional coordinates of a certain point with an arbitrary point as the origin (and consider the aspect of the earth) ** is useful. This article is very helpful.

Implementing mutual conversion between latitude / longitude and plane rectangular coordinates in Python --Qiita

What is the WGS-84 coordinate system?

The calculation used this time is a three-dimensional Cartesian coordinate system called ** WGS-84 ** used by GPS satellites. It is a kind of Earth fixed Cartesian coordinate system whose origin is the center of gravity of the earth. It means that ** it spins around with the earth **. The parameters of the WGS-84 reference ellipsoid are shown below ... I thought so, but if necessary, please google each one.

How to ask

First, place the variables as shown in the figure below. People who come up with this kind of thing are amazing. Latitude φ, Longitude λ, Ellipse height h, Coordinates of point Px, y, z. The ellipsoid in the center is the earth. GPS uses WGS-84 for this ellipsoid.

** Edited on 2020.04.01 ** The English translations of ellipsoidal height and altitude may fluctuate, so I deleted them. If you have any knowledge, I would appreciate it if you could comment.

Supplement: What is elliptical body height?

As the name suggests, it is the altitude ** measured from the surface of an ellipsoid that looks like the earth. The altitude returned by GPS is the height of this ellipse. It should be noted here that the elliptical body height is often referred to as ** different from the altitude **. Altitude is ** altitude measured from the average sea level. The line that extends the average sea level to land is called a geoid, and ** altitude ** is the distance between the target and the geoid. And the distance between the ellipsoid and the geoid is called ** geoid height **. In other words, ** elliptical body height can be expressed by the sum of altitude and geoid height **. It's complicated. When I illustrated it, it became a tremendous color scheme. The geoid height varies depending on the point on the earth and can be obtained from the calculation site of the Japan Meteorological Agency. https://vldb.gsi.go.jp/sokuchi/surveycalc/geoid/calcgh/calc_f.html

ceremony

Since it is difficult to write the entire derivation process, the coordinates of the point P can be expressed as follows by ** folding all edges **. Here, ʻa is the semimajor axis, ʻe is the eccentricity, and in the case of a WGS-84 reference ellipsoid, the following constants are used (e is actually calculated from the flattening).

a = 6378137 [m]
e = 0.0818191908426 [-]

The above formula seems to be ** famous (?) ** among those on the road, and if you look for it, you will find various derivations. Please google if necessary. Don't use this article as a reference.

When implemented in Python

I decided to make it a function like this. Arguments: latitude φ [deg], longitude λ [deg], ellipse height h [m] Return value: x coordinate [m], y coordinate [m], z coordinate [m]

import numpy as np

def xyz(phi, lamb, h):  #Latitude and longitude of the argument are in decimal notation
    
    #Convert latitude and longitude to radians
    phi = np.deg2rad(phi)
    lamb = np.deg2rad(lamb)

    #constant(WGS-For 84 reference ellipsoids)
    a = 6378137          #Semimajor axis[m]
    e = 0.0818191908426  #Eccentricity[-]
    
    N = a / np.sqrt(1 - (e ** 2 * np.sin(phi) ** 2))
    x = (N + h) * np.cos(phi) * np.cos(lamb)
    y = (N + h) * np.cos(phi) * np.sin(lamb)
    z = (N * (1 - e ** 2) + h) * np.sin(phi)
    
    return x, y, z  #[m]

that's all. please excuse my poor writing.

[PYTHON] Convert latitude, longitude, GPS altitude to 3D Cartesian coordinates