1

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.

2
  • 3
    Are you just looking for Haversine? It might help if you expanded your question to explain what constraints (e.g. distance limits, accurancy requirements, particular areas you do / don't care about, efficiency issues, platform issues) that you have to deal with, what you've already tried, and what did/didn't work for you. Just click edit above to update your question with these details.
    – BradHards
    Jun 21, 2013 at 11:57
  • @BradHards i have edited please check Jun 21, 2013 at 12:01

3 Answers 3

2

Background

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?.)

Formulas

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.

Cautionary note

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).

Data

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.)

Worked Example

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.

Limitations

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.

2
  • can you give the formula of calulating angle between equator to a latlng Jun 21, 2013 at 14:07
  • Let f be the geodetic latitude: by definition, that's the angle between "up" and the Equatorial plane. It's what is usually meant by "latitude." Let t be the geocentric latitude, which is the angle subtended by your location and the Equatorial plane. They are connected by atan(t) = btan(f). So if by "angle between equator" and a point you mean t, use this formula; if you mean f, then it's already given to you by the latitude.
    – whuber
    Jun 21, 2013 at 14:13
0

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.

0

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')

Your Answer

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

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