Geohash与LBS服务

在现在的互联网服务中，基于位置的服务(Location-Based Services, LBS)是一个非常重要的功能。大部分app都需要获取当前的地理位置信息，以提供诸如导航，天气预报，地址选择，附近商家推荐等功能。

本文将对LBS中常用到的一些概念和技术做一些基本介绍。

经纬度坐标

最基本的经纬度就不说了，这是LBS服务的基点。经度以本初子午线为界，分为东经和西经，纬度以赤道为分界线，分为南纬和北纬。以经度和纬度组成的坐标系统就是经纬度系统，用来描述地球表面的位置。

基准经纬度坐标系统也称为地球坐标系，其标准为WGS84（World Geodetic System 1984）。搭载GPS芯片的设备直接获取的经纬度就是WGS84坐标系的真实经纬度。

在国内为了安全需求，国家测绘局制定了GCJ-02（GCJ即"guojia cehui ju"的首字母）坐标系，在原有WGS84坐标基础上增加了一定的偏移量而成，也称为火星坐标系。国家规定所有商用和民用的经纬度都必须使用GCJ-02坐标。

国内各种地图，包括google中国地图使用的均是GCJ-02坐标系。百度地图单独使用百度自己的BD09坐标系。

不同坐标系经纬度转换工具，代码：https://github.com/wandergis/coordtransform。网页工具：https://tool.lu/coordinate。

经纬度小数点与精度关系

经纬度坐标的小数点位数直接影响其表示的精度。通常情况下，纬度和经度的坐标精度可以通过增加小数点后的位数来提高。粗略来说，0位小数表示的精度在百公里级别，1位小数是10公里级别，2位小数是1公里级别，3位小数是100米级别，4位小数是10米级别，以次类推。通常来说，导航或地图服务需要用到5～6位小数，天气预报服务需要用到2或3位小数。

GeoHash算法原理与应用

假设现在有一个场景需要查询周围1公里以内所有的商家，在排除直接调用地图API接口的情况下，该如何实现呢？

假定所有的商家位置经纬度均已录入数据库，根据上文所述，经纬度小数点后第二位每相差0.01，则距离相差约1公里。那我们可以将当前位置的坐标点(lat, lon)作为原点，查询经纬度相差0.01的所有商家信息。查询语句如下：

select * from table where latitude between (lat - 0.01) and (lat + 0.01) and longitude between (lon - 0.01) and (lon + 0.01);

这样将会查出以(lat, lon)为中心，边长约为2公里的正方形内的所有的商家记录。当然严格来说，应该是一个矩形，该矩形在纬度方向上的边长略比2公里大，在经度方向上的边长略比2公里小。

在不严格要求附近1公里，而是允许一定的误差的情况下，到这一步就能达到要求了。但如果严格要求附近1公里，则还要剔除掉其中部分商家。因为附近1公里实际上是以当前位置为圆心，半径为1公里的圆。正方形和圆不相交的区域内商家的距离是超出1公里的，如下图所示：

如果只做简单的查询且数据量有限，只要有联合索引，速度应该还是很快的。但如果要根据距离进行筛选，那响应速度就会随着商家数量增多而下降，因为涉及了大量的计算。而如果数据库表的数据量非常大，接口的并发很高，大量的数据库查询可能会导致严重的性能问题。

有没有更高效的方法呢？实际上是有的，那就是使用Geohash方案。

Geohash编码

GeoHash是一种地址编码方法，能够把二维的空间经纬度数据编码成一个字符串。

如果以经纬度(0, 0)为坐标原点建立一个坐标系，以西、南为负，东、北为正，负为0，正为1，则可以将地球表面分为4块，每块编码如下：

对任一块而言，选择其中心点为坐标原点，可以继续将其分为4块。比如说"11"这一块，代表的经度范围是(0,180), 纬度的范围是(0, 90)。如果选择东经90度，北纬45度为中心点，则又可以分为00, 01, 11, 10四块，与前面的编码11结合起来新的编码如下：

依次对每一个小块递归划分，可以得到更多的区域和对应的编码。

Geohash算法就是基于这种思想，划分的次数更多，得到的区域就更多，每个区域的面积也就更小。将每个区域用0和1进行编码，就能给所有的地理位置进行分区和编码。

下面以经纬度为(116.3895,39.9281)为例，进行实际的编码。

首先看纬度，编码过程如下：

