全国县市边界数据的处理(1M2017版)

全国地理信息资源目录服务系统上可以下载最新版本(2017年,1:1百万比例尺)的全国各要素基础数据,具体包括:

  • 共77幅DLG数据集,包括全国陆地范围以及台湾岛、海南岛、钓鱼岛、南海诸岛在内的主要岛屿及其临近海域。
  • 包括水系、居民地及设施、公路铁路、行政边界、地名等5个数据集12个数据层。

我们现在就需要其中的行政边界数据,本来以为下载后可以很方便的使用(有点想当然了),经过一番艰苦的处理后,数据终于是可用状态了。

在mac系统下进行了数据处理,使用了gdal/ogr、QGIS3,主要处理过程记录如下。

  1. 将下载下来的77副gdb数据,全部解压到一个目录中,比如1M2017目录,并在其中建立一个shp子目录。
  2. 在shell下,使用ogr2ogr命令,将所有gdb文件转换为一个大的shape文件,注意需要将其中的NAME字段编码转换为UTF8编码(官方是GBK编码)。
  3. 使用ogr2ogr,对合并后的BOUA(面状边界数据)文件根据PAC字段进行再次图形合并(因为默认的分幅会保留大量的正方形),生成output1.shp文件。
  4. output1.shp文件包含了大量了国外数据、海洋数据,这些数据可以通过ogr2ogr再次剔除,即设置PAC字段>100000且不等于250100,然后生成output2.shp文件。
  5. 仔细观察output2.shp文件,发现其中存在大量的拓扑错误(分幅边线的交叉遗留),在QGIS3中使用v.clean工具进行错误检查和处理,生成output3.shp文件。
  6. 对output3.shp文件进行人工错误检查(可借助QGIS3的图形检查工具check validity辅助检查),使用编辑工具(以及高级数字化工具)进行错误修复,主要使用节点编辑工具(删除节点)、多边形合并工具、多边形part删除工具等,最终生成全国(含港澳台)的县级边界数据集boua_cnty.shp。
  7. 使用QGIS3在boua_cnty数据中添加一个虚拟字段pac_city,并将其值定义为floor(PAC/100),然后对其执行dissolve操作,若需要将直辖市也当做一个独立城市边界,则需要单独处理一下四个直辖市的属性(将县区进行合并,主要是天津和重庆),形成全国市级边界数据集boua_city.shp文件。
for i in `ls -d *.gdb`;do ogr2ogr -f "ESRI Shapefile" -lco ENCODING=UTF-8  -append  shp $i;done
ogr2ogr -f "ESRI Shapefile" output1.shp -lco ENCODING=UTF-8 -dialect sqlite -sql "select st_union(geometry) as geometry,PAC,NAME from 'BOUA' group by PAC,NAME order by PAC,NAME" shp/BOUA.shp
ogr2ogr -f "ESRI Shapefile" output2.shp -lco ENCODING=UTF-8 -where "PAC>100000 and PAC<>250100" output1.shp

但这样生成的城市边界仅包括陆地边界,对于内陆型城市,在使用中问题不大,但对于沿海城市,则很多人类活动轨迹可能会超出陆地边界(比如舟山),因此还需要一个融合近海边界的城市边界。手头上有一个2012年版本的含海洋边界的城市边界数据,需要将其海洋边界数据融合到新的全国边界数据集中。

处理过程为:

  1. 对含有海洋边界的数据集bou2012,和最新全国城市边界数据集boua_city,执行difference操作,在QGIS3进行编辑,删除内陆边界中有差异的部分,在沿海边界删除海洋中的小岛(内多边形),同时添加一个虚拟字段pac_city(和boua_city中的字段相同),最终生成海洋城市边界数据集。
  2. 将boua_city数据集和上一步中生成的海洋城市边界数据集执行union操作。
  3. 对union后的数据集,根据pac_city字段执行dissolve操作。
  4. 人工检查最终生成的数据集。

因数据问题,此数据集中不含港澳台的海洋边界。

利用shell命令合并分县土地利用数据为一个shape File

为了制图服务的需求,需要将分县的土地利用arcinfo coverage数据合并为一个shape File。
用ogrinfo查看一个县的数据简要信息:
ogrinfo -so id04632323 PAL
Had to open data source read-only.
INFO: Open of `id04632323'
using driver `AVCBin' successful.

Layer name: PAL
Geometry: Polygon
Feature Count: 4065
Extent: (-831725.500000, 2525475.000000) - (-704241.375000, 2614980.250000)
Layer SRS WKT:
PROJCS["unnamed",
GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",
DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",
SPHEROID["Krassowsky 1940",6378245,298.3,
AUTHORITY["EPSG","7024"]],
AUTHORITY["EPSG","6024"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4024"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",25],
PARAMETER["standard_parallel_2",47],
PARAMETER["latitude_of_center",12],
PARAMETER["longitude_of_center",110],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["METERS",1]]
ArcIds: IntegerList (0.0)
AREA: Real (12.3)
PERIMETER: Real (12.3)
ID04632323#: Integer (5.0)
ID04632323-ID: Integer (5.0)

用QGIS打开后,可以发现其ID04632323-ID字段是土地利用的分类代码。而要将这些分县的数据合并,则需要利用ogr2ogr来合并生成。
主要问题:
1. 需要重命名字段,如将这个分类代码重命名为一个共用的名称,比如code。
2. 目录名和分类代码字段名相关,但需要进行大小写转换。
3. 还需要将每个县的目录名称保存到对应的数据库中。

在bash中,大小写转换可以用这个函数,将变量$i转换为大写:
${i^^}

ogr2ogr可以用sql选项来提取不同的字段并重命名:
select AREA as area,PERIMETER AS perimeter,ID04632323-ID as code from PAL
但在测试中,发现其不认这种写法,需要将其用单引号来括起来:
select AREA as area,PERIMETER AS perimeter,‘ID04632323-ID’ as code from PAL

而目录名称也可以通过变量的方式传递到sql中来。

最终命令脚本,要注意其中的引号和转义:
for i in id*;do
ogr2ogr -update -append -sql "select AREA as area,PERIMETER AS perimeter,'${i^^}-ID' as code,\"$i\" as county FROM PAL" landuse.shp $i;
done