目前实验中有个需求，即给定一个经纬度，判断其是否处于某个城市内部。本来是想使用Google Earth的接口，然而谷歌一直不通过我的开发者账号申请，遂自己写了一个程序来实现该功能。
首先一个城市的轮廓必然是多边形的，那么得先获取该多边形的每个点的坐标，即城市多边形边界的表示。
怎么找城市边界的坐标呢？如下图所示，下图来自于：Getting polygon boundaries of City in JSON from Google Maps API?
打开图中链接：https://nominatim.openstreetmap.org/
然后按照图中指示的操作即可获取所需地区的边界，有txt和json两种格式的数据。
下面以纽约为例。
获取到的简化版的地图边界信息如下图所示（当然里面也有原始的边界信息，只是转折点比较多而已）：
数据的表示形式如下图所示：
这里我将无关信息去掉，只保留坐标信息，下图为另一个区域的一个简单的表示。可以看出第一个和最后一个节点是同一个，即这里的点是沿顺时针顺序进行的。
数据获取与预处理完成，该开始写核心算法了。
下面的核心问题就是：判断一个点是否在多边形区域内部
随便百度一下，就会知道有一个行之有效的方法叫“射线法”，但是关于射线法的使用详细情况却往往未做说明，这其中还是有很多坑的。这篇文章将射线法中所要遇到的问题进行了全面的分析与解释，也提出了解决办法，因此我不多赘述。判断一个坐标点是否在不规则多边形内部的算法
但是该文章并未给出代码实现，因此本文在最后基于上述所说的应用场景给出一个python实现版本的例子。
# region1是一个list，若以上图为例，则其格式为：[[-74.035, 40.695],[-74.036, 40.694]...]
# 该函数就是判断一个给定点是否在region1这个list所围成的多边形内
# 需要注意：region1中第一个元素代表经度，第二个代表纬度，而函数实际输入中第一个是纬度，第二个是经度
def isInRegion_1(lat, lon):
count = 0
lat = float(lat)
lon = float(lon)
for i in range(len(region_1)):
# 如果是最后一个元素，那么必然是和第一个一样的，就啥也不干
if i+1 == len(region_1):
break
# 如果不是最后一个元素，那么需要和后一个元素一起判断给定点是否在区域内
la_1, lo_1 = region_1[i]
la_2, lo_2 = region_1[i+1]
la_1, lo_1, la_2, lo_2 = float(la_1), float(lo_1), float(la_2), float(lo_2)
# 以纬度确定位置，沿纬度向右作射线，看交点个数
if lat < min(la_1, la_2):
continue
if lat > max(la_1, la_2):
continue
# 如果和某一个共点那么直接返回true
if (lat, lon) == (la_1, lo_1) or (lat, lon) == (la_2, lo_2):
return True
# 如果和两点共线
if lat == la_1 == la_2:
if lon < min(lo_1, lo_2) or lon > max(lo_1, lo_2):
return False
else:
return True
# 接下来只需考虑射线穿越的情况，该情况下的特殊情况是射线穿越顶点
# 求交点的精度
cross_lon = (lat - la_1) * (lo_2 - lo_1) / (la_2 - la_1) + lo_1
# 因为需要考虑特殊情况的射线穿越顶点，因此只判断是不是穿越左边这个顶点
if cross_lon == la_1:
count += 1
# 其他情况
elif cross_lon > lon:
count += 1
if count%2 == 0:
return False
return True