判断点是否在多边形内的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保存的边界数据也不难改。文中代码同步于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')

本文分享自微信公众号 - 蛰虫始航(lyns_sailing)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-09-03

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Android机动车

Android中JSON库性能比较

JSON不管是在Web开发还是服务器开发中是相当常见的数据传输格式,一般情况我们对于JSON解析构造的性能并不需要过于关心,除非是在性能要求比较高的系统。 目前...

22020
来自专栏林德熙的博客

VisualStudio 给任意字符串给定正则等格式

在写正则或 json 等字符串的时候,期望有智能提示和语法纠错,在 VisualStudio 可以在字符串上面添加一句注释表示这个字符串的功能是什么,然后 Vi...

10110
来自专栏Java技术栈

什么是 Java 对象深拷贝?面试必问!

浅拷贝只是拷贝了源对象的地址,所以源对象的值发生变化时,拷贝对象的值也会发生变化。

10430
来自专栏一猿小讲

7行代码搞定WEB服务

作为一个 Java 程序猿,写代码久了,各种技术也就都尝试了一个遍。先从 SSH1(Spring、Struts1、Hibernate)摸爬滚打转变到 SSH2(...

10220
来自专栏前端一会

前端学习与工作书籍推荐

13050
来自专栏Web行业观察

你所不知道的ndJSON:序列化与管道流

一直以为对JSON所有的语法都了如指掌,毕竟json的标准用一只手都数的过来,直到我发现了一个叫ndJSON的标准,简单说,以下2种语法都是合法的:

22140
来自专栏前端重点笔记

深拷贝和浅拷贝原来是这样?

为了让读者更好的理解深浅拷贝,在讲深浅拷贝之前要引入基本数据类型 , 引用数据类型 和 数据储存(栈和堆)这几个概念,如果已经理解,可直接跳过这一part。

10330
来自专栏IT码农

java之@Controller和@RestController以及@GetMapping和@PostMapping接收参数的格式使用

一、1.使用@Controller 注解,在对应的方法上,视图解析器可以解析return 的jsp,html页面,并且跳转到相应页面

35020
来自专栏IT码农

java之@RequestBody的使用

@RequestBody主要用来接收前端传递给后端的json字符串中的数据的(请求体中的数据的);GET方式无请求体,所以使用@RequestB...

36810
来自专栏程序员的成长之路

JSP还有必要学吗?这篇文章告诉你

前后端分离已成为互联网项目开发的业界标准使用方式,通过nginx+tomcat的方式(也可以中间加一个nodejs)有效的进行解耦,并且前后端分离会为以后的大型...

20140

扫码关注云+社区

领取腾讯云代金券

年度创作总结 领取年终奖励