Translating geocoordinates to local x,z coordinates

[domenico]

I am seeking to map geolocations lon,lat onto a 3D scan in local Three.js coordinates x,z. I have two reference points so geolocation and x,z coordinate space can be matched. For a ‘third’ lon, lat at any position I would like to calculate its x,z coordinate with the help of the reference points.

The image below shows two reference points (small spheres) and the hypothetical third point (the large sphere) which requires an x,z coordinate calculation. Nb. the lon,lat input will not be further than 1km from the reference points — so any complex calculation shouldn’t be required to account for the earth’s curvature.

I’m a novice with this kind of mathematics, but I imagine there has to be a straightforward way to translate geocordinates to local coordinates, if I have at least two reference points.

[ccammack]

I’m guessing, but you might need a 3rd known point in the real world to make this work. If you have 3 known points in the world, you could treat one of them as an origin and subtract it from the other 2 points to create the X and Z axes of a world basis. Find those same known points in the image and do the same subtractions to create an image basis. Generate the matrix that transforms image space to world space and transform the unknown image point through that to get the world point of the unknown. The distance between longitude lines changes depending on latitude, so you’ll need to find some geo code to properly calculate the Haversine distance in meters between two (lon, lat) points and use that to create the world basis.

If you get it working, please come back here and report how you did it.

[Attila_Schroeder]

Two points alone are not enough to determine the location of a third point, because how does any third point differ from another? What you want to do looks like triangulation. That means you have two points. You connect these two points with a straight line or on a sphere with a geodesic (a great circle that goes through both reference points and has its origin at the center of the earth). If you now draw straight lines from both reference points, you automatically know their angles to the straight line connecting the two reference points. Then you can completely determine the triangle in the plane using the law of sines. Things look different on a sphere. The shortest path on a sphere between two points is always a great circle, i.e. a circle that goes through both points and has its origin in the center of the sphere. Two great circles on a sphere will always intersect. Is that what you want? Go from two given points in a given direction and then find out where these paths intersect?

At 1km you can work with the sine theorem, but you can’t avoid the transformations between spherical coordinates and Euclidean coordinates. First you have to convert the geo-coordinates into Euclidean coordinates, then use these to find your third point using triangulation and then convert this into spherical coordinates.

If long 3 and lat 3 are given as I assume from your picture and you are looking for the corresponding x, y, z coordinates, that’s easy:

theta: lat (-90 to 90)
phi: long (-180 to 180) or the same (0 to 360)

Remember that this is a Euclidean coordinate system that lies at the center of the earth and not flat on the surface at the respective point. So you will get very large x, y, z values. The question that immediately comes to mind is what is your reference system?

If you choose a coordinate system according to your own ideas, one that seems practical to you in relation to your local environment in order to work with x, z coordinates, you must necessarily know the orientation of this local system to the global one, to which theta and phi refer. If you have two reference points you will most likely view their connecting line as one of the x or z axes.

Then I unfortunately have to say that it doesn’t work without formula magic. To be more precise, tensor calculation.

[manthrax]

ChatGPT says:

function latLongToWebMercator(lat, lon) {
    const R = 6378137; // Earth radius in meters
    const d2r = Math.PI / 180; // Degrees to radians

    let x = R * lon * d2r;
    let y = R * Math.log(Math.tan((Math.PI / 4) + (lat * d2r / 2)));

    return { x, y };
}

For:

4. Web Mercator Projection
This is commonly used in web mapping applications like Google Maps:

[domenico]

@ccammack @Attila_Schroeder @manthrax thank you all for the replies. Based on all this information it’s going to take me a bit of time to process and come back with a solution. In the meantime I should clarify:

All values are known except for x3, z3. I would presume the three geolocations and two x,z values are enough. Even if it’s just a linear transformation because we’re dealing with small distances.

What I would like to be able to try, although I can’t yet think of the maths to do it (even with ChatGPT! I’ve tried and it seems to struggle remembering all the parameters), is to draw a triangle between each point, calculate all the angles including the direction of the triangle (using what we know about the position of the geolocations) and then switch out the lon, lat values for x,z and calculate the unknown x,z. This method would likely use the Haversine Formula to calculate the side lengths with some transformation to chord lengths to take out the curvature. Surely this would yield the result

