地圖幾何計算

百度地圖 JavaScript 開源庫 

類 BMapLib.GeoUtils

源文件: GeoUtils.js.

http://api.map.baidu.com/library/GeoUtils/1.2/docs/symbols/BMapLib.GeoUtils.html

http://api.map.baidu.com/library/GeoUtils/1.2/docs/symbols/src/BMapLib_GeoUtils.js.html

 示例 http://api.map.baidu.com/library/GeoUtils/1.2/examples/simple.html

 /**
  * @fileoverview GeoUtils類提供若干幾何算法，用來幫助用戶判斷點與矩形、
  * 圓形、多邊形線、多邊形面的關係,並提供計算折線長度和多邊形的面積的公式。 
  * 主入口類是<a href="symbols/BMapLib.GeoUtils.html">GeoUtils</a>，
  * 基於Baidu Map API 1.2。
  *
  * @author Baidu Map Api Group 
  * @version 1.2
  */
 
 /** 
  * @namespace BMap的全部library類均放在BMapLib命名空間下
  */
 var BMapLib = window.BMapLib = BMapLib || {};
 (function() { 
     
     /**
      * 地球半徑
      */
     var EARTHRADIUS = 6370996.81; 
 
     /** 
      * @exports GeoUtils as BMapLib.GeoUtils 
      */
     var GeoUtils =
     /**
      * GeoUtils類，靜態類，勿需實例化便可使用
      * @class GeoUtils類的<b>入口</b>。
      * 該類提供的都是靜態方法，勿需實例化便可使用。     
      */
     BMapLib.GeoUtils = function(){
         
     }
     
     /**
      * 判斷點是否在矩形內
      * @param {Point} point 點對象
      * @param {Bounds} bounds 矩形邊界對象
      * @returns {Boolean} 點在矩形內返回true,不然返回false
      */
     GeoUtils.isPointInRect = function(point, bounds){
         //檢查類型是否正確
         if (!(point instanceof BMap.Point) || 
             !(bounds instanceof BMap.Bounds)) {
             return false;
         }
         var sw = bounds.getSouthWest(); //西南腳點
         var ne = bounds.getNorthEast(); //東北腳點
         return (point.lng >= sw.lng && point.lng <= ne.lng && point.lat >= sw.lat && point.lat <= ne.lat);
     }
     
     /**
      * 判斷點是否在圓形內
      * @param {Point} point 點對象
      * @param {Circle} circle 圓形對象
      * @returns {Boolean} 點在圓形內返回true,不然返回false
      */
     GeoUtils.isPointInCircle = function(point, circle){
         //檢查類型是否正確
         if (!(point instanceof BMap.Point) || 
             !(circle instanceof BMap.Circle)) {
             return false;
         }
 
         //point與圓心距離小於圓形半徑，則點在圓內，不然在圓外
         var c = circle.getCenter();
         var r = circle.getRadius();
 
         var dis = GeoUtils.getDistance(point, c);
         if(dis <= r){
             return true;
         } else {
             return false;
         }
     }
     
     /**
      * 判斷點是否在折線上
      * @param {Point} point 點對象
      * @param {Polyline} polyline 折線對象
      * @returns {Boolean} 點在折線上返回true,不然返回false
      */
     GeoUtils.isPointOnPolyline = function(point, polyline){
         //檢查類型
         if(!(point instanceof BMap.Point) ||
             !(polyline instanceof BMap.Polyline)){
             return false;
         }
 
         //首先判斷點是否在線的外包矩形內，若是在，則進一步判斷，不然返回false
         var lineBounds = polyline.getBounds();
         if(!this.isPointInRect(point, lineBounds)){
             return false;
         }
 
         //判斷點是否在線段上，設點爲Q，線段爲P1P2 ，
         //判斷點Q在該線段上的依據是：( Q - P1 ) × ( P2 - P1 ) = 0，且 Q 在以 P1，P2爲對角頂點的矩形內
         var pts = polyline.getPath();
         for(var i = 0; i < pts.length - 1; i++){
             var curPt = pts[i];
             var nextPt = pts[i + 1];
             //首先判斷point是否在curPt和nextPt之間，即：此判斷該點是否在該線段的外包矩形內
             if (point.lng >= Math.min(curPt.lng, nextPt.lng) && point.lng <= Math.max(curPt.lng, nextPt.lng) &&
                 point.lat >= Math.min(curPt.lat, nextPt.lat) && point.lat <= Math.max(curPt.lat, nextPt.lat)){
                 //判斷點是否在直線上公式
                 var precision = (curPt.lng - point.lng) * (nextPt.lat - point.lat) - 
                     (nextPt.lng - point.lng) * (curPt.lat - point.lat);                
                 if(precision < 2e-10 && precision > -2e-10){//實質判斷是否接近0
                     return true;
                 }                
             }
         }
         
         return false;
     }
     
     /**
      * 判斷點是否多邊形內
      * @param {Point} point 點對象
      * @param {Polyline} polygon 多邊形對象
      * @returns {Boolean} 點在多邊形內返回true,不然返回false
      */
     GeoUtils.isPointInPolygon = function(point, polygon){
         //檢查類型
         if(!(point instanceof BMap.Point) ||
             !(polygon instanceof BMap.Polygon)){
             return false;
         }
 
         //首先判斷點是否在多邊形的外包矩形內，若是在，則進一步判斷，不然返回false
         var polygonBounds = polygon.getBounds();
         if(!this.isPointInRect(point, polygonBounds)){
             return false;
         }
 
         var pts = polygon.getPath();//獲取多邊形點
         
         //下述代碼來源：http://paulbourke.net/geometry/insidepoly/，進行了部分修改
         //基本思想是利用射線法，計算射線與多邊形各邊的交點，若是是偶數，則點在多邊形外，不然
         //在多邊形內。還會考慮一些特殊狀況，如點在多邊形頂點上，點在多邊形邊上等特殊狀況。
         
         var N = pts.length;
         var boundOrVertex = true; //若是點位於多邊形的頂點或邊上，也算作點在多邊形內，直接返回true
         var intersectCount = 0;//cross points count of x 
         var precision = 2e-10; //浮點類型計算時候與0比較時候的容差
         var p1, p2;//neighbour bound vertices
         var p = point; //測試點
         
         p1 = pts[0];//left vertex        
         for(var i = 1; i <= N; ++i){//check all rays            
             if(p.equals(p1)){
                 return boundOrVertex;//p is an vertex
             }
             
             p2 = pts[i % N];//right vertex            
             if(p.lat < Math.min(p1.lat, p2.lat) || p.lat > Math.max(p1.lat, p2.lat)){//ray is outside of our interests                
                 p1 = p2; 
                 continue;//next ray left point
             }
             
             if(p.lat > Math.min(p1.lat, p2.lat) && p.lat < Math.max(p1.lat, p2.lat)){//ray is crossing over by the algorithm (common part of)
                 if(p.lng <= Math.max(p1.lng, p2.lng)){//x is before of ray                    
                     if(p1.lat == p2.lat && p.lng >= Math.min(p1.lng, p2.lng)){//overlies on a horizontal ray
                         return boundOrVertex;
                     }
                     
                     if(p1.lng == p2.lng){//ray is vertical                        
                         if(p1.lng == p.lng){//overlies on a vertical ray
                             return boundOrVertex;
                         }else{//before ray
                             ++intersectCount;
                         } 
                     }else{//cross point on the left side                        
                         var xinters = (p.lat - p1.lat) * (p2.lng - p1.lng) / (p2.lat - p1.lat) + p1.lng;//cross point of lng                        
                         if(Math.abs(p.lng - xinters) < precision){//overlies on a ray
                             return boundOrVertex;
                         }
                         
                         if(p.lng < xinters){//before ray
                             ++intersectCount;
                         } 
                     }
                 }
             }else{//special case when ray is crossing through the vertex                
                 if(p.lat == p2.lat && p.lng <= p2.lng){//p crossing over p2                    
                     var p3 = pts[(i+1) % N]; //next vertex                    
                     if(p.lat >= Math.min(p1.lat, p3.lat) && p.lat <= Math.max(p1.lat, p3.lat)){//p.lat lies between p1.lat & p3.lat
                         ++intersectCount;
                     }else{
                         intersectCount += 2;
                     }
                 }
             }            
             p1 = p2;//next ray left point
         }
         
         if(intersectCount % 2 == 0){//偶數在多邊形外
             return false;
         } else { //奇數在多邊形內
             return true;
         }            
     }
 
     /**
      * 將度轉化爲弧度
      * @param {degree} Number 度     
      * @returns {Number} 弧度
      */
     GeoUtils.degreeToRad =  function(degree){
         return Math.PI * degree/180;    
     }
     
     /**
      * 將弧度轉化爲度
      * @param {radian} Number 弧度     
      * @returns {Number} 度
      */
     GeoUtils.radToDegree = function(rad){
         return (180 * rad) / Math.PI;       
     }
     
     /**
      * 將v值限定在a,b之間，緯度使用
      */
     function _getRange(v, a, b){
         if(a != null){
           v = Math.max(v, a);
         }
         if(b != null){
           v = Math.min(v, b);
         }
         return v;
     }
     
     /**
      * 將v值限定在a,b之間，經度使用
      */
     function _getLoop(v, a, b){
         while( v > b){
           v -= b - a
         }
         while(v < a){
           v += b - a
         }
         return v;
     }
 
     /**
      * 計算兩點之間的距離,兩點座標必須爲經緯度
      * @param {point1} Point 點對象
      * @param {point2} Point 點對象
      * @returns {Number} 兩點之間距離，單位爲米
      */
     GeoUtils.getDistance = function(point1, point2){
         //判斷類型
         if(!(point1 instanceof BMap.Point) ||
             !(point2 instanceof BMap.Point)){
             return 0;
         }
 
         point1.lng = _getLoop(point1.lng, -180, 180);
         point1.lat = _getRange(point1.lat, -74, 74);
         point2.lng = _getLoop(point2.lng, -180, 180);
         point2.lat = _getRange(point2.lat, -74, 74);
         
         var x1, x2, y1, y2;
         x1 = GeoUtils.degreeToRad(point1.lng);
         y1 = GeoUtils.degreeToRad(point1.lat);
         x2 = GeoUtils.degreeToRad(point2.lng);
         y2 = GeoUtils.degreeToRad(point2.lat);
 
         return EARTHRADIUS * Math.acos((Math.sin(y1) * Math.sin(y2) + Math.cos(y1) * Math.cos(y2) * Math.cos(x2 - x1)));    
     }
     
     /**
      * 計算折線或者點數組的長度
      * @param {Polyline|Array<Point>} polyline 折線對象或者點數組
      * @returns {Number} 折線或點數組對應的長度
      */
     GeoUtils.getPolylineDistance = function(polyline){
         //檢查類型
         if(polyline instanceof BMap.Polyline || 
             polyline instanceof Array){
             //將polyline統一爲數組
             var pts;
             if(polyline instanceof BMap.Polyline){
                 pts = polyline.getPath();
             } else {
                 pts = polyline;
             }
             
             if(pts.length < 2){//小於2個點，返回0
                 return 0;
             }
 
             //遍歷全部線段將其相加，計算整條線段的長度
             var totalDis = 0;
             for(var i =0; i < pts.length - 1; i++){
                 var curPt = pts[i];
                 var nextPt = pts[i + 1]
                 var dis = GeoUtils.getDistance(curPt, nextPt);
                 totalDis += dis;
             }
 
             return totalDis;
             
         } else {
             return 0;
         }
     }
     
     /**
      * 計算多邊形面或點數組構建圖形的面積,注意：座標類型只能是經緯度，且不適合計算自相交多邊形的面積
      * @param {Polygon|Array<Point>} polygon 多邊形面對象或者點數組
      * @returns {Number} 多邊形面或點數組構成圖形的面積
      */
     GeoUtils.getPolygonArea = function(polygon){
         //檢查類型
         if(!(polygon instanceof BMap.Polygon) &&
             !(polygon instanceof Array)){
             return 0;
         }
         var pts;
         if(polygon instanceof BMap.Polygon){
             pts = polygon.getPath();
         }else{
             pts = polygon;    
         }
         
         if(pts.length < 3){//小於3個頂點，不能構建面
             return 0;
         }
         
         var totalArea = 0;//初始化總面積
         var LowX = 0.0;
         var LowY = 0.0;
         var MiddleX = 0.0;
         var MiddleY = 0.0;
         var HighX = 0.0;
         var HighY = 0.0;
         var AM = 0.0;
         var BM = 0.0;
         var CM = 0.0;
         var AL = 0.0;
         var BL = 0.0;
         var CL = 0.0;
         var AH = 0.0;
         var BH = 0.0;
         var CH = 0.0;
         var CoefficientL = 0.0;
         var CoefficientH = 0.0;
         var ALtangent = 0.0;
         var BLtangent = 0.0;
         var CLtangent = 0.0;
         var AHtangent = 0.0;
         var BHtangent = 0.0;
         var CHtangent = 0.0;
         var ANormalLine = 0.0;
         var BNormalLine = 0.0;
         var CNormalLine = 0.0;
         var OrientationValue = 0.0;
         var AngleCos = 0.0;
         var Sum1 = 0.0;
         var Sum2 = 0.0;
         var Count2 = 0;
         var Count1 = 0;
         var Sum = 0.0;
         var Radius = EARTHRADIUS; //6378137.0,WGS84橢球半徑 
         var Count = pts.length;        
         for (var i = 0; i < Count; i++) {
             if (i == 0) {
                 LowX = pts[Count - 1].lng * Math.PI / 180;
                 LowY = pts[Count - 1].lat * Math.PI / 180;
                 MiddleX = pts[0].lng * Math.PI / 180;
                 MiddleY = pts[0].lat * Math.PI / 180;
                 HighX = pts[1].lng * Math.PI / 180;
                 HighY = pts[1].lat * Math.PI / 180;
             }
             else if (i == Count - 1) {
                 LowX = pts[Count - 2].lng * Math.PI / 180;
                 LowY = pts[Count - 2].lat * Math.PI / 180;
                 MiddleX = pts[Count - 1].lng * Math.PI / 180;
                 MiddleY = pts[Count - 1].lat * Math.PI / 180;
                 HighX = pts[0].lng * Math.PI / 180;
                 HighY = pts[0].lat * Math.PI / 180;
             }
             else {
                 LowX = pts[i - 1].lng * Math.PI / 180;
                 LowY = pts[i - 1].lat * Math.PI / 180;
                 MiddleX = pts[i].lng * Math.PI / 180;
                 MiddleY = pts[i].lat * Math.PI / 180;
                 HighX = pts[i + 1].lng * Math.PI / 180;
                 HighY = pts[i + 1].lat * Math.PI / 180;
             }
             AM = Math.cos(MiddleY) * Math.cos(MiddleX);
             BM = Math.cos(MiddleY) * Math.sin(MiddleX);
             CM = Math.sin(MiddleY);
             AL = Math.cos(LowY) * Math.cos(LowX);
             BL = Math.cos(LowY) * Math.sin(LowX);
             CL = Math.sin(LowY);
             AH = Math.cos(HighY) * Math.cos(HighX);
             BH = Math.cos(HighY) * Math.sin(HighX);
             CH = Math.sin(HighY);
             CoefficientL = (AM * AM + BM * BM + CM * CM) / (AM * AL + BM * BL + CM * CL);
             CoefficientH = (AM * AM + BM * BM + CM * CM) / (AM * AH + BM * BH + CM * CH);
             ALtangent = CoefficientL * AL - AM;
             BLtangent = CoefficientL * BL - BM;
             CLtangent = CoefficientL * CL - CM;
             AHtangent = CoefficientH * AH - AM;
             BHtangent = CoefficientH * BH - BM;
             CHtangent = CoefficientH * CH - CM;
             AngleCos = (AHtangent * ALtangent + BHtangent * BLtangent + CHtangent * CLtangent) / (Math.sqrt(AHtangent * AHtangent + BHtangent * BHtangent + CHtangent * CHtangent) * Math.sqrt(ALtangent * ALtangent + BLtangent * BLtangent + CLtangent * CLtangent));
             AngleCos = Math.acos(AngleCos);            
             ANormalLine = BHtangent * CLtangent - CHtangent * BLtangent;
             BNormalLine = 0 - (AHtangent * CLtangent - CHtangent * ALtangent);
             CNormalLine = AHtangent * BLtangent - BHtangent * ALtangent;
             if (AM != 0)
                 OrientationValue = ANormalLine / AM;
             else if (BM != 0)
                 OrientationValue = BNormalLine / BM;
             else
                 OrientationValue = CNormalLine / CM;
             if (OrientationValue > 0) {
                 Sum1 += AngleCos;
                 Count1++;
             }
             else {
                 Sum2 += AngleCos;
                 Count2++;
             }
         }        
         var tempSum1, tempSum2;
         tempSum1 = Sum1 + (2 * Math.PI * Count2 - Sum2);
         tempSum2 = (2 * Math.PI * Count1 - Sum1) + Sum2;
         if (Sum1 > Sum2) {
             if ((tempSum1 - (Count - 2) * Math.PI) < 1)
                 Sum = tempSum1;
             else
                 Sum = tempSum2;
         }
         else {
             if ((tempSum2 - (Count - 2) * Math.PI) < 1)
                 Sum = tempSum2;
             else
                 Sum = tempSum1;
         }
         totalArea = (Sum - (Count - 2) * Math.PI) * Radius * Radius;
 
         return totalArea; //返回總面積
     }
    
 })();//閉包結束

 參考文獻：http://erich.realtimerendering.com/ptinpoly/

博客地址：	https://www.cnblogs.com/smartstar/
博客版權：	本文以學習、研究和分享爲主，歡迎轉載，但必須在文章頁面明顯位置給出原文鏈接。若是文中有不妥或者錯誤的地方還望高手的你指出，以避免誤人子弟。若是以爲本文對你有所幫助不如【推薦】一下！若是你有更好的建議，不如留言一塊兒討論，共同進步！再次感謝您耐心的讀完本篇文章。