经纬度转城市反地理查询系统搭建

想知道某个经纬度属于哪个城市，通常可以通过地图 API 的接口实现。但是地图服务商的 API 通常会有配额限制。问题来了，能否搭建自己的经纬度反查系统呢？

GADM

GADM 主页：https://gadm.org/

GADM，全称 Database of Global Administrative Areas，是一个高精度的全球行政区划数据库。其包含了全球所有国家和地区的国界、省界、市界、区界等多个级别的行政区划边界数据。

警告：GADM 提供的中国国界数据不符合中国的领土主张，省界、市界、区界等数据也不一定是最新的版本。此方案仅供系统测试，此部分使用不建议使用。

GADM 提供了两种下载方式：

• 下载全球所有国家和地区的所有数据https://gadm.org/download_world.html
• 按国家下载https://gadm.org/download_country_v3.html

由于全球数据量巨大，建议根据需要按照国家下载数据。需要说明的是，GADM 中对 country 的定义为“any entity withan ISO country code”。因而如果想要下载完整的中国数据，实际上需要下载 China、Hong Kong、Macao 和 Taiwan 四个数据。

对于每个数据，GADM 提供了 5 种不同的格式：

• Geopackage：文件后缀是 gpkg，矢量格式，官网对其评价是非常好的空间地理存储格式的文件。可以被 GDAL/OGR、ArcGIS、QGIS 等软件读取。但是在国内用的很少，网上介绍处理文件的方法，大多数是使用 shp 文件。
• Shapefile：可直接用于 ArcGIS 等软件，把官网的文件下载后解压后，得到 (.shp, .shx, .dbf, .prj) 组成的文件。
• KMZ：可直接在 Google Earth 中打开
• R (sp)：可直接用于 R 语言绘图
• R (sf)：可直接用于 R 语言绘图

PostGIS

PostGIS 是 PostgreSQL 关系数据库的空间操作扩展。它为 PostgreSQL 提供了存储、查询和修改空间关系的能力。空间数据库像存储和操作数据库中其他任何对象一样去存储和操作空间对象。

PostGIS 特性与功能：

• PostGIS 支持所有的空间数据类型。这些类型包括：点（POINT）、线（LINESTRING）、多边形（POLYGON）、多点（MULTIPOINT）、多线（MULTILINESTRING）、多多边形（MULTIPOLYGON）和集合对象集（GEOMETRYCOLLECTION）等。
• PostGIS 支持所有的对象表达方法。比如 WKT 和 WKB。
• PostGIS 支持所有的数据存取和构造方法。如 GeomFromText()、AsBinary()，以及 GeometryN() 等。
• PostGIS 提供简单的空间分析函数。如 Area 和 Length，同时也提供其他一些具有复杂分析功能的函数。比如 Distance。
• PostGIS 提供了对于元数据的支持。如 GEOMETRY_COLUMNS 和 SPATIAL_REF_SYS，同时，PostGIS 也提供了相应的支持函数，如 AddGeometryColumn 和 DropGeometryColumn。
• PostGIS 提供了一系列的二元谓词（如 Contains、Within、Overlaps 和 Touches）用于检测空间对象之间的空间关系，同时返回布尔值来表征对象之间符合这个关系。
• PostGIS 提供了空间操作符（如 Union 和 Difference）用于空间数据操作。比如，Union 操作符融合多边形之间的边界。两个交迭的多边形通过 Union 运算就会形成一个新的多边形，这个新的多边形的边界为两个多边形中最大边界。

相关链接：

• https://www.postgresql.org/
• https://postgis.net/

使用 GADM 搭建反地理系统

1、下载 GADM 数据，这里下载的是 Shapefile 文件 gadm36_shp.zip。解压后获得如下文件：

• gadm36.cpg
• gadm36.dbf
• gadm36.prj
• gadm36.shp
• gadm36.shx

2、创建数据库

CREATE DATABASE mygis
WITH
OWNER = postgres
ENCODING = 'UTF8'
LC_COLLATE = 'Chinese(Simplified)_China.936'
LC_CTYPE = 'Chinese(Simplified)_China.936'
TABLESPACE = pg_default
CONNECTION LIMIT = -1;
CREATE extension postgis;