If I figure this out I’ll post it.

[domenico]

@ccammack @Attila_Schroeder @manthrax Here’s the code I have so far. I feel like I’m almost there — I know what I need but I don’t know the correct math.

ChatGPT has been valuable in writing the formulas (it’s like a translator!) and I’ve validated that the distances outputted for the given geocoordinates are correct. The Haversine to chord conversion is negligible, but since in theory the x,z space is not curved I’ll keep it in.

The logic that’s missing is how to translate the distance from point 1 to 3 into an x,z coordinate and use the bearings I have calculated to choose the right side n.b. if this part isn’t added, the third point could fall on either side of the line created by points 1 & 2. The logic below “below is a work in progress” is incorrect.

// converts degrees into radians
function toRadians(degrees) {
    return degrees * Math.PI / 180;
}

// converts radians into degrees
function toDegrees(radians) {
    return radians * 180 / Math.PI;
}

// calculates the distance between two geocoordinates in meters
function haversineDistance(lat1, lon1, lat2, lon2) {
    const R = 6371000; // Radius of the Earth in meters
    const dLat = toRadians(lat2 - lat1);
    const dLon = toRadians(lon2 - lon1);
    const a = 
        Math.sin(dLat / 2) * Math.sin(dLat / 2) +
        Math.cos(toRadians(lat1)) * Math.cos(toRadians(lat2)) * 
        Math.sin(dLon / 2) * Math.sin(dLon / 2);
    const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
    const d = R * c;
    return d;
}

// converts curved haversine distance to the chord distance i.e. straight line
// the conversion is for practical purposes negligible when used over short distances 
function haversineToChord(haversineDist) {
    const R = 6371000; // Radius of the Earth in meters
    return 2 * R * Math.sin(haversineDist / (2 * R));
}

// calculate the bearing of two geolocations with respect to true north
// the result is the degrees off true north and helps with determining direction
function calculateBearing(lat1, lon1, lat2, lon2) {
    const dLon = toRadians(lon2 - lon1);
    const y = Math.sin(dLon) * Math.cos(toRadians(lat2));
    const x = Math.cos(toRadians(lat1)) * Math.sin(toRadians(lat2)) -
        Math.sin(toRadians(lat1)) * Math.cos(toRadians(lat2)) * Math.cos(dLon);
    let brng = Math.atan2(y, x);
    brng = toDegrees(brng);
    return (brng + 360) % 360;
}

// calculate the x,z value of the 'third' unknown coordinate
// using two control points with known lon,lat and x,z
function thirdPointXZ(lon3, lat3, controlPoints) {
    const { lon1, lat1, x1, z1, lon2, lat2, x2, z2 } = controlPoints;//provide control points as constants to access

    const dist12 = haversineDistance(lat1, lon1, lat2, lon2);// haversine distance between control points
    const chord12 = haversineToChord(dist12);// chord distance between control points

    const dist13 = haversineDistance(lat1, lon1, lat3, lon3);// haversine distance between control point 1 and the 3rd point
    const chord13 = haversineToChord(dist13);//chord distance between control point 1 and the 3rd point

    const bearing12 = calculateBearing(lat1, lon1, lat2, lon2);// true north bearing of control points
    const bearing13 = calculateBearing(lat1, lon1, lat3, lon3);// true north bearing of point 1 to 3

    const bearingDelta = bearing13 - bearing12;

    console.log('1 to 2 dist ' + dist12 + ', 1 to 3 dist ' + dist13)
    console.log('1 to 2 dist chord ' + chord12 + ', 1 to 3 dist chord ' + chord13)
    console.log('1 to 2 bearing ' + bearing12 + ', 1 to 3 bearing ' + bearing13)

    // BELOW IS WORK IN PROGRESS
    // need to apply x,z in position of lon3, lat3
    // compare the 1 to 2 vs 1 to 3 bearings if the difference is greater than 180deg

    const dx = x2 - x1;
    const dz = z2 - z1;

    const distRatio = chord13 / chord12;

    const x3 = x1 + dx * distRatio * Math.cos(toRadians(bearingDelta));
    const z3 = z1 + dz * distRatio * Math.sin(toRadians(bearingDelta));
    
    return { x: x3, z: z3 };

}

