[GIS] How to convert longitude & latitude to a 3d (geocentric x,y,z) coordinate

centroidscoordinate systemqgisr

I'm interested in finding the centroid of a point cloud. I understand that to do this I need to convert latitude and longitude values to 3d coordinates, average them and then convert the resulting coordinate back to latitude and longitude (From a previous answer: https://gis.stackexchange.com/a/6026/8891).

How would I go about doing this (preferably in either QGIS or R)? Is there an appropriate formula which I can recreate in R or a function already available in QGIS?

The main aim of this is to identify the following: If each point in the point cloud was a person, and all the people had to travel (as the crow flies) to a single point, which point would minimise the total distance travelled by the whole group?

Best Answer

if you can use python, you can convert coordinates from WGS84 to geocentric coordinate system (EPSG:4978) using this snippet of python code.

import pyproj

wgs84 = pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')##4326
geocentric= pyproj.Proj('+proj=geocent +datum=WGS84 +units=m +no_defs') ##4978

x2 , y2  , z2= pyproj.transform( wgs84, geocentric, x, y, z)

Related Solutions

Local Coordinate to Geocentric Conversion – How to

The question asks how to convert local coordinates into geocentric coordinates. As an example,

What would be the geocentric coordinates of a point displaced 250 meters west and 250 meters south of -108.619987 degrees longitude, 36.234600 degrees latitude (at sea level)?

