1

Sample coordinates in EPSG:4326 (WGS84)

I am trying to find the latitude and longitude of a 3rd point on a map with EPSG:4326 (WGS84) projection. I know the coordinates of two other points and the distance in meters to the 3rd point. Ideally, I am trying to create a function to return said coordinates of the 3rd point.

4
  • 4
    Is the azimuth not missing ? Or are you trying to find all solutions ? Nov 28, 2019 at 8:24
  • Yes I do not have the azimuth, I guess finding all solutions and then figuring out the suitable answer is the only approach?
    – Brad w
    Nov 28, 2019 at 8:33
  • 2
    Then, I think to use a library (like shapely but project your coordinates before) for build a circular ring for each point and find intersection points of the rings. Nov 28, 2019 at 9:05
  • Depending on the software you use, @J. Monticolo's method is the way to go except that you will need to figure out which intersection point is the one you want as there will be 2 points.
    – GforGIS
    Nov 28, 2019 at 9:13

2 Answers 2

2

Note that in your example, it is no solution because the two points are distanced by 30 kilometers.

Here a solution with pyproj for manage projection from EPSG: 4326 WGS84 to EPSG: 3857 Pseudo-Mercator, and shapely :

from pyproj import Proj, transform
from shapely.geometry import Point

crs_in = Proj(init='epsg:4326')
crs_out = Proj(init='epsg:3857')
x1, y1 = 140.1, 35.4  # your first point
distance1 = 200  # in meters
# x2, y2 = 140.2, 35.6  # your second point
x2, y2 = 140.1, 35.402  # a second point with 2 solutions
distance2 = 240  # in meters

x1p, y1p = transform(crs_in, crs_out, x1, y1)
x2p, y2p = transform(crs_in, crs_out, x2, y2)

circle1 = Point(x1p, y1p).buffer(distance1).exterior
circle2 = Point(x2p, y2p).buffer(distance2).exterior

inter12 = circle1.intersection(circle2)

if inter12.is_empty:
    print("No solutions")
elif inter12.type == "Point":
    xsol, ysol = transform(crs_out, crs_in, inter12.x, inter12.y)
    print(f"Solution: (x: {xsol}, y: {ysol})")
elif inter12.type == "MultiPoint":
    i = 0
    for geom in inter12.geoms:
        i += 1
        xsol, ysol = transform(crs_out, crs_in, geom.x, geom.y)
        print(f"Solution {i}: (x: {xsol}, y: {ysol})")

The result is :

Solution 1: (x: 140.09846905765485, y: 35.40076316724516)
Solution 2: (x: 140.10153094234514, y: 35.40076316724516)
2
  • Thanks for the feedback. Does this solution takes the spherical distance into consideration? For example, like hubeny's formula etc.
    – Brad w
    Dec 2, 2019 at 1:38
  • It's just a cartesian distance as long as I projected the coordinates to pseudo-mercator (cartesian, units in meters). And after, I retransform coordinates in latitude / longitude. Dec 2, 2019 at 8:17
0

If you do not have the bearing or azimuth it will not be possible to project the points given a starting point and distance unless if you create programmatically 2 circles with the initial point as the center and use the distance as a radius. The intersection of the two circles would return 2 points of which one is the one of interest. If no azimuth is given then both are the results. The intersection could return nil if one circle is located inside another.

1
  • Will this approach take the spherical distance of the earth into consideration? And can I find out the bearing or azimuth from an image, and use it for the calculations?
    – Brad w
    Dec 2, 2019 at 1:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.