搭建自己的反地理查询系统

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

GADM

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

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

警告:GADM提供的中国国界数据不符合中国的领土主张,省界、市界、区界等数据也不一定是最新的版本。在正式刊物中发表使用此类数据的图件时需格外谨慎。

GADM提供了两种下载方式:

由于全球数据量巨大,建议根据需要按照国家下载数据。需要说明的是,GADM 中对country 的定义为“any entity with an 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运算就会形成一个新的多边形,这个新的多边形的边界为两个多边形中最大边界。

相关链接:

项目实战

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"的字符0x0x82 0x27在编码"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);

可得:

我来评几句
登录后评论

已发表评论数()

相关站点

+订阅
热门文章