"Turn" the method of calculating distances by latitude and longitude coordinates (latitude distance calculation)

Source: Internet
Author: User
Tags asin

Recently, the online search for the "method of calculating distances by latitude and longitude coordinates" reveals that most of the online code is as follows:

#define PI 3.14159265

static double Rc = 6378137; Equatorial radius

static double Rj = 6356725; Polar radius

Class JWD

{

Public

Double M_longitude, M_latitude;

Double M_radlo, M_radla;

Double Ec;

Double Ed;

Public

JWD (double longitude, double latitude)

{

M_longitude = longitude;

M_latitude = Latitude;

M_radlo = longitude * pi/180.;

M_radla = Latitude * pi/180.;

Ec = Rj + (RC-RJ) * (90.-m_latitude)/90.;

Ed = Ec * cos (M_RADLA);

}

~jwd () {};

};

Static JWD getjwdb (Jwd A, double x,double y)

{

Double dx=x;

Double dy=y;

Double BJD = (Dx/a.ed + a.m_radlo) * 180./PI;

Double BWD = (dy/a.ec + a.m_radla) * 180./PI;

JWD B (BJD, BWD);

return B;

}

void Main ()

{

Double referla=30.0;

Double referlo=60.0;

Double dx=500.0;

Double dy=60.0;

JWD A (Referla,referlo), B (0.0,0.0);

B=GETJWDB (A,dx,dy);

cout < < "LA =" < < B.m_latitude < < "lo=" < < B.m_longitude < < Endl;

}

The above code, similar to this one, is easiest to find through search engines. Most are reproduced with each other, never annotated and annotated. Cause a lot of trouble for beginners. Especially:

Ec = Rj + (RC-RJ) * (90.-m_latitude)/90.;

Ed = Ec * cos (M_RADLA);

Ec, ed these 2 parameters, some people also posted on csdn asked "ec,ed what is the meaning of the function". But no one has ever seen an answer.

Links to questions:

http://www.gisforum.net/bbs/TopicOther.asp?t=5&BoardID=33&id=155609

http://bbs.csdn.net/topics/320024634

I have just started to get into this problem, in the Chinese world can only search for these Ctrl-C to CTRL + V things. There is no answer on the world's largest Chinese internet. People who already know it is simple, but beginners must be puzzled. What's more, lbs applications are already very much. There must be a lot of people who have found a better way to answer it. But for me, with the most commonly used keywords, the easiest to find the copy of the answer, but a lot of doubt. In the curious heart, after some understanding, I will give you the results of analysis. Here's the code for C #:

        Public Const Double Ea = 6378137;   Equator Radius Public          const Double Eb = 6356725;   Polar radius         private static void getjwdb (double Glat, double Glon, double distance, double angle, out double BJD, out do Uble BWD)        {            Double dx = distance * + math.sin (angle * math.pi/180.0);            Double dy = distance * Math.Cos (angle * math.pi/180.0);             Double EC = 6356725 + 21412 * (90.0-glat)/90.0;            21412 is the difference between the equatorial radius and the polar radius            Double ec = Eb + (ea-eb) * (90.0-glat)/90.0;             Double ed = EC * Math.Cos (Glat * math.pi/180);            BJD = (dx/ed + glon * math.pi/180.0) * 180.0/MATH.PI;            BWD = (Dy/ec + glat * math.pi/180.0) * 180.0/MATH.PI;        }

Above this function a look is the person who understands Chinese, the name is called Getjwdb, obtains latitude and longitude degree. He obtains another latitude and longitude coordinate according to the longitude, latitude, distance and one angle of the input, and the output parameter is BJD, BWD.

1) Angle * math.pi/180.0 effect is to convert the angle to radians, longitude coordinates is the angle value, the calculation needs to be converted to radians. All calculations here are in radians.

2) The function is in the direction of the northern (due north), which is the compass, to 0 degrees, and the clockwise direction increases. For example, if the distance distance is D, the DX is the length of the x-axis direction, that is, the length of the longitude longitude direction, and dy is the length of the y-axis direction, that is, the length of the latitude latitude direction.