// Real data is passed to generate x3,z3
// the control points act as a reference to scale points on a model with real world  geocoordinates
// controlPoints = {
// lon1: 0,// longitude for CP1
// lat1: 0,// latitude for CP1
// x1: 0,// x value in Three.js space for CP1
// z1: 0,// y value in Three.js space for CP1
// lon2: 0,// longitude for CP2
// lat2: 0,// latitude for CP2
// x2: 0,// x value in Three.js space for CP2
// z2:0// y value in Three.js space for CP2
// };
//
// thirdPointXZ(lon3, lat3, controlPoints) 

[PavelBoytchev]

I wrote a function, but it is not tested with real data. With fake data (lines 26-30) it gives the correct result:

https://codepen.io/boytchev/pen/qBLjVKZ?editors=0011

I’m using an algorithm, that does not require circles and intersections. However, it requires that visually 1lon = 1lat (which means you may need to scale input and output data if the location is closer to the poles, or if the view point is not from the top).

Edit: Here is the algorithm: left is the known triangle, right is the triangle with missing vertex C. Imagine you are in point A, you walk m meters to B, then turn right α degrees and walk another n meters to point C.

Now, repeat the same with the other triangle. You start from point A’, you walk M meters to B’, then turn right α degrees and walk another N meters to point C’ – N is unknown, but it is proportional to n (i.e. N/M = n/m):

[domenico]

@PavelBoytchev Pavel this is amazing! When my brain is fresh tomorrow morning (Australian Eastern time) I am going to unpack what you’ve made and see how it works… It’s really interesting to see that you’re using Three.js vectors to run calculations.

[domenico]

@PavelBoytchev Unfortunately I’m not getting consistent results — anything within the reference points area works well, but outside of that I’m not seeing the expected x,z for given lon,lat. Below is sample data, with the expected x,z and the resulting x,z. I also have all the points visualised here Third point x,z calc debug - Google My Maps

Sample Point 2 provides the most curious result. Some variance is OK because a model might not be 100% spatially accurate. However I cannot understand what about sample point 2 would produce that kind of result

// Reference Points
// lat1: -38.1495921086427
// lat2: -38.1492866772213
// lon1: 144.30724288969
// lon2: 144.307338593531
// x1: -2.34
// x2: -3.28
// z1: 12.62
// z2: 11.2

// Sample Point 1
// lon3: 144.3073613926176
// lat3: -38.14901617699974

var P = geoLocal(
	new THREE.Vector3(144.30724288969,0,-38.1495921086427),// lon1, 0, lat1
	new THREE.Vector3(144.307338593531,0,-38.1492866772213),// lon2, 0, lat2
	new THREE.Vector3(144.3073613926176,0,-38.14901617699974),// lon3, 0, lat3
	new THREE.Vector3(-2.34,0,12.62),// x1, 0, z1
	new THREE.Vector3(-3.28,0,11.2)// x2, 0, z2
)

// Expected (approx): -4.21, 10.16
// Result: -3.7957815693929535 9.850952647135689


// Sample Point 2
// lon3: 144.30823108398587
// lat3: -38.14923334727532

var P = geoLocal(
	new THREE.Vector3(144.30724288969,0,-38.1495921086427),// lon1, 0, lat1
	new THREE.Vector3(144.307338593531,0,-38.1492866772213),// lon2, 0, lat2
	new THREE.Vector3(144.30823108398587,0,-38.14923334727532),// lon3, 0, lat3
	new THREE.Vector3(-2.34,0,12.62),// x1, 0, z1
	new THREE.Vector3(-3.28,0,11.2)// x2, 0, z2
)

// Expected (approx): -0.88, 8.24
// Result: -7.920793570380081 12.24465589317603


// Sample Point 3
// lon3: 144.30748132466096
// lat3:  -38.14995687807887

var P = geoLocal(
	new THREE.Vector3(144.30724288969,0,-38.1495921086427),// lon1, 0, lat1
	new THREE.Vector3(144.307338593531,0,-38.1492866772213),// lon2, 0, lat2
	new THREE.Vector3(144.30748132466096,0,-38.14995687807887),// lon3, 0, lat3
	new THREE.Vector3(-2.34,0,12.62),// x1, 0, z1
	new THREE.Vector3(-3.28,0,11.2)// x2, 0, z2
)