(For the reason why only six decimal places are used here, please see https://gis.stackexchange.com/a/8674.)

This encompasses two geometrical questions:

What are the geocentric coordinates corresponding to a given location given in terms of latitude and longitude (at a nominal height of zero)?
What rotation will bring the local frame (in which x points east, y points north, and z points up) into the geocentric frame (in which x points towards the Prime Meridian at the Equator, z points at the North Pole, and x-y-z form a positively oriented orthogonal coordinate system)?

Let's solve these geometrically. Because the math for the first one has been done elsewhere on this site, the new part is the second.

To begin with, consider the geocentric frame at the North Pole:

Frame 1

In this and subsequent images, the blue arrow represents the local "x" direction, the red arrow the local "y" direction, and the green arrow the local "z" direction. At the outset, the three local directions correspond with the geocentric axes as shown.

(Sharp-eyed readers will notice the relatively extreme eccentricity of this representation of the ellipsoid, especially in the third and fourth images: the Earth looks somewhat "squashed." This is by design: the software that made these images uses the equations worked out below. Exaggerating the eccentricity emphasizes potentially small errors that can creep into ellipsoidal calculations, making them easier to spot and rectify.)

For reasons that will soon be clear, first rotate the local frame 90 degrees counterclockwise around its z-axis:

Frame 2

Next, rotate this frame around the geocentric y-axis to bring it to the desired latitude. (For this illustration we're headed toward Indonesia, in the southern hemisphere.)

Frame 3

(Notice how the red arrow (local y) is pointing due north and the blue arrow (local x) is pointing due east: it was to establish this x=east, y=north correspondence that the first rotation was done.)

Finally, rotate this frame once more around the geocentric z-axis to bring it to the desired longitude.

Frame 4

These three rotations provide the Euler angles in the Z-Y'-Z'' convention. The amounts of rotation were:

First, 90 degrees around the z axis from x towards y.
Second, 90 - 'latitude' degrees around the y axis from z towards x.
Third, 'longitude' degrees around the z axis from x towards y.

(It should now be evident that these transformations were worked out in reverse: starting with a local east-north-up frame at a given point, a rotation around the earth's axis (geocentric z) can take us to the Prime Meridian, and thence a rotation around the geocentric y axis can take us to the North Pole, where we discover a simple 90 degree rotation around the pole will align the local frame with the geocentric one. The key idea is to use rotations around the geocentric axes to line up the frames: although there are many ways to accomplish this, they will all produce the same answer.)

It remains to convert the geometry into arithmetic. The 3 x 3 matrix representing a rotation of t from axis i towards axis j will have values of cos(t) in its (i,i) and (j,j) positions, sin(t) in the (j,i) position, -sin(t) in the (i,j) position, a 1 in the other diagonal position, and zeros elsewhere. Successive rotations are computed by multiplying the matrices (from right to left). Writing f for the longitude and q for the latitude, the resulting rotation matrix is

-sin(f), -cos(f) sin(q), cos(f) cos(q)
 cos(f), -sin(f) sin(q), sin(f) cos(q)
      0,         cos(q),        sin(q)

For the values of f and q in the question, the matrix is

 0.94765709  0.1887300   -0.2575404
-0.3192899   0.56015335  -0.76438353
 0           0.80660351   0.59109287

The first column is a unit-length vector giving the local x (east) direction in terms of the geocentric coordinates (as a check, notice it has no third component: x is entirely east-west). The second column is a unit vector giving the local y (north) direction and the third column is the local outward unit normal, or "up" direction.

As far as the origin of the local system goes, we can work it out from properties of the Earth's ellipsoid. The WGS 84 ellipsoid, for instance, has a major axis of a = 6378137 meters (that's the equatorial radius) and a minor axis of b = 6356752.3142 meters (that's the polar radius). (The ratio of these axes equals 0.996647189335, which is less than 1 by the "flattening," equal to 1 / 298.257223563.) Consequently, the geocentric coordinates of the origin are given by

a cos(q') cos(f) 
a cos(q') sin(f) 
b sin(q')

where q' = Atan(b/a tan(q)) converts the geodetic latitude q to the angle between the origin and the equator. In the example, these coordinates (in meters) are

-1,644,552.1
-4,881,054.5
 3,749,220.2

Finally, to find the geocentric coordinates of a local displacement from this origin, we form the linear combination in terms of the geocentric frame. E.g., a displacement of 250 meters west is -250 in the local x direction, equivalent to a geocentric displacement of

-250 * (0.94765709, -0.3192899, 0) = (-236.9, 79.8, 0)

and a displacement of 250 meters south is -250 in the local y direction, equivalent to a geocentric displacement of

-250 * (0.1887300, 0.56015335, 0.80660351) = (-47.2, -140.0, -201.7).

Applying these two displacements to the geocentric origin (by adding them component-by-component) yields geocentric coordinates of

(-1644836.2, -4881114.7, 3749018.6).

As the question notes, all this can be compactly represented as a four by four transformation matrix: its upper left three by three block is the rotation matrix, the first three entries in its rightmost column are the geocentric coordinates for the given longitude and latitude, and its bottom row is (0,0,0,1). Local coordinates (x',y',z') are converted to geocentric coordinates by left-multiplying the column vector (x',y',z',1)' by this matrix.

Incidentally, the earth's surface is usually not at sea level at land locations like this one. No problem: the surface elevation at the origin represents a displacement in the local z direction. Simply add that multiple of the local z vector to obtain the final value. In the example, the surface elevation at the origin is approximately 2000 meters. The corresponding geocentric displacement therefore equals

2000 * (-0.2575404, -0.76438353, 0.59109287) = (-515.081, -1528.77, 1182.19)

Adding this to the previous result gives

(-1645351.3, -4882643.5, 3750200.7)

Because non-sea-level contributions are considerable, elevation cannot be ignored when computing geocentric coordinates.

[GIS] Converting latitude longitude coordinate to x y coordinate

I assume your data are stored inside a PostGIS database. You can use PostGIS functions to do the transformations directly. No need for external libraries.

    select 
        st_x(st_transform(geom,<SRID>)) as X, 
        st_y(st_transform(geom,<srid>)) as Y 
    from 
        mTable 
    where 
        ...

Best Answer

Related Solutions

Local Coordinate to Geocentric Conversion – How to

[GIS] Converting latitude longitude coordinate to x y coordinate

Related Question