如何判斷一個指定的位置點座標（GPS上的經緯度）是否落在一個多邊形區域內？

天府雲創發表於2018-06-11

業務場景舉例：快遞選擇收穫區域、車輛電子圍欄、運動軌跡路線、地理位置資訊檢測範圍和地圖等過濾等等。

比方說地圖上有一塊區域（抽象成多邊形），然後裡面每一個位置點(畫素點)都有對應的GPS的經緯度座標值，題目要求的就是判斷任意點（使用者輸入的資訊）與多邊形的位置關係（是否在裡面還是在圖形區域外面）。


具體有一個需求為：每一個店維護多個可配送的地址，配送地址為地圖中的多邊形區域，使用者選擇收貨地址的時候需要判斷該收貨地址在不在多邊形區域內。（給定一個點的座標以及一個多邊形的所有頂點座標。要求能夠判斷這個點是在多邊形內，還是在多邊形外？ ）


驗證地址：Map Polygon/Polyline Tool https://www.keene.edu/campus/maps/tool/


以上的ABCDE,分別是以下陣列裡面的資料

[java] view plain copy
  1. Point[] ps = new Point[] { new Point(120.2043 , 30.2795), new Point(120.2030 , 30.2511), new Point(120.1810 , 30.2543), new Point(120.1798 , 30.2781), new Point(120.1926,30.2752) };  

那麼問題來了：如何判斷一個指定的經緯度點是否落在一個多邊形區域內？

這個是演算法和判斷的依據也是解答本題的關鍵。（保證準確率和速度的關鍵）

網路上有很多演算法和第三庫來實現這類功能，但本文著重講原理和自己實現寫程式。

話不多說，直接丟擲寫程式通常用的流程模型圖，如下；


經過在網上的一番搜尋，

發現目前比較通用的就是射線法和座標軸法，而我採用的就是X軸射線法。

主要理論來源於西安交大的一篇論文（即參考文獻的第二條） 

判斷一個點向左的射線跟一個多邊形的交叉點有幾個，如果結果為奇數的話那麼說明這個點落在多邊形中，反之則不在。

1、理論支援：如果從需要判斷的點出發的一條射線與該多邊形的焦點個數為奇數，則該點在此多邊形內，否則該點在此多邊形外。（射線不能與多邊形頂點相交）
2、程式設計思路：
該程式的思路是從A點出發向左做一條水平射線（平行於x軸，向X軸的反方向），判斷與各邊是否有焦點。
dLon1, dLon2, dLat1, dLat2分別表示邊的起點和終點的經度和緯度（x軸和y軸）。
先判斷A點是否在邊的兩端點d1和d2的水平平行線之間，不在則不可能有交點，繼續判斷下一條邊。
在之間則說明可能與A點向左的射線有交點，接下來利用幾何方法得到A點的水平直線與該邊交點的x座標。
然後判斷交點的x座標在A點的左側還是右側，左側則總交點數加一，右側則不在A點左射線上，繼續判斷下一條邊。
程式碼講解： 主要的類有兩個↓↓

一個是座標點的抽象類（點的座標位置），另一個是位置關係判斷工具類（判斷點和多邊形的關係）。 

因為座標描點要涉及到座標以及小數點的經緯關係，故要用到浮點型運算，也就是要使用到double雙精度。而地圖涉及到很多個畫素點構成的圖形區域，故要用到list陣列。結果要展示需要用到js程式碼。

def IsPtInPoly(aLon, aLat, pointList):  
    ''''' 
    :param aLon: double 經度 
    :param aLat: double 緯度 
    :param pointList: list [(lon, lat)...] 多邊形點的順序需根據順時針或逆時針，不能亂 
    '''  
      
    iSum = 0  
    iCount = len(pointList)  
      
    if(iCount < 3):  
        return False  
      
      
    for i in range(iCount):  
          
        pLon1 = pointList[i][0]  
        pLat1 = pointList[i][1]  
          
        if(i == iCount - 1):  
              
            pLon2 = pointList[0][0]  
            pLat2 = pointList[0][1]  
        else:  
            pLon2 = pointList[i + 1][0]  
            pLat2 = pointList[i + 1][1]  
          
        if ((aLat >= pLat1) and (aLat < pLat2)) or ((aLat>=pLat2) and (aLat < pLat1)):  
              
            if (abs(pLat1 - pLat2) > 0):  
                  
                pLon = pLon1 - ((pLon1 - pLon2) * (pLat1 - aLat)) / (pLat1 - pLat2);  
                  
                if(pLon < aLon):  
                    iSum += 1  
  
    if(iSum % 2 != 0):  
        return True  
    else:  
        return False
