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.

nmtoken
  • 13,355
  • 5
  • 38
  • 87
Brad w
  • 11
  • 1
  • 4
    Is the azimuth not missing ? Or are you trying to find all solutions ? – J. Monticolo Nov 28 '19 at 08: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 '19 at 08: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. – J. Monticolo Nov 28 '19 at 09: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 '19 at 09:13

2 Answers2

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)
J. Monticolo
  • 15,695
  • 1
  • 29
  • 64
  • Thanks for the feedback. Does this solution takes the spherical distance into consideration? For example, like hubeny's formula etc. – Brad w Dec 02 '19 at 01: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. – J. Monticolo Dec 02 '19 at 08: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.

GforGIS
  • 3,126
  • 3
  • 19
  • 38
  • 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 02 '19 at 01:40