判斷點是否在多邊形內的Python實現及小應用（射線法）

判斷一個點是否在多邊形內是處理空間資料時經常面對的需求，例如GIS中的點選功能、根據多邊形邊界篩選出位於多邊形內的點、求交集、篩選不在多邊形內的點等等。
判斷一個點是否在多邊形內有幾種不同的思路，相應的方法（感覺還談不上演算法）有：

• 射線法：從判斷點向某個統一方向作射線，依交點個數的奇偶判斷；
• 轉角法：按照多邊形頂點逆時針順序，根據頂點和判斷點連線的方向正負（設定角度逆時針為正）求和判斷；
• 夾角和法：求判斷點與所有邊的夾角和，等於360度則在多邊形內部。
• 面積和法：求判斷點與多邊形邊組成的三角形面積和，等於多邊形面積則點在多邊形內部。

面積和法涉及多個面積的計算，比較複雜，夾角和法以及轉角法用到角度計算，會涉及反三角函式，計算開銷比較大，而射線法主要涉及迴圈多邊形的每條邊進行求交運算，但大部分邊可以通過簡單座標比對直接排除，因此這是比較好的方法。

射線法的實現

射線法就是以判斷點開始，向右（或向左）的水平方向作一射線，計算該射線與多邊形每條邊的交點個數，如果交點個數為奇數，則點位於多邊形內，偶數則在多邊形外。該演算法對於複合多邊形也能正確判斷。

射線法的關鍵是正確計算射線與每條邊是否相交。並且規定線段與射線重疊或者射線經過線段下端點屬於不相交。首先排除掉不相交的情況，下圖的情況都是需要排除掉的：

排除掉這些情況的函式如下：

def isRayIntersectsSegment(poi,s_poi,e_poi): #[x,y] [lng,lat]
    #輸入：判斷點，邊起點，邊終點，都是[lng,lat]格式陣列
    if s_poi[1]==e_poi[1]: #排除與射線平行、重合，線段首尾端點重合的情況
        return False
    if s_poi[1]>poi[1] and e_poi[1]>poi[1]: #線段在射線上邊
        return False
    if s_poi[1]<poi[1] and e_poi[1]<poi[1]: #線段在射線下邊
        return False
    if s_poi[1]==poi[1] and e_poi[1]>poi[1]: #交點為下端點，對應spoint
        return False
    if e_poi[1]==poi[1] and s_poi[1]>poi[1]: #交點為下端點，對應epoint
        return False
    if s_poi[0]<poi[0] and e_poi[1]<poi[1]: #線段在射線左邊
        return False

    xseg=e_poi[0]-(e_poi[0]-s_poi[0])*(e_poi[1]-poi[1])/(e_poi[1]-s_poi[1]) #求交
    if xseg<poi[0]: #交點在射線起點的左側
        return False
    return True  #排除上述情況之後

排除掉上述情況真正需要求交點來判斷的情況只有兩種：

函式isRayIntersectsSegment()裡求交的部分就是利用兩個三角形的比例關係求出交點在起點的左邊還是右邊；用圖去理解如下：

最後判斷的程式碼如下：

def isPoiWithinPoly(poi,poly):
    #輸入：點，多邊形三維陣列
    #poly=[[[x1,y1],[x2,y2],……,[xn,yn],[x1,y1]],[[w1,t1],……[wk,tk]]] 三維陣列

    #可以先判斷點是否在外包矩形內 
    #if not isPoiWithinBox(poi,mbr=[[0,0],[180,90]]): return False
    #但算最小外包矩形本身需要迴圈邊，會造成開銷，本處略去
    sinsc=0 #交點個數
    for epoly in poly: #迴圈每條邊的曲線->each polygon 是二維陣列[[x1,y1],…[xn,yn]]
        for i in range(len(epoly)-1): #[0,len-1]
            s_poi=epoly[i]
            e_poi=epoly[i+1]
            if isRayIntersectsSegment(poi,s_poi,e_poi):
                sinsc+=1 #有交點就加1

    return True if sinsc%2==1 else  False

我們取一個比較複雜的多邊形進行測試，多邊形和一些點如圖：

用isPoiWithinPoly()的測試結果如下：

點在多邊形內的應用

上面第一段已經描述了一些應用場景，下面給出一個應用的例子：有一堆點資料存在csv檔案裡，如何檢索位於某個城市的點出來，檢索出來之後的分析（例如加標籤、改屬性、做統計還是其他）這裡不討論，檢索的結果統一寫到新檔案裡。點輸入的格式如下：

id,name,wgslng,wgslat,score,adds
1,沃美,116.3309,40.0706,4.3,昌平回龍觀同成街華聯購物中心4樓
2,星美國際,116.446,39.916,5,金匯路8號世界城E座
3,……

城市邊界為geojson格式，就是加了一些限定條件的json格式資料，如果需要詳細瞭解geojson格式，可以參考本人之前的文章：GEOJSON標準格式學習。形如：

{
  "type": "FeatureCollection",
  "features": [{
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Polygon",
        "coordinates":
         [
            [
              [108.71658325195312,34.231106222010531],
              [108.96240234375,34.168635904722734],
              [109.00222778320313,34.354774165387568],
              [108.80172729492186,34.35023911062779],
              [108.71658325195312,34.231106222010531]
            ]
          ]
        }
      }
  ]
}

下面的程式碼只考慮了Polygon的情況，對於MultiPolygon也是比較容易改的，要改為處理kml儲存的邊界資料也不難改。文中程式碼同步於本人GitHub。

import json
import csv
def pointInPolygon():
    gfile = './beijing_poly_wgs84.geojson' #utf-8編碼
    cin_path = './poi_cinema_wgs84.csv'
    out_path = './beijing_poi_cinema_wgs84.csv' #輸出檔案

    pindex = [2, 3]  # wgslng,wgslat 在的位置

    with open(out_path, 'w', newline='') as cout_file:
        fin = open(cin_path, 'r', encoding='gbk') #出現編碼錯誤就改編碼 utf-8
        gfn = open(gfile, 'r', encoding='utf-8')
        gjson = json.load(gfn)
        polygon = gjson["features"][0]["geometry"]['coordinates'] #提取多邊形,如果是4維陣列需要相應的處理
        filewriter = csv.writer(cout_file, delimiter=',')
        w = 0
        for line in csv.reader(fin, delimiter=','):
            if w == 0: #寫入表頭 id,name,… 如果沒有就去掉if語句
                filewriter.writerow(line)
                w = 1
                continue
            point = [float(line[pindex[0]]), float(line[pindex][1])]
            if isPoiWithinPoly(point, polygon): #在多邊形內，寫入新表
                filewriter.writerow(line)
            else:
                continue
        fin.close()
        gfn.close()
    print('done')