次数	纬度范围	划分区间0	划分区间1	39.9281
1	(-90,90)	(-90,0)	(0,90)	1
2	(0,90)	(0,45)	(45,90)	0
3	(0,45)	(0,22.5)	(22.5,45)	1
4	(22.5,45)	(22.5,33.75)	(33.75,45)	1
5	(33.75,45)	(33.75,39.375)	(39.375,45)	1
6	(39.375,45)	(39.375,42.1875)	(42.1875,45)	0
7	(39.375,42.1875)	(39.375,40.7812)	(40.7812,42.1875)	0
8	(39.375,40.7812)	(39.375,40.0781)	(40.0781,40.7812)	0
9	(39.375,40.0781)	(39.375,39.7265)	(39.7265,40.0781)	1
10	(39.7265,40.0781)	(39.7265,39.9023)	(39.9023,40.0781)	1
11	(39.9023,40.0781)	(39.9023,39.9902)	(39.9902,40.0781)	0
12	(39.9023,39.9902)	(39.9023,39.9462)	(39.9462,39.9902)	0
13	(39.9023,39.9462)	(39.9023,39.9243)	(39.9243,39.9462)	0
14	(39.9023,39.9243)	(39.9023,39.9133)	(39.9133,39.9243)	1
15	(39.9133,39.9243)	(39.9133,39.9188)	(39.9188,39.9243)	1

划分15次时精度范围在500米左右，也方便接下来的base32编码。

则纬度编码为：1 0 1 1 1 0 0 0 1 1 0 0 0 1 1，同理将经度编码为：1 1 0 1 0 0 1 0 1 1 0 0 0 1 0。

将经度在前，纬度在后，两者错位合并得到字符串为：

经度	1		1		0		1		0		0		1		0		1		1		0		0		0		1		0
纬度		1		0		1		1		1		0		0		0		1		1		0		0		0		1		1
合并	1	1	1	0	0	1	1	1	0	1	0	0	1	0	0	0	1	1	1	1	0	0	0	0	0	0	1	1	0	1

接下来按5位一组得到字符串：11100 11101 00100 01111 00000 01101，每组转成10进制数为: 28 29 4 15 0 13。由于每组数字最大为31，因此按照以下规则对其使用base32编码：

0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
0	1	2	3	4	5	6	7	8	9	b	c	d	e	f	g
16	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31
h	j	k	m	n	p	q	r	s	t	bu	v	w	x	y	z

得到编码后的字符串为：wx4g0e，这就是经纬度坐标(116.3895,39.9281)的geohash值。如果对经纬度划分次数更多，则得到的geohash位数也更多，精度更高。

Geohash算法特性

Geohash表示的并不是一个点，而是一个矩形区域。其位数越多，代表的矩形区域就越小，精度越高。

Geohash的编码前缀可以表示更大的区域，这个从上文的编码过程就能分析出来。例如将经纬度(116.3895,39.9281)各划分10次，得到的二进制数合并起来只有20位，经过base32编码之后只有4位，也就是wx4g，此时的纬度范围为：39.9023～40.0781，而wx4g0e表示的纬度范围为：39.9188～39.9243。

下表是geohash长度与实际距离误差：

geohash长度	lat长度	lon长度	lat误差	lon误差	km误差
1	2	3	±23	±23	±2500
2	5	5	±2.8	±5.6	±630
3	7	8	±0.7	±0.7	±78
4	10	10	±0.087	±0.18	±20
5	12	13	±0.022	±0.022	±2.4
6	15	15	±0.0027	±0.0055	±0.61
7	17	18	±0.00068	±0.00068	±0.076
8	20	20	±0.000086	±0.000172	±0.01911
9	22	23	±0.000021	±0.000021	±0.00478
10	25	25	±0.00000268	±0.00000536	±0.0005971

lat, lon长度是指编码后的二进制长度

对上表稍微解释下：

当geohash长度为1时，对应的二进制字符串位数为5，因为经度范围比纬度大，所以经度取3位，纬度取2位。根据上面介绍的编码过程可知，经度3位表示对经度划分了3次，从经度方向来看，将地球表面划成了2^3 = 8块。

同样的，纬度2位表示划分了2次，从纬度方向来看将地球表面划成了2^2 = 4块。两者结合起来，一共将地球表面划分成了32块。每一块的geohash码如下：

