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

  • 2019 年 10 月 10 日
  • 筆記

判斷一個點是否在多邊形內是處理空間數據時經常面對的需求，例如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保存的邊界數據也不難改。文中代碼同步於https://github.com/ QLWeilcf/ LcfGeoProject/blob/ master/poiWithinPolygon.py 中。

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')

VirMach 便宜 VPS

QNews

QNews