1、是座標點的抽象類 
Java程式碼  
  1. package com.niux.crm.core.common.bmap;  
  2.   
  3. /** 
  4.  * 用於構造百度地圖中的經緯度點 
  5.  *  
  6.  * @author zhengtian 
  7.  * @date 2013-8-5 下午02:54:41 
  8.  */  
  9. public class BmapPoint {  
  10.     private double lng;// 經度  
  11.     private double lat;// 緯度  
  12.   
  13.     public BmapPoint() {  
  14.   
  15.     }  
  16.   
  17.     public BmapPoint(double lng, double lat) {  
  18.         this.lng = lng;  
  19.         this.lat = lat;  
  20.     }  
  21.   
  22.     @Override  
  23.     public boolean equals(Object obj) {  
  24.         if (obj instanceof BmapPoint) {  
  25.             BmapPoint bmapPoint = (BmapPoint) obj;  
  26.             return (bmapPoint.getLng() == lng && bmapPoint.getLat() == lat) ? true : false;  
  27.         } else {  
  28.             return false;  
  29.         }  
  30.     }  
  31.   
  32.     public double getLng() {  
  33.         return lng;  
  34.     }  
  35.   
  36.     public void setLng(double lng) {  
  37.         this.lng = lng;  
  38.     }  
  39.   
  40.     public double getLat() {  
  41.         return lat;  
  42.     }  
  43.   
  44.     public void setLat(double lat) {  
  45.         this.lat = lat;  
  46.     }  
  47. }  


2、位置關係判斷工具類 

點與多邊形的位置關係的判定規則：1、，根據多邊形的座標，虛擬出一個外包矩形，主要是為了提前過濾不相關的點，減少運算量。2、然後判斷是否有重合的點。3、判斷點與斜線的交點。4、判斷點過頂點的情況。5、判斷點與邊重合的情況。6、判斷點在邊上的情況。 

其中點過頂點，以及點與邊重合的情況，主要採用了加權邊的思想，論文與程式碼中有註釋。 