3、将 shapefile 格式数据导入到 PostgreSQL 中

PostgreSQL 中自带了 shp2pgsql.exe 可帮忙完成此操作。shp2pgsql.exe 的使用方法如下：

USAGE: shp2pgsql [OPTIONS] shapefile [schema.]table
OPTIONS:
-s <srid>
 设置 srid，缺省为 -1
(-d|a|c|p) 互斥选项:
-d 重新建立表，并插入数据。
-a 在同一个表中增加数据
-c 建立新表，并插入数据。(缺省)
-p 只创建表
-g <geocolumn> 指定要创建的表的空间字段名称(在追加数据时有用)
-D 使用 dump 方式，缺省是生成 sql
-G Use geography type (requires lon/lat data).
-k 保持 PostgreSQL 标识符方式
-i 使用 int4 类型 dbf 文件里的 integer 类型
-I 在空间字段上建立索引
-S Generate simple geometries instead of MULTI geometries.
-W <encoding> shape 文件属性列的字符格式。缺省是 ascII
-N <policy> 指定 geometries 为空时的操作(insert, skip, abort)
-n 只导入 dbf 文件
-? 显示帮助

参考链接：https://www.bostongis.com/pgsql2shp_shp2pgsql_quickguide.bqg

具体命令：

shp2pgsql.exe -s 4326 -I .\gadm36.shp public.gadm36 > gadm36.sql
psql.exe -h localhost -p 5432 -U postgres -d mygis -f .\gadm36.sql

其中 SRID 4326 可通过http://prj2epsg.org/ 查询获得。

执行后报如下错误：

psql:./gadm36.sql:15:错误: 编码 "GBK" 的字符 0x0x820x27 在编码 "UTF8" 没有相对应值
psql:./gadm36.sql:21:错误: 语法错误 在 "," 或附近的
第 1 行, 00FCF0C41D152, CF6FD0FF, 00C301B, B6DBE3D0DF000080CD46A, 0CC, 3,...

尝试进行修改：

shp2pgsql.exe -s 4326 -I -W GBK .\gadm36.shp public.gadm36 > gadm36.sql

报错内容变为：

Unable to convert data value to UTF-8 (iconv reports "Illegal byte sequence"). Current encoding is "GBK". Try "LATIN1" (Western European), or one of the values described at http://www.postgresql.org/docs/current/static/multibyte.html.

修改为：

shp2pgsql.exe -s 4326 -I -W LATIN1 .\gadm36.shp public.gadm36 > gadm36.sql

报错内容消失，数据可正常导入.

尝试通过经纬度获取所在城市：

SELECT * FROM public.gadm36 WHERE ST_Within(ST_SetSRID(ST_MakePoint(121.55223, 38.86758), 4326), geom);

SQL 能正常执行，但数据库中的一些字段由于先前选择了 LATIN1 编码而导致了乱码：

尝试了很多方法，无法解决上面的编码问题，PostgreSQL 在建立数据库时也无法选择 GBK 和 GB13080。最终方案：

更换方式，使用 ogr2ogr 导入 Geopackage 数据：

ogr2ogr.exe -f PostgreSQL PG:'dbname=mygis host=localhost user=postgres' .\gadm36.gpkg

报如下错误：

ERROR 1: Unable to find driver `PostgreSQL'.

解决方案，不使用 PostgreSQL 下的 ogr2ogr，更换 QGIS 下能顺利执行：

& 'D:\Program Files\QGIS 3.10\bin\ogr2ogr.exe' -f PostgreSQL PG:'dbname=mygis host=localhost user=postgres' .\gadm36.gpkg

执行完后会生成一张 gadm 表。再次进行查询：

SELECT * FROM public.gadm WHERE ST_Within(ST_SetSRID(ST_MakePoint(121.55223, 38.86758), 4326), geom);

可得：

使用阿里 GeoJson 数据搭建反地理系统

前面讲到的国外地图数据存在的问题，不建议使用。今天要介绍的方法是使用国内阿里巴巴提供的全国城市边界数据。方法非常简单便捷，具体流程：