#coding=UTF-8
import csv
import json
# 点是否在外包矩形内
def isPoiWithinBox(poi, sbox, toler=0.0001):
# sbox=[[x1,y1],[x2,y2]]
# 不考虑在边界上，需要考虑就加等号
if poi[0] > sbox[0][0] and poi[0] < sbox[1][0] and poi[1] > sbox[0][1] and poi[1] < sbox[1][1]:
return True
if toler > 0:
pass
return False
# 射线与边是否有交点
def isRayIntersectsSegment(poi, s_poi, e_poi): # [x,y] [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]:
return False
if e_poi[1] == poi[1] and s_poi[1] > poi[1]:
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
#只适用简单多边形
def isPoiWithinSimplePoly(poi, simPoly, tolerance=0.0001):
# 点；多边形数组；容限
# simPoly=[[x1,y1],[x2,y2],……,[xn,yn],[x1,y1]]
# 如果simPoly=[[x1,y1],[x2,y2],……,[xn,yn]] i循环到终点后还需要判断[xn,yx]和[x1,y1]
# 先判断点是否在外包矩形内
if not isPoiWithinBox(poi, [[0, 0], [180, 90]], tolerance):
return False
polylen = len(simPoly)
sinsc = 0 # 交点个数
for i in range(polylen - 1):
s_poi = simPoly[i]
e_poi = simPoly[i + 1]
if isRayIntersectsSegment(poi, s_poi, e_poi):
sinsc += 1
return True if sinsc % 2 == 1 else False
def isPoiWithinPoly(poi, poly, tolerance=0.0001):
# 点；多边形三维数组；容限
# poly=[[[x1,y1],[x2,y2],……,[xn,yn],[x1,y1]],[[w1,t1],……[wk,tk]]] 三维数组
# 先判断点是否在外包矩形内
if not isPoiWithinBox(poi, [[0, 0], [180, 90]], tolerance):
return False
sinsc = 0 # 交点个数
for epoly in poly: #逐个二维数组进行判断
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
return sinse %2 !=0 #这样更简洁些
#return True if sinsc % 2 == 1 else False
### 比较完备的版本
def pointInPolygon(cin_path,out_path,gfile,t=0):
pindex = [2,3] # wgslng,wgslat 2,3
with open(out_path, 'w', newline='') as cut_file:
fin = open(cin_path, 'r', encoding='gbk')
gfn = open(gfile, 'r', encoding='utf-8')
gjson = json.load(gfn)
if gjson["features"][0]["geometry"]['type']=="MultiPolygon":
plist=gjson["features"][0]["geometry"]['coordinates'] #四维
filewriter = csv.writer(cut_file, delimiter=',')
w = 0
for line in csv.reader(fin, delimiter=','):
if w == 0:
filewriter.writerow(line)
w = 1
continue
point = [float(line[pindex[0]]), float(line[pindex[1]])]
for polygon in plist: #
if isPoiWithinPoly(point, polygon):
filewriter.writerow(line)
break
fin.close()
gfn.close()
elif gjson["features"][0]["geometry"]['type']=="Polygon":
polygon=gjson["features"][0]["geometry"]['coordinates'] #三维
filewriter = csv.writer(cut_file, delimiter=',')
w = 0
for line in csv.reader(fin, delimiter=','):
if w == 0:
filewriter.writerow(line)
w = 1
continue
point = [float(line[pindex[0]]), float(line[pindex[1]])]
if isPoiWithinPoly(point, polygon):
filewriter.writerow(line)
fin.close()
gfn.close()
else:
print('check',gfile)
print('end')
#调用
def baTch():
import os
import glob
wpath="D:/DigitalC_data/coordConvert" #文件路径
sname="D:/DigitalC_data/coordConvertOut"
gpath='D:/cityBoundaryJson/guangzhou_wgs84.json'
for input_file in glob.glob(os.path.join(wpath, '*.csv')):
fname=input_file.split('\\')
pointInPolygon(input_file,os.path.join(sname,fname[1]),gpath)
print(fname[1])
## 应用
### 应用3
'''
对一个csv数据，如果点在多边形内，给某一列(tag)（或者加一列）加上这个多边形的属性（有多个多边形）
'''
def givePolyTag():
pass
### 应用方式1
def pointInPolygon1():
gfile = 'E:/ComputerGraphicsProj/Geojson2kml/J2K_data/深圳poly_bd09.geojson'
cin_path = 'L:/OtherSys/DigitalCityData/深圳特征图层/city_site_poi_sec_shenzhen.csv'
out_path = 'L:/OtherSys/DigitalCityData/深圳特征图层/city_site_poi_sec_shenzhen_out1.csv'
pindex = [4, 5] # wgslng,wgslat 2,3
with open(out_path, 'w', newline='') as cut_file:
fin = open(cin_path, 'r', encoding='gbk')
gfn = open(gfile, 'r', encoding='utf-8')
gjson = json.load(gfn)
polygon = gjson["features"][0]["geometry"]['coordinates'][0]
filewriter = csv.writer(cut_file, delimiter=',')
w = 0
for line in csv.reader(fin, delimiter=','):
if w == 0:
filewriter.writerow(line)
w = 1
continue
point = [float(line[pindex[0]]), float(line[pindex[1]])]
if isPoiWithinSimplePoly(point, polygon):
filewriter.writerow(line)
fin.close()
gfn.close()
print('done')
#pointInPolygon1()
def csvToDArrary(csvpath):#csv文件转二维数组
cdata = []
with open(csvpath, 'r', encoding='gbk') as rfile:
for line in csv.reader(rfile, delimiter=','):
cdata.append(line)
return cdata
### 应用方式2
def pointInPolygon2():
gfile = 'E:/ComputerGraphicsProj/Geojson2kml/J2K_data/深圳poly_bd09.geojson'
cin_path = 'L:/OtherSys/DigitalCityData/深圳特征图层/shenzhen_tAllNotIn.csv'
out_path = 'L:/OtherSys/DigitalCityData/深圳特征图层/shenzhen_tAllNotIn_out2.csv'
pindex = [4, 5] # wgslng,wgslat 2,3
with open(out_path, 'w', newline='') as cut_file:
gfn = open(gfile, 'r', encoding='utf-8')
gjson = json.load(gfn)
polygon = gjson["features"][0]["geometry"]['coordinates'][0]
filewriter = csv.writer(cut_file, delimiter=',')
w = 0
cdata = csvToDArrary(cin_path)
for line in cdata:
if w == 0:
filewriter.writerow(line)
w = 1
continue
point = [float(line[pindex[0]]), float(line[pindex][1])]
if isPoiWithinSimplePoly(point, polygon):
filewriter.writerow(line)
gfn.close()
print('done')
pointInPolygon()