此时每块区域表示的纬度范围为180/4=45度，如果将其中心点看作是该块区域的代表位置，则其边界点到中心点的纬度误差为±22.5度。经度范围是360/8=45度，边界点到中心点的经度误差为±22.5度，表中以±23计。最后一列的距离误差实际上是经度方向上的距离误差。上文说过，经度每相差一度，距离大约相差100km，所以±23表示经度方向上距离误差为±2300km，这是一个大概值，实际上更精确一点就是表中所示的2500km。

当geohash的长度增加时，划分的区域更多，对应的经纬度误差更小，距离误差也越小。当geohash长度为10时，距离误差在0.5米左右，已经足够导航级别的应用了。

GeoHash缺陷

上面介绍geohash编码时分别给经纬度编码然后合并，为什么要这样做，原理是什么？

回顾下上文介绍过的将经纬度每块区域按照中心点划分然后编码的方式：

这种类似Z型的空间填充曲线称为Peano曲线，由意大利数学家朱塞佩·皮亚诺 Giuseppe Peano 在1890年首次提出。对Peano曲线来说，大部分情况下编码相似的区域距离相近，这是Peano曲线最大的特点。但也有少数区域编码相邻但距离却相差较远，比如0111与1000，编码是相邻的，距离相差却很大，也就是所谓的矩形编码突变。

事实上除Peano空间填充曲线外，还有很多其他的空间填充曲线。其中效果公认较好是Hilbert空间填充曲线。

解决Geohash的边界问题

Geohash表示的是一个矩形区域，知道当前位置的geohash，则可以查出当前矩形内的所有的POI信息。

POI(Point of Interest)即兴趣点。比如外卖商家，咖啡店，加油站等等，都可以是POI。

但这样还是存在一个问题。看下图：

假设当前位置是X1，如果用X1对应的Geohash值来查询的话，只能查到X2这一个POI, 但是明明附近有一个y1距离更近。

解决这一问题的思路很简单，那就是除了定位点的GeoHash之外，还使用周围8个区域的GeoHash编码进行查询，可以避免这个问题。

逆地理编码

理解了Geohash的基本原理之后，就可以应用到多种实际场景。

比如根据经纬度查询所在位置的行政区，也就是所谓的逆地理编码。

经纬度转geohash

经纬度转geohash很简单，先引入maven库：

<dependency>
    <groupId>ch.hsr</groupId>
    <artifactId>geohash</artifactId>
    <version>{version}</version>
</dependency>

然后直接调用经纬度转geohash方法即可：

public void cvt2GeoHash() {
    double latitude = 31.9749;  // 示例纬度
    double longitude = 118.7594;  // 示例经度
    GeoHash geoHash = GeoHash.withCharacterPrecision(latitude, longitude, 6);
    String geoHashString = geoHash.toBase32();
    System.out.println("Geohash: " + geoHashString);
}
// 输出为：Geohash: wtsmyk

行政区域与Geohash的映射

如果使用长度为2的geohash，则中国地图是可以用29块矩形块来覆盖的，如下所示：

参见geohash地图：https://github.com/hanks-tan/geohashMap

如果使用6位或者7位长度的geohash呢？从geohash长度2精确到长度6时，经度和纬度分别增长2位，对应的二进制都是10位，即经度和纬度各继续划分了10次，最终生成的矩形块数量为：29 * 2¹⁰ * 2¹⁰ = 30408704。 如果使用7位的geohash，那数量就更多了，达到了10亿级别。

地理围栏

中国行政区编码和地理边界坐标数据可以从官方或者地图服务提供商获取。这里以github项目⁵获取的数据为例，南京市雨花台区地理围栏的经纬度示例如下：

118.588269 31.903692,118.592949 31.90721,118.594536 31.909511,118.596855 31.911992,118.59862 31.915976,118.603309 31.921371,118.607655 31.931118,118.608228 31.931264,118.613317 31.937454,118.626557 31.948969,118.641814 31.967451,118.643892 31.969662,118.653251 31.979575,118.655451 31.974806,118.656755 31.971974,118.657955 31.969369,118.658433 31.968311,118.658764 31.967621,118.659347 31.966642,118.660312 31.965526,118.660982 31.964918,118.661782 31.964336,118.662948 31.963841,118.663922 31.963537,118.664487 31.963444,118.666201 31.963133,118.670403 31.962455,118.678563 31.961115,118.685888 31.959907,118.688881 31.95941,118.69303 31.958742,118.696744 31.958115,118.700989 31.957418,118.705494 31.95672,118.711633 31.955861,118.711729 31.956664,118.711642 31.957517,118.711651 31.957669,118.711712 31.957773,118.714312 31.958224,118.714912 31.958435,118.715755 31.958916,118.717886 31.960396,118.718755 31.961211,118.719477 31.961999,118.720894 31.964038,118.721876 31.965721,118.722198 31.966136,118.722824 31.966896,118.724571 31.968584,118.725327 31.969181,118.728004 31.971293,118.730386 31.973011,118.731281 31.973687,118.732028 31.974224,118.73235 31.974527, ... ...