// Expected (approx): -0.16, 13.2
// Result: -3.020415227359035 14.836481585990997

[PavelBoytchev]

That difference between expected and calculated results for points 1 and 3 is huge. I would not be happy with such precision. As for point 2, the difference indicates something is very wrong.

The online map is very helpful!

I tried to debug the coordinates that you get, but I was unable to guess the origin and the scale of the local coordinate system that you used. When I introduce my own system (the yellow one in the image below), all coordinates and calculations for the three sample points look OK:

• the red labels are (lon,lat) coordinates taken from Google maps URLs – they differ slightly from the ones that you use; for shorter numbers I kept only 6 digits after the decimal point, I feel such precision is sufficient

• the yellow labels are (x,z) coordinates that use the yellow coordinate system

• I added scaling of input data (lines 9-13) is the source code, because the points are far from the Equator

• the coordinates of the 3 test cases are also included in the source code (it is the same CodePen)

Could you check whether you agree with the results that I get?

If you cannot find why we get different results for correctness, could you send me how is your local coordinate system, so that I could check it too (attached is PPTX file for the image above: TheMap.pptx (251.7 KB) )

[domenico]

Cheers Pavel for all of this. Attached is a screenshot of the model in Three.js. What I’ve done is laid out markers at intervals to give you a sense of the real x,z coordinates. The green markers are the reference points (the bigger one is x2,z2). I’ve superimposed the x,z values. I’ve tried out your new code, and the same issues are happening, just different output positions.

If none of this information helps — I wonder if there’s anything salvageable from my code above as a starting point? I can certainly validate everything up to the point where it says ‘work in progress’ is working as it should. I have accurate distances and bearings, and if the bearing between AB and AC is greater or less than 180deg the side on which x3,z3 is can be determined. I just don’t know how to transform said information from lat, lon to x,z — or more specifically: turning the calculated lengths (and what we know about the x1,z1 & x2,z2 values of the reference points) into the x3,z3 calculation (based on the bearing & bearing evaluation).

[PavelBoytchev]

I noticed, that our local coordinate systems were flipped – i.e. your (x,z) corresponds to my (z,x). When I swap the local coordinates in the input, the results are much closer.

Here is what I get (showing swapped coordinates):

Point       Expected         Result
  #1      10.16 -4.16     10.14 -4.29
  #2       8.24 -0.88      8.21 -0.92
  #3      13.20 -0.16     13.21 -0.20

[Attila_Schroeder]

I wanted to do the math myself because I was curious. What are the x and z units of measurement? If I calculate the distance from point 1 to point 2 using the longitude and latitude I get 35 m.
In your local coordinates, P1 and P2 are just about 1.7 apart. What kind of length scale is that?

[PavelBoytchev]

In my calculations x and z are measured in whatever units you like. It should not matter if you keep sufficiently far from Planck length.

[domenico]

@PavelBoytchev this is an amazing result! Just so I don’t make any error on my end could you please specify exactly how to make this switch. I also wonder whether the same result applies in both v1 and v2 of your code? Comments in the code would be super helpful if you can add them? Especially for anyone looking at this in future.

@Attila_Schroeder x,z are completely arbitrary. The reference points provide meaning.

[PavelBoytchev]

I just exchanged the X and Z local coordinates in the input vectors:
https://codepen.io/boytchev/pen/qBLjVKZ

As of comments, do you mean comments in the body of the geoLocal function?

Edit: the function is now commented.

[Attila_Schroeder]

Since the two reference points automatically create a scale and a fixed relationship to the geo-coordinates, the mathematical clarity is guaranteed. I was just curious because I don’t know the scaling.
@PavelBoytchev Sounds like you have the solution.

[domenico]

Thank you for the solution Pavel.

It’s curious to see too that you have the Y coordinate space empty.

In my context, I haven’t referenced altitude, but it is possible to get altitude information for the reference points and user. So a geo lon, altitude and lat could be marked in x,y,z respectively — unless there’s an extra operation required to process for Y.

[PavelBoytchev]

All the vectors are 3D, because when the signed angle is calculated, the calculation temporarily goes into 3D (lines 37-40). Apart from this calculation, y is not needed and it is 0. I have never thought of using altitude, but my expectation is that a similar calculation should be possible, but it will need two angles, not just one.