DX and dy can also be calculated in the direction of due east (0 degrees). That: dx=distance* Cos (θ), dy=distance*sin (θ). The difference is that cos is interchangeable with sin. (not marked on the chart)

Figure 1

3) The EA represents the equatorial radius, EB represents the polar radius, the earth is an approximate sphere, and the EA has a slight gap with the EB. The EC's role is to correct the length of the ball radius due to the changing latitude. If ec=eb+ (EA-EB) * (90-0)/90=ea at glat=0, i.e. at the equator, the EC is just the equator radius Ea; If at the Pole glat=90,ec=eb+ (ea-eb) * (90-90)/90=eb, then the EC It's just the polar radius EB.

4) Ed is the radius of the latitude circle where the Glat is located, such as:

Figure 2

Cross-section of the sphere, at this time the area of the largest section, the circle called the ball of the large circle (great Cycle), along the Meridian section, the resulting is a great circle (great Cycle). The sphere is not through the sphere. The Circle of the cross section. It's called a small circle. The circle in which the latitude circle is located is a small circle. Earth radius r, average r=6371.0km, ed if using the earth radius r, that is Ed=r*cos (θ), you can see

Based on the 2 latitude and longitude points, calculate the distance between the 2 latitude and longitude points (distance through latitude) "The radius of the circle (AO ') where A and D are located in the latitude" mentioned here. Put it in the function above, because it constantly corrects the Earth radius ec, that is ed = EC * Math.Cos (Glat * math.pi/180);

5) According to the latitude and longitude of the Earth, such as:

The length of the longitude line between the latitudes of each equal is certain. Section A, the length of segment B is the same.

Between the longitude of each dividing line, the length of the latitude segment is reduced from the equator to the 2 pole. C end, the length of the D segment is not the same.

Figure 3

See, that dx/ed is equivalent to the ratio of DX length to total length in the latitude of Glat, calculated as a longitude span. If the longitude span plus the starting given longitude is the final longitude.

The same dy/r is the ratio of the length of dy in the longitude of Glon to the average radius r of the earth, which is calculated to be a latitude span. If this latitude span plus the starting given latitude is the final latitude. This uses R, which takes the average radius of the Earth.

The DY/EC is replacing the average radius r with an ever-revised EC radius.

(Dy/ec + Glat * math.pi/180.0) is the final latitude of the starting latitude plus the distance distance, and the result needs to be converted to an angle. The conversion angle method is: radians * 180.0/math.pi.

BWD = (Dy/ec + glat * math.pi/180.0) * 180.0/MATH.PI;

(dx/ed + Glon * math.pi/180.0) is the final longitude of the starting longitude plus the distance distance, and the result needs to be converted to an angle.

BJD = (dx/ed + glon * math.pi/180.0) * 180.0/MATH.PI;

