#!/usr/bin/env python # -*- coding: utf-8 -*- # @Date : 2018-10-07 15:49:37 # @Author : Sheldon ([email protected]) # @Blog : 谢耳朵的派森笔记 # @Link : https://www.cnblogs.com/shld/ # @Version : 0.0.1 def isinpolygon(point,vertex_lst:list, contain_boundary=True): #检测点是否位于区域外接矩形内 lngaxis, lataxis = zip(*vertex_lst) minlng, maxlng = min(lngaxis),max(lngaxis) minlat, maxlat = min(lataxis),max(lataxis) lng, lat = point if contain_boundary: isin = (minlng<=lng<=maxlng) & (minlat<=lat<=maxlat) else: isin = (minlng<lng<maxlng) & (minlat<lat<maxlat) return isin def isintersect(poi,spoi,epoi): #输入:判断点,边起点,边终点,都是[lng,lat]格式数组 #射线为向东的纬线 #可能存在的bug,当区域横跨本初子午线或180度经线的时候可能有问题 lng, lat = poi slng, slat = spoi elng, elat = epoi if poi == spoi: #print("在顶点上") return None if slat==elat: #排除与射线平行、重合,线段首尾端点重合的情况 return False if slat>lat and elat>lat: #线段在射线上边 return False if slat<lat and elat<lat: #线段在射线下边 return False if slat==lat and elat>lat: #交点为下端点,对应spoint return False if elat==lat and slat>lat: #交点为下端点,对应epoint return False if slng<lng and elat<lat: #线段在射线左边 return False #求交点 xseg=elng-(elng-slng)*(elat-lat)/(elat-slat) if xseg == lng: #print("点在多边形的边上") return None if xseg<lng: #交点在射线起点的左侧 return False return True #排除上述情况之后 def isin_multipolygon(poi,vertex_lst, contain_boundary=True): # 判断是否在外包矩形内,如果不在,直接返回false if not isinpolygon(poi, vertex_lst, contain_boundary): return False sinsc = 0 for spoi, epoi in zip(vertex_lst[:-1],vertex_lst[1::]): intersect = isintersect(poi, spoi, epoi) if intersect is None: return (False, True)[contain_boundary] elif intersect: sinsc+=1 return sinsc%2==1 if __name__ == '__main__': vertex_raw = ['113.77253623662,22.734971291583', '113.77255514955,22.73509716357', '113.77655989881,22.736350557773', '113.77676822305,22.73640965307', '113.77736723203,22.736534992321', '113.77940615018,22.736878669244', '113.77947030985,22.736860953444', '113.77949943202,22.736797638175', '113.77996334483,22.734671733626', '113.78003353483,22.734638924801', '113.78109183843,22.734901179433', '113.78111784887,22.734964846247', '113.78115384936,22.735067381603', '113.78114374702,22.735157490465', '113.78107244941,22.73530034639', '113.78095676023,22.735766684331', '113.78094658962,22.735932778191', '113.78103664924,22.736114636979', '113.78116482396,22.736273022099', '113.78134419131,22.736359769194', '113.78211401039,22.73651290222', '113.78222325602,22.736555467145', '113.7823555342,22.73663071491', '113.78248580377,22.736711977153', '113.78259294826,22.736865523263', '113.7826400563,22.736883891204', '113.78273032433,22.73685368898', '113.78357214618,22.735114513862', '113.78288857701,22.734768949082', '113.78456713092,22.731111070805', '113.78499799148,22.730246917863', '113.78507426121,22.730178805366', '113.7855499059,22.72982978801', '113.7855439672,22.729710894129', '113.78547284386,22.729577977108', '113.78485039809,22.72898526075', '113.7844013172,22.728642818037', '113.78407147738,22.728484550981', '113.78369248319,22.728371904397', '113.7832943991,22.728348425846', '113.78279101654,22.728381265282', '113.78223548106,22.728490624284', '113.78120871157,22.728746702049', '113.78070843445,22.728835826604', '113.78028539185,22.728892874896', '113.77960479337,22.728941739674', '113.77917585519,22.728917533458', '113.77850550189,22.728784767411', '113.77800554816,22.728664945674', '113.77713301773,22.728385531131', '113.77575487604,22.727878867397', '113.77366621063,22.727101561255', '113.77346692784,22.72714784993', '113.77333268596,22.727276687378', '113.77266852558,22.733599777451', '113.77264730132,22.73384585746', '113.77253623662,22.734971291583'] vertex_lst = [[float(strnum) for strnum in poistr.split(",")] for poistr in vertex_raw] poi = [113.7726,22.735] print(isin_multipolygon(poi,vertex_lst, contain_boundary=True))
判断点是否在区域的python实现(射线法)
猜你喜欢
转载自www.cnblogs.com/shld/p/9758303.html
今日推荐
周排行