Java程式碼  
  1. package com.niux.crm.core.common.bmap;  
  2.   
  3. import java.util.Arrays;  
  4.   
  5.   
  6. /** 
  7.  * 用於點與多邊形位置關係的判斷 
  8.  *  
  9.  * @author zhengtian 
  10.  * @date 2013-8-5 上午11:59:35 
  11.  */  
  12. public class GraphUtils {  
  13.   
  14.     /** 
  15.      * 判斷點是否在多邊形內(基本思路是用交點法) 
  16.      *  
  17.      * @param point 
  18.      * @param boundaryPoints 
  19.      * @return 
  20.      */  
  21.     public static boolean isPointInPolygon(BmapPoint point, BmapPoint[] boundaryPoints) {  
  22.         // 防止第一個點與最後一個點相同  
  23.         if (boundaryPoints != null && boundaryPoints.length > 0  
  24.                 && boundaryPoints[boundaryPoints.length - 1].equals(boundaryPoints[0])) {  
  25.             boundaryPoints = Arrays.copyOf(boundaryPoints, boundaryPoints.length - 1);  
  26.         }  
  27.         int pointCount = boundaryPoints.length;  
  28.   
  29.         // 首先判斷點是否在多邊形的外包矩形內，如果在，則進一步判斷，否則返回false  
  30.         if (!isPointInRectangle(point, boundaryPoints)) {  
  31.             return false;  
  32.         }  
  33.   
  34.         // 如果點與多邊形的其中一個頂點重合，那麼直接返回true  
  35.         for (int i = 0; i < pointCount; i++) {  
  36.             if (point.equals(boundaryPoints[i])) {  
  37.                 return true;  
  38.             }  
  39.         }  
  40.   
  41.         /** 
  42.          * 基本思想是利用X軸射線法，計算射線與多邊形各邊的交點，如果是偶數，則點在多邊形外，否則在多邊形內。還會考慮一些特殊情況，如點在多邊形頂點上 
  43.          * ， 點在多邊形邊上等特殊情況。 
  44.          */  
  45.         // X軸射線與多邊形的交點數  
  46.         int intersectPointCount = 0;  
  47.         // X軸射線與多邊形的交點權值  
  48.         float intersectPointWeights = 0;  
  49.         // 浮點型別計算時候與0比較時候的容差  
  50.         double precision = 2e-10;  
  51.         // 邊P1P2的兩個端點  
  52.         BmapPoint point1 = boundaryPoints[0], point2;  
  53.         // 迴圈判斷所有的邊  
  54.         for (int i = 1; i <= pointCount; i++) {  
  55.             point2 = boundaryPoints[i % pointCount];  
  56.   
  57.             /** 
  58.              * 如果點的y座標在邊P1P2的y座標開區間範圍之外，那麼不相交。 
  59.              */  
  60.             if (point.getLat() < Math.min(point1.getLat(), point2.getLat())  
  61.                     || point.getLat() > Math.max(point1.getLat(), point2.getLat())) {  
  62.                 point1 = point2;  
  63.                 continue;  
  64.             }  
  65.   
  66.             /** 
  67.              * 此處判斷射線與邊相交 
  68.              */  
  69.             if (point.getLat() > Math.min(point1.getLat(), point2.getLat())  
  70.                     && point.getLat() < Math.max(point1.getLat(), point2.getLat())) {// 如果點的y座標在邊P1P2的y座標開區間內  
  71.                 if (point1.getLng() == point2.getLng()) {// 若邊P1P2是垂直的  
  72.                     if (point.getLng() == point1.getLng()) {  
  73.                         // 若點在垂直的邊P1P2上，則點在多邊形內  
  74.                         return true;  
  75.                     } else if (point.getLng() < point1.getLng()) {  
  76.                         // 若點在在垂直的邊P1P2左邊，則點與該邊必然有交點  
  77.                         ++intersectPointCount;  
  78.                     }  
  79.                 } else {// 若邊P1P2是斜線  
  80.                     if (point.getLng() <= Math.min(point1.getLng(), point2.getLng())) {// 點point的x座標在點P1和P2的左側  
  81.                         ++intersectPointCount;  
  82.                     } else if (point.getLng() > Math.min(point1.getLng(), point2.getLng())  
  83.                             && point.getLng() < Math.max(point1.getLng(), point2.getLng())) {// 點point的x座標在點P1和P2的x座標中間  
  84.                         double slopeDiff = 0.0d;  
  85.                         if (point1.getLat() > point2.getLat()) {  
  86.                             slopeDiff = (point.getLat() - point2.getLat()) / (point.getLng() - point2.getLng())  
  87.                                     - (point1.getLat() - point2.getLat()) / (point1.getLng() - point2.getLng());  
  88.                         } else {  
  89.                             slopeDiff = (point.getLat() - point1.getLat()) / (point.getLng() - point1.getLng())  
  90.                                     - (point2.getLat() - point1.getLat()) / (point2.getLng() - point1.getLng());  
  91.                         }  
  92.                         if (slopeDiff > 0) {  
  93.                             if (slopeDiff < precision) {// 由於double精度在計算時會有損失，故匹配一定的容差。經試驗，座標經度可以達到0.0001  
  94.                                 // 點在斜線P1P2上  
  95.                                 return true;  
  96.                             } else {  
  97.                                 // 點與斜線P1P2有交點  
  98.                                 intersectPointCount++;  
  99.                             }  
  100.                         }  
  101.                     }  
  102.                 }  
  103.             } else {  
  104.                 // 邊P1P2水平  
  105.                 if (point1.getLat() == point2.getLat()) {  
  106.                     if (point.getLng() <= Math.max(point1.getLng(), point2.getLng())  
  107.                             && point.getLng() >= Math.min(point1.getLng(), point2.getLng())) {  
  108.                         // 若點在水平的邊P1P2上，則點在多邊形內  
  109.                         return true;  
  110.                     }  
  111.                 }  
  112.                 /** 
  113.                  * 判斷點通過多邊形頂點 
  114.                  */  
  115.                 if (((point.getLat() == point1.getLat() && point.getLng() < point1.getLng()))  
  116.                         || (point.getLat() == point2.getLat() && point.getLng() < point2.getLng())) {  
  117.                     if (point2.getLat() < point1.getLat()) {  
  118.                         intersectPointWeights += -0.5;  
  119.                     } else if (point2.getLat() > point1.getLat()) {  
  120.                         intersectPointWeights += 0.5;  
  121.                     }  
  122.                 }  
  123.             }  
  124.             point1 = point2;  
  125.         }  
  126.   
  127.         if ((intersectPointCount + Math.abs(intersectPointWeights)) % 2 == 0) {// 偶數在多邊形外  
  128.             return false;  
  129.         } else { // 奇數在多邊形內  
  130.             return true;  
  131.         }  
  132.     }  
  133.   
  134.     /** 
  135.      * 判斷點是否在矩形內在矩形邊界上，也算在矩形內(根據這些點，構造一個外包矩形) 
  136.      *  
  137.      * @param point 
  138.      *            點物件 
  139.      * @param boundaryPoints 
  140.      *            矩形邊界點 
  141.      * @return 
  142.      */  
  143.     public static boolean isPointInRectangle(BmapPoint point, BmapPoint[] boundaryPoints) {  
  144.         BmapPoint southWestPoint = getSouthWestPoint(boundaryPoints); // 西南角點  
  145.         BmapPoint northEastPoint = getNorthEastPoint(boundaryPoints); // 東北角點  
  146.         return (point.getLng() >= southWestPoint.getLng() && point.getLng() <= northEastPoint.getLng()  
  147.                 && point.getLat() >= southWestPoint.getLat() && point.getLat() <= northEastPoint.getLat());  
  148.   
  149.     }  
  150.   
  151.     /** 
  152.      * 根據這組座標，畫一個矩形，然後得到這個矩形西南角的頂點座標 
  153.      *  
  154.      * @param vertexs 
  155.      * @return 
  156.      */  
  157.     private static BmapPoint getSouthWestPoint(BmapPoint[] vertexs) {  
  158.         double minLng = vertexs[0].getLng(), minLat = vertexs[0].getLat();  
  159.         for (BmapPoint bmapPoint : vertexs) {  
  160.             double lng = bmapPoint.getLng();  
  161.             double lat = bmapPoint.getLat();  
  162.             if (lng < minLng) {  
  163.                 minLng = lng;  
  164.             }  
  165.             if (lat < minLat) {  
  166.                 minLat = lat;  
  167.             }  
  168.         }  
  169.         return new BmapPoint(minLng, minLat);  
  170.     }  
  171.   
  172.     /** 
  173.      * 根據這組座標，畫一個矩形，然後得到這個矩形東北角的頂點座標 
  174.      *  
  175.      * @param vertexs 
  176.      * @return 
  177.      */  
  178.     private static BmapPoint getNorthEastPoint(BmapPoint[] vertexs) {  
  179.         double maxLng = 0.0d, maxLat = 0.0d;  
  180.         for (BmapPoint bmapPoint : vertexs) {  
  181.             double lng = bmapPoint.getLng();  
  182.             double lat = bmapPoint.getLat();  
  183.             if (lng > maxLng) {  
  184.                 maxLng = lng;  
  185.             }  
  186.             if (lat > maxLat) {  
  187.                 maxLat = lat;  
  188.             }  
  189.         }  
  190.         return new BmapPoint(maxLng, maxLat);  
  191.     }  
  192.   
  193. }  