This is based on a latitude and longitude coordinates, distance and then seek another latitude and longitude coordinate function, mainly is to determine a minimum outsourcing rectangle (Minimum bounding rectangle, referred to as MBR). For example, I'm looking for all the merchant information, attraction information, etc. within a 5-kilometer radius of a coordinate point (Lat,lon). This MBR is the largest range, and this rectangle is a minimum rectangle that contains all of these valid information in the range of 5 km. Using the formula, the MBR can be obtained by finding 270 coordinate points in the direction of 0 degrees, 90 degrees, 180 degrees and four degrees in four directions.

        Public Const Double Ea = 6378137;     Equator radius Public Const Double EB = 6356725; Polar radius private static void Getlatlon (double LAT, double LON, double distance, double angle, out double newLon,            Out double Newlat) {Double dx = distance * + * Math.sin (angle * math.pi/180.0);            Double dy = distance * Math.Cos (angle * math.pi/180.0);            Double EC = Eb + (ea-eb) * (90.0-lat)/90.0;            Double ed = EC * Math.Cos (LAT * math.pi/180);            NewLon = (dx/ed + LON * math.pi/180.0) * 180.0/MATH.PI;         Newlat = (Dy/ec + LAT * math.pi/180.0) * 180.0/MATH.PI; public static void Getrectrange (double centorlatitude, double centorlogitude, double distance, out double Maxlat            Itude, out double minlatitude, out double maxlongitude, out double minlongitude) {Double temp = 0.0; Getlatlon (Centorlatitude, centorlogitude, distance, 0, out temp, out MAxlatitude);            Getlatlon (Centorlatitude, centorlogitude, distance, N, out temp, out minlatitude);            Getlatlon (Centorlatitude, centorlogitude, distance, N, out Maxlongitude, out temp);        Getlatlon (Centorlatitude, centorlogitude, distance, +, out minlongitude, out temp); }

Here you get the Maxlatitude, Minlatitude, Maxlongitude, minlongitude that make up the four vertices of the rectangle.

Another method on the web,

Http://www.movable-type.co.uk/scripts/latlong.html

The section "Destination Point given distance and bearing from start point" is described here. I paste the code directly:

Here Getrectrange This function is also to due north as the starting angle, clockwise direction, to obtain Maxlatitude, Minlatitude, Maxlongitude, minlongitude such an MBR. The error of 2 methods is very small. I think it's all a formula that can be used.

         <summary>//Whereφis Latitude,λis Longitude,θis The bearing (clockwise from north), Δis the angular distance d/r;        D being the distance travelled, R the earth ' s radius//bearing azimuth 0,90,180,270//</summary>        private static void Getnewlatlon (double lat, double lon, double D, double bearing, out double lat2, out double lon2)            {lat2 = 0.0;            Lon2 = 0.0;            Double R = 6378.137;            varφ1 = Convertdegreestoradians (LAT);            varλ1 = Convertdegreestoradians (lon);             varθ= Convertdegreestoradians (bearing); varφ2 = Math.asin (Math.sin (φ1) * Math.Cos (D/R) + math.cos (φ1) * Math.sin (D/R) * Math.Cos (            θ));  Varλ2 =λ1 + math.atan2 (math.sin (θ) * Math.sin (D/R) * Math.Cos (φ1), Math.Cos (D/R)            -Math.sin (φ1) * Math.sin (φ2)); λ2 = (λ2 + 3 * Math.PI)% (2 * math.pi)- Math.PI;            Normalise to-180..+180°lat2 = Convertradianstodegrees (φ2);         Lon2 = Convertradianstodegrees (λ2); }

If there is an application, the table holds 1 million of the data, the data contains a lat, Lon latitude and longitude information. You can first get an MBR based on the latitude and longitude of the input, and then through a similar

SELECT Id

From idinfotable

WHERE latitude >= Minlat and Latitude < Maxlat

and longitude >= Minlon and longitude < Maxlon

The way to filter out most of the data, and then through the "2 latitude and longitude points, calculate the distance between the 2 latitude and longitude points (distance by latitude)" The method mentioned here is fine filtering.

Full code:

Using system;using system.collections.generic;using system.linq;using system.text;        Namespace getjwd{public class Getmbr {public double maxlatitude;        public double minlatitude;        public double maxlongitude;         public double minlongitude;        public double MaxLatitude2;        public double MinLatitude2;        public double MaxLongitude2;        public double MinLongitude2; Public getmbr (double centorlatitude, double centorlogitude, double distance) {Getrectrange (Centorlatitu            De, centorlogitude, distance, out-maxlatitude, out-minlatitude, out-maxlongitude, out-minlongitude); GetRectRange2 (Centorlatitude, centorlogitude, distance, out MaxLatitude2, out MinLatitude2, out MaxLongitude2, out        MINLONGITUDE2);     } public Const Double Ea = 6378137;     Equator radius Public Const Double EB = 6356725; Polar radius private static void Getlatlon (double LAT, double LON, double distance, double angle, out Double NewLon, out double newlat) {Double dx = distance * + * Math.sin (angle * math.pi/180.0);            Double dy = distance * Math.Cos (angle * math.pi/180.0);            Double EC = Eb + (ea-eb) * (90.0-lat)/90.0;            Double ed = EC * Math.Cos (LAT * math.pi/180);            NewLon = (dx/ed + LON * math.pi/180.0) * 180.0/MATH.PI;         Newlat = (Dy/ec + LAT * math.pi/180.0) * 180.0/MATH.PI;                                      public static void Getrectrange (double centorlatitude, double centorlogitude, double distance,        Out double-maxlatitude, out double-minlatitude, out double-maxlongitude, out double minlongitude)            {Double temp = 0.0;            Getlatlon (Centorlatitude, centorlogitude, distance, 0, out temp, out maxlatitude);            Getlatlon (Centorlatitude, centorlogitude, distance, N, out temp, out minlatitude); Getlatlon (Centorlatitude, centorlogitude, distance, Maxlongitude, out, temp);        Getlatlon (Centorlatitude, centorlogitude, distance, +, out minlongitude, out temp);                                      public static void GetRectRange2 (double centorlatitude, double centorlogitude, double distance,        Out double-maxlatitude, out double-minlatitude, out double-maxlongitude, out double minlongitude)            {Double temp = 0.0;            Getnewlatlon (Centorlatitude, centorlogitude, distance, 0, out maxlatitude, out temp);            Getnewlatlon (Centorlatitude, centorlogitude, distance, N, out Minlatitude, out temp);            Getnewlatlon (Centorlatitude, centorlogitude, distance, N, out temp, out maxlongitude);        Getnewlatlon (Centorlatitude, centorlogitude, distance, he, + temp, out minlongitude);         }///<summary>//Whereφis Latitude,λis Longitude,θis The bearing (clockwise from north), Δis the angular distance d/r; D being the Distance travelled, R the earth ' s radius///bearing azimuth 0,90,180,270//</summary> Private St            atic void Getnewlatlon (double lat, double lon, double D, double bearing, out double lat2, out double lon2) {            LAT2 = 0.0;            Lon2 = 0.0;            Double R = 6378.137;            varφ1 = Convertdegreestoradians (LAT);            varλ1 = Convertdegreestoradians (lon);             varθ= Convertdegreestoradians (bearing); varφ2 = Math.asin (Math.sin (φ1) * Math.Cos (D/R) + math.cos (φ1) * Math.sin (D/R) * Math.Cos (            θ));  Varλ2 =λ1 + math.atan2 (math.sin (θ) * Math.sin (D/R) * Math.Cos (φ1), Math.Cos (D/R)            -Math.sin (φ1) * Math.sin (φ2)); λ2 = (λ2 + 3 * Math.PI)% (2 * math.pi)-Math.PI;            Normalise to-180..+180°lat2 = Convertradianstodegrees (φ2);         Lon2 = Convertradianstodegrees (λ2); } public static Double CoNvertdegreestoradians (double degrees) {return degrees * math.pi/180;         } public static double Convertradianstodegrees (double radian) {return radian * 180.0/MATH.PI; }} class Test {static void Main (string[] args) {Double latorg = 22.54587746             , lonorg = 114.12873077;             var Gpsdis = new Getmbr (latorg, lonorg, 5); Console.WriteLine ("Maxlat:{0}, Minlat:{1}, maxlon:{2}, Minlon:{3}", Gpsdis. Maxlatitude, Gpsdis. Minlatitude, Gpsdis. Maxlongitude, Gpsdis.            Minlongitude); Console.WriteLine ("Maxlat:{0}, Minlat:{1}, maxlon:{2}, Minlon:{3}", Gpsdis. MaxLatitude2, Gpsdis. MinLatitude2, Gpsdis. MaxLongitude2, Gpsdis.        MINLONGITUDE2); }    }}

This article transferred from: http://www.cnblogs.com/softfair/p/lat_lon_distance_bearing_new_lat_lon.html

"Turn" the method of calculating distances by latitude and longitude coordinates (latitude distance calculation)

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.