Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
I have a latlng i want to find four points with equal distance from this latlangs in all directions east, west, south an north is there any formula .
I am trying to build polygons around higways so draw polygon around poly line the distance will be less than 1 km.
For these purposes, the Earth is assumed to have the shape of an oblate ellipsoid of revolution around its axis. The major axis, a, is the radius of the Earth at its equator. The minor axis, b, is the radius of the Earth at either pole. The squared eccentricity (measuring the "squishedness") of this ellipse is e^2 = 1 - (b/a)^2.
For simpler but approximate formulas, we may ignore the eccentricity and take the earth to be a perfect sphere. In the formulas the follow, set e^2 = 0 and a = b = (2a+b)/3, which I will write as R. (The justification for this is found at How accurate is approximating the Earth as a sphere?.)
Distances east-west lie along a parallel of latitude, which is a perfect circle. At the geodetic latitude f, the radius of this circle is
r = a * cos(f) / sqrt(1 - e^2 sin(f)^2).
(The spherical approximation is r = R * cos(f).) Therefore, a distance x in the east-west direction corresponds to an angle of
dl = x / r
radians. Convert that to degrees (or grads or whatever you prefer to measure longitude in) and add that value to the longitude (for an east direction) or subtract it (for a west direction).
Laying distances off in the north-south direction requires computing elliptic integrals. However, for the very short distances contemplated in this question (just a few kilometers, or less than 0.001 times a), we may assume the meridian is approximately circular. The best approximating circle ("osculating circle") has radius
s = a * (1 - e^2) / sqrt(1 - e^2 * sin(f)^2)^3.
(The spherical approximation is s = R.)
Therefore, a distance y in the north-south direction corresponds to an angle of
df = y / s
radians. Add that value to the latitude (for a north direction) or subtract it (for a south direction). If the result takes you over a pole, you will need then to subtract it from Pi radians (or 180 degrees or whatever) and add or subtract Pi radians (180 degrees) to the longitude.
The east-west points are actually closer to the original point than you might think. This is because the desired distance is laid off along a circle of latitude rather than a straight line. For instance, a point 10,000 kilometers to the east of a location at 40 degrees latitude is located only 9100 kilometers away. These differences are inconsequential for short distances (less than about 1000 kilometers).
In the World Geodetic System, the major axis has length a = 6378137 meters, minor axis b = 6356752.3142 meters, and inverse flattening f = 1/298.257223563. This is equivalent to a squared eccentricity e^2 = 0.00669438000426. (These values are connected by e^2 = 1 - (b/a)^2 = 2f - f^2.)
Consider the location at (39.9522° N, 75.1642° W). We will find four points exactly one kilometer away in the east, west, north, and south directions. First compute
sin(f) = sin(39.9522°) = 0.642148
sqrt(1 - e^2 sin(f)^2) = sqrt(1 - 0.00669438000426 * 0.642148^2) = 0.998619
Plugging these into the equations gives
r = 4,896,117.456 meters
s = 6,361,763.203 meters
(These values are over-precise for this application, whose coordinates are given only to six significant figures: the extra precision will help when checking a software implementation.) The offsets are
dl = 1000 / r = 0.0002042434662 radian = 0.0117023 degree
df = 1000 / s = 0.0001571891264 radian = 0.00900627 degree
That places the four points at
(39.9522000, 75.1759023), (39.9522000, 75.1524977) (east and west)
(39.9612063, 75.1642000), (39.9431937, 75.1642000) (north and south).
In the spherical approximation, the four points are
(39.9522000, 75.1759316), (39.9522000, 75.1524684) (east and west)
(39.9611601, 75.1642000), (39.9432399, 75.1642000) (north and south).
The approximate locations are 2.5 meters too far in the east-west directions and 1.45 meters too close in the north-south directions. This is a typical level of error (it is approximately the size of the flattening, which is around one part per 300, or 3.3 meters out of 1000 meters). When this much error is acceptable, the simpler spherical formulas are fine; otherwise, you need the ellipsoidal formulas.
The east-west results for the ellipsoid are completely accurate, but the north-south results depend on approximating an elliptical meridian by a circle. They will degrade in accuracy over longer distances, but will still be pretty good out to a few thousand kilometers. For instance--to continue the example--the point the formula places a thousand kilometers north is actually 1000.79 kilometers away, an error of less than one part in 1200.
Here's a solution in Octave or MATLAB. It requires MATLAB package 39108.
octave:1> format long octave:2> [lat, lon] = geodreckon(39.9522, -75.1642, 1000, 90*[0:3]); octave:3> [lat; lon]' ans = 39.9612062665079 -75.1642000000000 39.9521994093929 -75.1524977114628 39.9431937194577 -75.1642000000000 39.9521994093929 -75.1759022885372
This computes the 4 points 1km NESW from (39.9522°N, 75.1642°W). The results are essentially exact. Note that the points follow geodesics from the starting point, so that the EW points end up at a slightly lower latitude.
Here is a solution I used to use
colab for this https://colab.research.google.com/drive/1-i2UDQ8DWRRMRmqGzdCvt6qh2gu9oX7f?usp=sharing
def create_random_point(x0,y0,distance):
"""
Utility method for simulation of the points
"""
r = distance/ 111300
u = np.random.uniform(0,1)
v = np.random.uniform(0,1)
w = r * np.sqrt(u)
t = 2 * np.pi * v
x = w * np.cos(t)
x1 = x / np.cos(y0)
y = w * np.sin(t)
return (x0+x1, y0 +y)
fig = plt.figure()
ax = host_subplot(111, axes_class=AA.Axes)
#ax.set_ylim(76,78)
#ax.set_xlim(13,13.1)
ax.set_autoscale_on(True)
longitude1,latitude1 = -81.81269,39.33015
ax.plot(latitude1,longitude1,'ro')
