判斷一個(gè)點(diǎn)是否在多邊形內(nèi)是處理空間數(shù)據(jù)時(shí)經(jīng)常面對的需求，例如GIS中的點(diǎn)選功能初澎、根據(jù)多邊形邊界篩選出位于多邊形內(nèi)的點(diǎn)举农、求交集、篩選不在多邊形內(nèi)的點(diǎn)等等腌且。
判斷一個(gè)點(diǎn)是否在多邊形內(nèi)有幾種不同的思路梗肝，相應(yīng)的方法（感覺還談不上算法）有：
- 射線法：從判斷點(diǎn)向某個(gè)統(tǒng)一方向作射線，依交點(diǎn)個(gè)數(shù)的奇偶判斷铺董；
- 轉(zhuǎn)角法：按照多邊形頂點(diǎn)逆時(shí)針順序巫击，根據(jù)頂點(diǎn)和判斷點(diǎn)連線的方向正負(fù)（設(shè)定角度逆時(shí)針為正）求和判斷；
- 夾角和法：求判斷點(diǎn)與所有邊的夾角和，等于360度則在多邊形內(nèi)部喘鸟。
- 面積和法：求判斷點(diǎn)與多邊形邊組成的三角形面積和匆绣，等于多邊形面積則點(diǎn)在多邊形內(nèi)部。
面積和法涉及多個(gè)面積的計(jì)算什黑，比較復(fù)雜崎淳，夾角和法以及轉(zhuǎn)角法用到角度計(jì)算，會(huì)涉及反三角函數(shù)愕把，計(jì)算開銷比較大拣凹，而射線法主要涉及循環(huán)多邊形的每條邊進(jìn)行求交運(yùn)算，但大部分邊可以通過簡單坐標(biāo)比對直接排除恨豁，因此這是比較好的方法嚣镜。
射線法的實(shí)現(xiàn)
射線法就是以判斷點(diǎn)開始，向右（或向左）的水平方向作一射線橘蜜，計(jì)算該射線與多邊形每條邊的交點(diǎn)個(gè)數(shù)菊匿，如果交點(diǎn)個(gè)數(shù)為奇數(shù)，則點(diǎn)位于多邊形內(nèi)计福，偶數(shù)則在多邊形外跌捆。該算法對于復(fù)合多邊形也能正確判斷。
射線法的關(guān)鍵是正確計(jì)算射線與每條邊是否相交象颖。并且規(guī)定線段與射線重疊或者射線經(jīng)過線段下端點(diǎn)屬于不相交佩厚。首先排除掉不相交的情況，下圖的情況都是需要排除掉的：
排除掉這些情況的函數(shù)如下：
def isRayIntersectsSegment(poi,s_poi,e_poi): #[x,y] [lng,lat]
#輸入：判斷點(diǎn)说订，邊起點(diǎn)抄瓦，邊終點(diǎn)，都是[lng,lat]格式數(shù)組
if s_poi[1]==e_poi[1]: #排除與射線平行陶冷、重合钙姊，線段首尾端點(diǎn)重合的情況
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]: #交點(diǎn)為下端點(diǎn)，對應(yīng)spoint
return False
if e_poi[1]==poi[1] and s_poi[1]>poi[1]: #交點(diǎn)為下端點(diǎn)埃叭，對應(yīng)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]: #交點(diǎn)在射線起點(diǎn)的左側(cè)
return False
return True #排除上述情況之后
排除掉上述情況真正需要求交點(diǎn)來判斷的情況只有兩種：
函數(shù)isRayIntersectsSegment()里求交的部分就是利用兩個(gè)三角形的比例關(guān)系求出交點(diǎn)在起點(diǎn)的左邊還是右邊摸恍；用圖去理解如下：
最后判斷的代碼如下：
def isPoiWithinPoly(poi,poly):
#輸入：點(diǎn)，多邊形三維數(shù)組
#poly=[[[x1,y1],[x2,y2],……,[xn,yn],[x1,y1]],[[w1,t1],……[wk,tk]]] 三維數(shù)組
#可以先判斷點(diǎn)是否在外包矩形內(nèi)
#if not isPoiWithinBox(poi,mbr=[[0,0],[180,90]]): return False
#但算最小外包矩形本身需要循環(huán)邊赤屋，會(huì)造成開銷，本處略去
sinsc=0 #交點(diǎn)個(gè)數(shù)
for epoly in poly: #循環(huán)每條邊的曲線->each polygon 是二維數(shù)組[[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 #有交點(diǎn)就加1
return True if sinsc%2==1 else False
我們?nèi)∫粋€(gè)比較復(fù)雜的多邊形進(jìn)行測試壁袄，多邊形和一些點(diǎn)如圖：
用isPoiWithinPoly()的測試結(jié)果如下：
點(diǎn)在多邊形內(nèi)的應(yīng)用
上面第一段已經(jīng)描述了一些應(yīng)用場景类早，下面給出一個(gè)應(yīng)用的例子：有一堆點(diǎn)數(shù)據(jù)存在csv文件里，如何檢索位于某個(gè)城市的點(diǎn)出來嗜逻，檢索出來之后的分析（例如加標(biāo)簽涩僻、改屬性、做統(tǒng)計(jì)還是其他）這里不討論，檢索的結(jié)果統(tǒng)一寫到新文件里逆日。點(diǎn)輸入的格式如下：
id,name,wgslng,wgslat,score,adds
1,沃美,116.3309,40.0706,4.3,昌平回龍觀同成街華聯(lián)購物中心4樓
2,星美國際,116.446,39.916,5,金匯路8號(hào)世界城E座
3,……
城市邊界為geojson格式嵌巷，就是加了一些限定條件的json格式數(shù)據(jù)，如果需要詳細(xì)了解geojson格式室抽，可以參考本人之前的文章：GEOJSON標(biāo)準(zhǔn)格式學(xué)習(xí)搪哪。形如：
{
"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保存的邊界數(shù)據(jù)也不難改晓折。文中代碼同步于本人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') #出現(xiàn)編碼錯(cuò)誤就改編碼 utf-8
gfn = open(gfile, 'r', encoding='utf-8')
gjson = json.load(gfn)
polygon = gjson["features"][0]["geometry"]['coordinates'] #提取多邊形,如果是4維數(shù)組需要相應(yīng)的處理
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): #在多邊形內(nèi)兽泄，寫入新表
filewriter.writerow(line)
else:
continue
fin.close()
gfn.close()
print('done')
歡迎關(guān)注本人公眾號(hào)漓概，有更多有趣內(nèi)容和資料：