地图服务商一般都会提供限额的免费api调用次数，可用来作为测试使用。比如高德地图，行政区查询接口提供每日5000次免费调用。具体可以注册官方开放平台之后看到。

有了区级行政区域地理围栏后，就可以将geohash与行政区域的映射关系计算出来了。

映射关系计算

Geohash与行政区域映射关系计算过程如下：

1. 遍历所有长度为1的geohash块，判断每一块是否与国家地理围栏相交，如果是则将这一块继续划分为32块，得到长度为2的geohash块，以此类推，直到geohash长度为5，如果该geohash块与地理围栏相交，则将其加入到geohashList列表中。
2. 遍历geohashList，判断每一块是否与某个行政区地理围栏相交，如果是，则取出所有与其相交的行政区地理围栏：
   2.1 遍历与该geohash块相交的地理围栏，判断该geohash块是否在该地理围栏内，如果是，则将geohash与行政区映射关系放到一个blockingQueue中；如果不是，则继续将该geohash块划分成32个geohash子块，继续执行此步骤，直到geohash块与行政区映射关系全部放入blockingQueue中，或者geohash长度为7，此时强制将geohash与行政区的映射关系放入blockingQueue中。
3. 从blockingQueue中取出gehash与行政区的映射关系将其存到数据库或者redis中。

下面对这个过程进行具体说明。

根据前文介绍可知，如果geohash长度为1，则地球表面可划分为32块，其中中国陆地地图与t、v、y、w四块相交。每一块继续划分为32块，得到长度为2的geohash，则有29块相交。当geohash长度为5时，中国行政区域占据多少geohash块呢？

在Java中可以通过Geometry库的Polygon类的intersects方法来计算出来。

该方法可以判断两个多边形是否相交，只要有至少一个相同的点，则返回true。因此，此方法判断的多边形相交包含两种情况，一种是普通意义上的相交，另一种是包含。

下面是伪代码实现：

