如何通过经纬度反查城市行政区域

发表于 分类于 阅读次数：
本文主要介绍如何通过经纬度反查城市行政区域，主要分为两个步骤：1. 采集省市区边界数据，导入到 MongoDB、PostgreSQL 等支持地理空间索引的数据库中；2. 对于给定的经纬度，可以通过边界数据查询所在的行政区域。边界数据来源于开源项目，需要经过数据处理和清洗。另外，本文还介绍了行政区划代码、坐标系转换，并给出了源码实现和一键化部署的仓库地址。

前言

业务需求中需要根据经纬度反查具体的城市行政区域，需要精确到区县级别，调用频率正常情况下是 100 次/秒，但是峰值情况下可能会达到 100K 次/秒。如果是业余项目，可以使用 根据经纬度查城市 这个工具网站获取，但是频率和次数都有限制，不适合生产环境。

一开始是打算使用地图服务商提供的开放坐标拾取 API，但是在细致调研后放弃了这个方案，主要原因有：

  1. 地图服务商提供的免费 API 有调用次数和频率限制，如果使用收费版本，就带来了额外的成本，我们的业务场景并不需要地图服务商提供的其他功能，感觉有点浪费。
  2. 我们的业务场景要求精度其实不高，只需要到区县级别，而且对于数据的有效性和实时性并没有特别高的要求。
  3. 数据保密性，地图服务商提供的 API 会将经纬度信息传递给服务商，这可能涉及到用户隐私问题，这点也是我们最为担心的。

在经过相关资料的调研后，笔者发现根据经纬度反查城市行政区域的方法其实并不复杂，主要分为两个步骤：

  1. 首先采集省市区边界数据，导入到 MongoDB、PostgreSQL 等支持地理空间索引的数据库中。
  2. 对于给定的经纬度，可以通过边界数据查询所在的行政区域，MongoDB 中对应的是 geoIntersects 操作符，PostgreSQL 中对应的是 ST_Contains 函数。

最大困难其实在于第一步，如何拿到最新的省市区边界数据呢？我在 github 上找到了一个项目 GitHub - xiangyuecn/AreaCity-JsSpider-StatsGov，里面免费提供了最新的国家省市区边界数据（四级行政区域需要收费），那接下来的步骤就很简单了，将数据导入到 MongoDB 数据库中，然后创建地理空间索引，就可以实现经纬度反查城市行政区域了。

实现仓库地址：convexwf/geo-inverse-query: This is a simple Python script that queries the region of a given latitude and longitude.

开箱即用，只需要 docker 和 docker-compose 环境，一键启动服务，即可实现经纬度反查城市行政区域，具体细节可以看项目的 README 文档，原理和实现流程如下。

省市区边界数据采集

源数据来自 github 上的 GitHub - xiangyuecn/AreaCity-JsSpider-StatsGov，当前最新版本：2022.230704.230813版，更新于 2023-08-13，整合了 统计局2022-12-29、民政部2023-04-23、腾讯地图行政区划2023-07-04、高德地图行政区划采集当天 数据。

数据下载链接：最新2023年省市区县乡镇街道行政区划数据 - AreaCity-JsSpider-StatsGov。

数据字段包括 id、pid、deep(级别，级别为 1 表示省级，级别为 3 表示区县级)、name(当前行政区域名称)、ext_path(当前行政区域全称)、geo(行政区域中心，一维数组，经纬度)、polygon(行政区域边界)。

polygon 行政区域边界可能包含多个多边形，比如北京市东城区就是由两个多边形组成的，不同的多边形用 ; 符号间隔。多边形由多个点组成，不同点用 , 符号间隔，点的表示形式为 {longitude} {latitude}。

样例数据如下（polygon 精简化）：

1
2
3
4
5
6
7
8
9
{
  "id": "110101",
  "pid": "110100",
  "deep": 3,
  "name": "东城区",
  "ext_path": "北京市 北京市 东城区",
  "geo": [116.416357, 39.928353],
  "polygon": "116.416357 39.928353,116.416357 39.928353;116.416357 39.928353,116.416357 39.928353"
}

需要从原始数据中抽取出需要的字段，然后将数据导入到 MongoDB 数据库中。在此过程需要注意几点：

  1. 由于 MongoDB 中的地理空间索引只支持 GeoJSON 格式，因此需要将原始数据中的经纬度坐标转换为 GeoJSON 格式。原始数据中的 polygon 对应 MongoDB 中的 MultiPolygon 类型，具体类型介绍见 GeoJSON Objects MultiPolygon— MongoDB Manual。
  2. MongoDB 要求 MultiPolygon 对象至少包含一个多边形，并且每个多边形必须至少有 4 个点（起始点和结束点相同）。由于原始数据中的 polygon 边界数据中并不是闭合形状，因此需要手动将边界数据首尾相连，形成闭合形状后再导入到 MongoDB 中。

导入 MongoDB 数据库后，需要针对边界数据创建地理空间索引，以便于后续查询。

1
db.city.createIndex({ "city_boundary": "2dsphere" })

查询：根据经纬度查询所在的行政区域。

1
2
3
4
5
6
7
8
9
10
db.city.find({
  city_boundary: {
    $geoIntersects: {
      $geometry: {
        type: "Point",
        coordinates: [113.832874, 22.647417]
      }
    }
  }
});

行政区划代码

上面的数据字段里还包含了 id 和 pid，其实就是行政区划代码。当然也可以自己给不同的行政区域分配不同的城市 ID，但是现有的行政区划代码已经包含了省市区县的信息，因此可以直接使用行政区划代码来查询。

数据来源：2022年中华人民共和国县以上行政区划代码，共有 3208 条记录，详细记录了全国三级行政区域的名称和行政区划代码，经过简单处理后得到数据如下：

1
2
3
4
5
6
7
8
9
10
11
单位名称,行政区划代码
北京市 北京市 北京市,110000
北京市 北京市 东城区,110101
北京市 北京市 西城区,110102
北京市 北京市 朝阳区,110105
北京市 北京市 丰台区,110106
北京市 北京市 石景山区,110107
北京市 北京市 海淀区,110108
北京市 北京市 门头沟区,110109
北京市 北京市 房山区,110111
......

观察可知，行政区划代码的前两位表示省级行政区划代码，中间两位表示地级行政区划代码，后两位表示县区级行政区划代码，其中直辖市（北京、上海等）是没有下辖地级市的。在涉及到城市区域相关的查询时，可以通过行政 ID 前缀匹配的方式快速定位到下辖区域的集合，比如说需要查找广东省的所有地级市，可以通过前两位 44 进行前缀匹配。

坐标系转换

在实际应用中，经纬度坐标系一般分为两类：第一类是 WGS84 (WorldGeodetic System 1984) 坐标系，是由美国国防部制定的一种坐标系，用于 GPS 系统，谷歌地球、谷歌地图等地图服务商使用的坐标系；第二类是 GCJ02 (国测局坐标系)，是中国国家测绘局制定的一种坐标系，用于国内地图服务，高德地图、腾讯地图等地图服务商使用的坐标系。

另外还有一种坐标系是 BD09 (百度坐标系)，是百度地图使用的坐标系，是在 GCJ02 坐标系的基础上进行了加密处理，此处不做详细介绍。

在实际应用中，经常会遇到 WGS84 和 GCJ02 之间的坐标转换，以上面的数据为例，省市区边界数据采集的坐标系是 GCJ02，如果我们需要查询的是 WGS84 坐标系下的经纬度，就需要进行坐标转换。

坐标系转换的方法有很多开源库实现，本文采用的是 GitHub - wandergis/coordTransform_py 的 Python 版本实现。

Reference