【關於3D地圖圖形請參考】

IOS高德3D地圖畫多邊形，以及判斷某一經緯度是否在該多邊形內 - 簡書 https://www.jianshu.com/p/fb1177cac1ec

翻看了一下高德提供的技術文件，沒找到3D地圖下該怎麼完成此功能。無奈自己翻找高德的SDK，發現了這兩個類：MAPolygon MAOverlayRenderer

MAPolygon用於定義一個由多個點組成的閉合多邊形,點與點之間按順序尾部相連,第一個點與最後一個點相連,通常MAPolygon是MAPolygonView的model

MAOverlayRenderer是地圖覆蓋物Renderer的基類,提供繪製overlay的介面但並無實際的實現（render相關方法只能在重寫後的glRender方法中使用）

【參考文獻】
1、兩條直線的關係 
http://www.cnblogs.com/devymex/archive/2010/08/19/1803885.html 
2、點與多邊形的關係 

http://wenku.baidu.com/view/5e3913a2b0717fd5360cdccf.html?qq-pf-to=pcqq.c2c 

3、https://en.wikipedia.org/wiki/Haversine_formula

4、百度map js程式碼 https://www.cnblogs.com/relax/p/3507014.html

5、關於經緯度得到的多邊形面積。 https://www.cnblogs.com/JeffController/p/5618742.html

【附錄座標軸函式公式圖】

直接上程式碼（兩個點）半正矢公式 計算（Haversine formula）：因為難度大暫時不考慮此類演算法。


關於指定的經緯度是否落在多邊形內 - https://blog.csdn.net/qq_22929803/article/details/46818009



相關文章