for each geohash1 in allGeohashLength1 {
    if (polygonIntersects(geohash1, contryPolygon)) {
        geohash2List = subdivide(geohash1)
        for each geohash2 in geohash2List {
            if (polygonIntersects(geohash2, contryPolygon)) {
                geohash3List = subdivide(geohash2)
                for each geohash3 in geohash3List {
                    if (polygonIntersects(geohash3, contryPolygon)) {
                        geohash4List = subdivide(geohash3)
                        for each geohash4 in geohash4List {
                            if (polygonIntersects(geohash4, contryPolygon)) {
                                geohash5List = subdivide(geohash4)
                                for each geohash5 in geohash5List {
                                    if (polygonIntersects(geohash5, contryPolygon)) {
                                        geohashList.add(geohash5)
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

为了方便理解，这里用5个嵌套的for循环来判断geohash块与行政围栏是否相交。事实上可以用stream操作来简化代码。

intersects方法判断的是两个多边形的相交关系，因此geohash需要转换成Polygon类型。一个geohash块是一个矩形，有四个点，将这四个点连起来就是一个Polygon。代码如下：

private boolean polygonIntersects(GeoHash geoHash, Polygon polygon) {
    GeometryFactory geometryFactory = new GeometryFactory();
    // Create JTS Polygon from GeoHash
    Coordinate[] coordinates = new Coordinate[5];
    coordinates[0] = new Coordinate(geoHash.westLng, geoHash.southLat);
    coordinates[1] = new Coordinate(geoHash.eastLng, geoHash.southLat);
    coordinates[2] = new Coordinate(geoHash.eastLng, geoHash.northLat);
    coordinates[3] = new Coordinate(geoHash.westLng, geoHash.northLat);
    // Closing the ring
    coordinates[4] = coordinates[0];
    Polygon geoHashPolygon = geometryFactory.createPolygon(coordinates);
    return polygon.intersects(geoHashPolygon);
}

另外，判断geoHash多边形是否与国家地理围栏多边形相交，等价于判断geoHash多边形是否与任意一个区(县)级行政区域多边形相交。如果无法获取国家地理围栏多边形，用区(县)级地理围栏多边形也是可以的。

最终得到一个包含5位长度的geohash列表，该列表内每个geohash块均与某个区(县)级行政区域相交。

接下来遍历geohashList，取出与当前geohash块相交的行政区域多边形，遍历这些行政区域，判断当前geohash是否在这个行政区域内，如果在，就得到其映射关系。否则将该geohash继续划分32块，遍历这32块geohash来做同样的操作，直到geohash长度为7为止。伪代码实现如下：

function intersectPolygons(geohash, polygons) {
	for each geohash in geohashList {
    	for each polygon in polygons {
        	if (intersects(geohash, polygon)) {
            	intersectPolygons.add(polygon);
            }
        }
    }
}

// 定义函数处理 geohash 块与行政区的映射关系
function processGeohash(geohash, intersectPolygons) {
	for each polygon in intersectPolygons {
		if (isWithin(geohash, polygon)) {
        	blockingQueue.add(new Mapping(geohash, district))
        } else if(geohash.length >= 7 ) {
        	blockingQueue.add(new Mapping(geohash, district))
        } else {
        	subGeohashes = subdivide(geohash)
        	for each subGeohash in subGeohashes {
        		processGeohash(subGeohash, intersectPolygons(subGeohash, polygons))
        	}
        }
	}
}

判断一个多边形在另一个多边形内，可以用Polygon的contains方法。

为什么处理geohash与行政区的映射关系时有三个条件分支呢？这是因为geohash块与区级行政区域的位置关系有两种情况，一种是包含关系，如下所示：

此时geohash与行政区域可以直接映射，例如: wtsgk -> {“province”: “xx省”, “city”: “xx市”, “admin”: “xx区”}。

还有一种就是geohash与多个行政区域相交，如下所示：

此时显然是不能直接将geohash映射到区域。因此需要将该geohash继续细分，从5位长度到6位，得到新的geohash块如下：

遍历新的geohash块，对每一块执行相同的操作逻辑，直到geohash长度到7为止。此时即使长度为7的geohash块仍然是与行政区域相交，而不是包含于某一个行政区，也将geohash与行政区域的映射关系保存下来。

为什么不继续划分了呢？根据上文geohash的长度和精度之间的关系可知，当geohash长度为7时，误差为±76m，这个误差对于判断区级行政区域来说是可以接受的。如果要更精确一点，还可以继续划分到8位长度。但与之对应的，geohash的数量就多了很多，需要更多的计算和存储资源。

行政区域的调整

通过上面的介绍可以看出，将geohash映射到区(县)级行政区域的计算量是比较大的，但好在通常情况下这种计算只要做一次就行了。如果说考虑到国家对区(县)级行政区域会有调整，则需要定期重新计算geohash与行政区域的映射关系。至于这个定期是什么频率，就取决于你的业务能忍受多长时间的数据延迟了。

重新计算映射关系时，如何更新数据呢？

考虑如下场景：

• 原来一个5位或6位长度的geohash块在某个行政区域内(包含的)，调整区域边界后，该geohash与行政区域是相交(非包含)的。如下图所示：

• 原来一个7位长度的geohash与多个行政区域相交(非包含)的，调整区域边界后，该geohash被包含到某个区域了。如下所示：

对于第一种情况，原geohash与行政区域的映射关系变成无效的了，它将被更高位数的geohash来代替。如wtx1b这个geohash块被划分成了6位的wtx1bj，wtx1bk，…。而wtx1d2这个geohash块被划分成了7位的wtx1d2b，…。这种情况下，新增了一些6位和7位的geohash，原对应的5位和6位geohash则是无效的，需要被移除。

对于第二种情况，原7位长度的geohash与行政区域的映射关系变成无效的了，它将被更低位数的geohash来代替。如wtx1dnf向上升级到了6位wtx1dn。如果此时这个升级后的geohash与另外的31个6位geohash都在同一个行政区域内，则可以继续升级为wtx1d。这种情况下，新增了6位或5位的geohash，而原对应的7位geohash需要被移除。

实际上这两种情况可以合并判断：当前新增的geohash长度为len，则其6或者5位的父geohash需要被移除，其6或者7位的子geohash需要被移除。伪代码如下：

for (String geohash : newAddedGeohashes) {
    int len = geohash.length();
    // parent geohashs
    for (int i = minPrecision; i < len; i++) {
        geoHashesToRemove.add(geohash.substring(0, i));
    }
    // children geohashes
    parent = geohash;
    for (int i = len + 1; i <= maxPrecision; i++) {
        List<Geohash> children = subdivide(parent);
        List<String> geocodes = children.stream().map(Geohash::getGeocode).toList();
        geoHashesToRemove.addAll(geocodes);
        parent = children;
    }
}

最后将geohash与行政区域的映射关系保存到数据库或者redis中，就能根据经纬度来查询对应的行政区域了，当然在边界上是有误差的，精度为±76m。

Redis-GeoHash

Redis从3.2版本开始增加了对Geohash的支持。输入一个经纬度，可以得到其对应的geohash值，还能计算不同位置的距离。

命令行操作

添加

127.0.0.1:6379> geoadd loc 116.3895 39.9281 loc1 116.3945 39.9123 loc2
(integer) 2

查询geohash

127.0.0.1:6379>  geohash loc loc1
"wx4g0s8npc0"
127.0.0.1:6379>  geohash loc loc1 loc2
"wx4g0s8npc0"
"wx4g09gfh20"

redis的geohash长度默认为11，不可修改

查询经纬度

127.0.0.1:6379> geopos loc loc1 loc2
"116.38950079679489136"
"39.92810001044335166"
"116.39450043439865112"
"39.91230109345719512"

计算距离

127.0.0.1:6379> geodist loc loc1 loc2 km
"1.8083"
127.0.0.1:6379> geodist loc loc1 loc2 m
"1808.2692"

根据member查询附近指定范围内的位置点

127.0.0.1:6379> georadiusbymember loc loc1 1 km withdist count 2 desc 
"loc1"
"0.0000"
127.0.0.1:6379> georadiusbymember loc loc1 2 km withdist count 3 asc 

"loc1"
"0.0000"
"loc2"
"1.8083"

radius表示距离是半径。注意，此命令获取的结果包含本身在内

根据经纬度查询附近指定范围内的位置点

127.0.0.1:6379> georadius loc 116.38 39.92 2 km withcoord withdist count 3 asc
"loc1"
"1.2118"
"116.38950079679489136"
"39.92810001044335166"

"loc2"
"1.5045"
"116.39450043439865112"
"39.91230109345719512"

点射法判断点是否在多边形内

点射法理论是：从当前点向右水平发出一条射线，判断射线与多边形的相交点的个数，如果为奇数则在多边形内，否则在多边形外。

代码如下非常简单如下所示：

public boolean isPointInsidePolygon(Point[] polygon, double lng, double lat) {
    boolean inside = false;
    int n = polygon.length;
    // j 指向上一个顶点（循环初始时取最后一个顶点）
    for (int i = 0, j = n - 1; i < n; j = i++) {
        double xi = polygon[i].lng;
        double yi = polygon[i].lat;
        double xj = polygon[j].lng;
        double yj = polygon[j].lat;
        // 判断点的纬度是否在这条边的两端纬度之间
        // 并计算射线与边的交点的经度，如果点的经度小于交点经度，则视为有一次交点
        // 通过直线边的参数方程(x,y)=(lngI,latI)+t×(lngJ−lngI,latJ−latI) 可计算交点经度
        boolean intersect = ((yi > lat) != (yj > lat)) &&
                            (lng < (xj - xi) * (lat - yi) / (yj - yi) + xi);
        if (intersect) {
            inside = !inside;  // 每遇到一次交点，inside 状态取反
        }
    }
    return inside;
}

至于多边形的坐标点，可以从地图服务商获得。这里就不展开了。

参考资料

[1]. https://www.cnblogs.com/aiweixiao/p/6188081.html
[2]. https://www.cnblogs.com/feiquan/p/11380461.html
[3]. https://en.wikipedia.org/wiki/Geohash
[4]. https://blog.csdn.net/qq_20128967/article/details/108340233
[5]. https://github.com/xiangyuecn/AreaCity-JsSpider-StatsGov.