使用
對上面的公式解釋如下：
- Lung1 Lat1表示A點經緯度， Lung2 Lat2表示B點經緯度；
- a=Lat1 – Lat2 為兩點緯度之差 b=Lung1 -Lung2 為兩點經度之差；
- 6378.137為地球半徑，單位為千米；
- 計算出來的結果單位為千米，若將半徑改為米為單位則計算的結果單位為米。
- 計算精度與谷歌地圖的距離精度差不多，相差范圍在0.2米以下。
參數說明：
lng：經度
lat：緯度
地球半徑：6378.137（千米）
一般地圖上顯示的坐標順序為，緯度在前（范圍-90 ~ 90），經度在后（范圍-180 ~ 180）
使用（計算直線距離）：
這種計算方式一般都是直線距離
sql語句：
SELECT *, 6378.138 * 2 * ASIN( SQRT( POW( SIN( ( '.$lat.' * PI() / 180 - lat * PI() / 180 ) / 2 ), 2 ) + COS('.$lat.' * PI() / 180) * COS(lat * PI() / 180) * POW( SIN( ( '.$lng.' * PI() / 180 - lng * PI() / 180 ) / 2 ), 2 ) ) ) *1000 AS distance FROM distance ORDER BY distance ASC
JS：
/** * 轉換弧度 * @param d * @returns {number} */ function getRad(d){ var PI = Math.PI; return d*PI/180.0; } /** * 根據經緯度計算兩點間距離 * @param lng1 * @param lat1 * @param lng2 * @param lat2 * @returns {number|*} * @constructor */ function CoolWPDistance(lng1,lat1,lng2,lat2){ var f = getRad((lat1 + lat2)/2); var g = getRad((lat1 - lat2)/2); var l = getRad((lng1 - lng2)/2); var sg = Math.sin(g); var sl = Math.sin(l); var sf = Math.sin(f); var s,c,w,r,d,h1,h2; var a = 6378137.0;//The Radius of eath in meter. var fl = 1/298.257; sg = sg*sg; sl = sl*sl; sf = sf*sf; s = sg*(1-sl) + (1-sf)*sl; c = (1-sg)*(1-sl) + sf*sl; w = Math.atan(Math.sqrt(s/c)); r = Math.sqrt(s*c)/w; d = 2*w*a; h1 = (3*r -1)/2/c; h2 = (3*r +1)/2/s; s = d*(1 + fl*(h1*sf*(1-sg) - h2*(1-sf)*sg)); if(s >= 1000 && s <= 99000){ var kilometer = s/1000; s = kilometer.toFixed(1) + 'km'; }else if(s > 99000){ s = '>99km'; }else{ s = Math.round(s) + 'm'; } // s = s/1000; // s = s.toFixed(2);//指定小數點后的位數。 return s; }
Java
第一步：添加Maven依賴
<dependency> <groupId>org.gavaghan</groupId> <artifactId>geodesy</artifactId> <version>1.1.3</version> </dependency>
第二步：代碼實現
public class DistanceUtil { public static void main(String[] args) { System.out.println("經緯度距離計算結果：" + getDistance(109.371319, 22.155406, 108.009758, 21.679011) + "米"); } public static double getDistance(double longitudeFrom, double latitudeFrom, double longitudeTo, double latitudeTo) { GlobalCoordinates source = new GlobalCoordinates(latitudeFrom, longitudeFrom); GlobalCoordinates target = new GlobalCoordinates(latitudeTo, longitudeTo); return new GeodeticCalculator().calculateGeodeticCurve(Ellipsoid.Sphere, source, target).getEllipsoidalDistance(); } }
GPS，經度，緯度，距離在線計算器：
https://www.osgeo.cn/app/s1884
http://www.box3.cn/page/lbs.html
算法實現
問題描述
給定一個景點的經緯度，給定距離，給定形狀，判斷其它點是否在某個區域內：
圓形方案
使用通用的地球上兩點距離函數，圓形只需要判斷距離，正方形需要計算兩次距離（指定經度與景點的經度一樣，計算是否在范圍內，然后指定緯度與景點的緯度一樣，計算是否在范圍內，如果都在范圍內，則代表該點在景區范圍內），其它形狀基於這個基礎類推
距離函數
Scala實現方式：
/** * Created on 2020/05/15. * * 求地球上兩點間的距離 返回的是 double 格式的 km * 第一點經緯度為（lat1,lng1），第二點經緯度為（lat2,lng2），地球平均半徑R=6378.137 * 按照0度經線的基准，東經取經度的正值(Longitude)，西經取經度負值(-Longitude)，北緯取90-緯度值(90- Latitude)，南緯取90+緯度值(90+Latitude)， * 則經過上述處理過后的兩點被計為(MLon1, MLat1)和(MLon2, MLat2)。那么根據三角推導，可以得到計算兩點距離的如下公式： * C = sin(MLat1)*sin(MLat2)*cos(MLon1-MLon2) + cos(MLat1)*cos(MLat2)，Distance = R*Arccos(C)*Pi/180 * 如果僅對經度作正負的處理，而不對緯度作90-Latitude(假設都是北半球，南半球只有澳洲具有應用意義)的處理，那么公式將是： * C = sin(LatA)*sin(LatB) + cos(LatA)*cos(LatB)*cos(MLonA-MLonB)，Distance = R*Arccos(C)*Pi/180 * 三角函數的輸入和輸出都采用弧度值，那么公式還可以寫作: * C = sin(Lat1*Pi/180)*sin(Lat2*Pi/180) + cos(Lat1*Pi/180)*cos(Lat2*Pi/180)*cos((MLon1-MLon2)*Pi/180)，Distance = R*Arccos(C)*Pi/180 * rad()函數求弧度，Distance（）函數求距離 * * 結果驗證工具地址 http://www.storyday.com/wp-content/uploads/2008/09/latlung_dis.html */ def distance(lat1: Double, lng1: Double, lat2: Double, lng2: Double): Double = { val EARTH_RADIUS = 6378.137 val radLat1 = rad(lat1) val radLat2 = rad(lat2) val a = rad(lat1) - rad(lat2) val b = rad(lng1) - rad(lng2) val distance = EARTH_RADIUS * 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin((b) / 2), 2))) // printLog.debug("lat1: " + lat1 + " lng1: " + lng1 + " lat2: " + lat2 + " lng2: " + lng2) printLog.debug("distance:" + distance) distance }
不規則圖形方案
scala實現函數：使用現成的算法PNPoly即可實現
/** * 多點位置判斷算法 判斷一個坐標點是否在不規則多邊形內部 * * * 在 GIS（地理信息管理系統）中，判斷一個坐標是否在多邊形內部是個經常要遇到的問題。乍聽起來還挺復雜。 * 根據 W. Randolph Franklin 提出的 PNPoly (http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html) 算法，只需區區幾行代碼就解決了這個問題: * * * 針對每一個點，算法遍歷多邊形相鄰的每兩個頂點（即一條邊），假如待判斷點滿足以下兩個條件即改變點是否在多邊形內的狀態標識c： * 待判斷點的Y坐標在點i和點j的Y坐標范圍之內 * 待判斷點的X坐標在點i和點j連線之下 * 遍歷所有的邊之后假如以上兩個條件同時滿足奇數次則該帶判斷點位於多邊形之內，否則位於多邊形之外。 * 算法復雜度為O(n)，其中n為多邊形的頂點個數。 * * @param vertexes * @param testPoint * @return */ def pNPoly(vertexes: Array[LocationPoint], testPoint: LocationPoint): Boolean = { var flag = false var flag0 = false var flag1 = false var flag2 = false var flag3 = false var j = vertexes.length - 1 val loop = new Breaks loop.breakable { for (i <- 0 until vertexes.length) { if (i != 0) { j = i - 1 } if ((vertexes(i).lat == testPoint.lat) && (vertexes(i).long == testPoint.long)) { flag0 = true } if (((vertexes(i).lat - testPoint.lat) * (vertexes(i).long - testPoint.long) * (vertexes(j).lat - testPoint.lat) * (vertexes(j).long - testPoint.long) == 0) && (((vertexes(i).lat > testPoint.lat) != (vertexes(j).lat > testPoint.lat)) && ((vertexes(i).long > testPoint.long) != (vertexes(j).long > testPoint.long)))) { flag1 = true } if (((vertexes(i).lat - testPoint.lat) * (vertexes(i).long - testPoint.long) * (vertexes(j).lat - testPoint.lat) * (vertexes(j).long - testPoint.long) != 0) && ((vertexes(i).lat - testPoint.lat) / (vertexes(i).long - testPoint.long) == (vertexes(j).lat - testPoint.lat) / (vertexes(j).long - testPoint.long))) { flag2 = true } if (((vertexes(i).lat > testPoint.lat) != (vertexes(j).lat > testPoint.lat)) && (testPoint.long < (vertexes(j).long - vertexes(i).long) * (testPoint.lat - vertexes(j).lat) / (vertexes(j).lat - vertexes(i).lat) + vertexes(i).long)) { flag3 = true } flag = flag0 || flag1 || flag2 || flag3 if (flag) { loop.break } } } flag }
在運用該算法之前，有個優化方案，就是先進行最小外接矩形范圍判斷，如果不在最小外接矩形中，則直接跳過
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) { // 點在多邊形之外的實現 }
首先判斷點是否在多邊形的最小外接矩形之內，該步驟不是必須的，但是可以有效避免不必要的計算。
判斷點是否在多邊形內算法實現
PNPoly
PNPoly是由W. Randolph Franklin提出的一種在二維空間較為高效地判斷點是否在多邊形內的算法，算法具有很好的通用性，對凸多邊形、凹多邊形以及含有空洞的多邊形均適用。
最小外接矩形范圍判斷
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) { // 點在多邊形之外 }
首先判斷點是否在多邊形的最小外接矩形之內，該步驟不是必須的，但是可以有效避免不必要的計算。
核心算法
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert-1; i < nvert; j = i++) { if ( ((verty[i]>testy) != (verty[j]>testy)) && (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) c = !c; } return c; }
針對每一個點，算法遍歷多邊形相鄰的每兩個頂點（即一條邊），假如待判斷點滿足以下兩個條件即改變點是否在多邊形內的狀態標識c：
- 待判斷點的Y坐標在點i和點j的Y坐標范圍之內
- 待判斷點的X坐標在點i和點j連線之下
遍歷所有的邊之后假如以上兩個條件同時滿足奇數次則該帶判斷點位於多邊形之內，否則位於多邊形之外。
算法復雜度為O(n)，其中n為多邊形的頂點個數。
算法主頁： PNPoly
相關出版物：Haines, Eric, “Point in Polygon Strategies,” Graphics Gems IV, ed. Paul Heckbert, Academic Press, p. 24-46, 1994.
基於掃描填充多邊形的方法
若針對同一多邊形有大量的點需要判斷是否在該多邊形的話上述方法計算量將會很大，可以借助計算機圖形學中的多邊形填充方法，采用影像輔助判斷點是否在多邊形內，基於voronoi圖的影像鑲嵌（正射影像制作）是該方法典型的運用。
首先計算多邊形的最小外接矩形，采用openCV繪制一張與外接矩形大小的黑色圖像（圖像分辨率可自定義，為適應之后的填充算法創建的影像的大小在掃描方向上需多加1個像素），在圖像上繪制多邊形的邊。
image = cv::Mat::zeros(height,width,CV_8UC1);
line(image,pointi,pointj,255);
繪制多邊形邊
對上面的影像進行填充：
for (int i = 0; i < height; i++) { int trans_count(0); bool is_in_pologon(false); int trans_index(0); for (int j = 0; j < width-1; j++) { //cv::Vec3b a = image.at<cv::uchar>(i,j); if (image.at<cv::uchar>(i,j)==white_pixl&&image.at<cv::uchar>(i,j+1)!=white_pixl) { is_in_pologon = !is_in_pologon; trans_index = j; trans_count++; } if (is_in_pologon) { image.at<cv::uchar>(i,j) = 255; } } if (is_in_pologon) { image.at<cv::uchar>(i,width-1); } if (trans_count == 1) { for (int j = width-1; j > trans_index; j--) { image.at<cv::uchar>(i,j) = 0; } } }
掃描填充影像將待判斷點的坐標換算到影像上若得到相應的像素值為255則該點在這個多邊形之內，若得相應像素值為0則該點在多邊形外。