最新数据处理 之 1:100万与1:25万全国基础地理数据库批量合并

背景

数据来源：全国地理信息资源目录服务系统

链接：全国地理信息资源目录服务系统

这个网站是由自然资源部主管， 国家基础地理信息中心进行运行维护的，因此数据是比较权威的。目前，网站上提供了以下三类数据：

• 30米全球地表覆盖数据，数据整体现势性达到2010年

• 1:100万全国基础地理数据库，数据整体现势性达到2015年

• 1:25万全国基础地理数据库，数据整体现势性达到2015年

个人可以进行注册，然后查询和下载，下载的时候，需要提供单位和下载的目的。提交后，很快就可以进行下载了。

但是，数据是基于块的，下载后，还需要进行逐一拼接，虽然用ArcGIS的工具可以实现，但手动点击操作比较繁琐。因此，本文主要介绍如何进行批量拼接，用到的是ArcPy工具。

数据组织形式如上图。由于1：100万的数据是采用gbd数据库的形式存储的，1：25万的数据是采用shape文件进行保存的，因此两个数据的拼接还需要注意。废话不多说，直接上代码。

1:100万全国基础地理数据库

import os
import arcpy
gdb_list = []
file_dir = "gdb"  # put the unzipped gdb file (such as A49.gdb) in this folder
dir_list = os.listdir(file_dir)
for dir in dir_list:
    if ".gdb" in dir:
        gdb_list.append(file_dir + "\\" + dir)
print(gdb_list)  # check the gdb
out_path = "dataset_100w"
out_gdb_name = "merge_100w"
merge_gdb = out_path + "\\" + out_gdb_name + ".gdb"
if arcpy.Exists(merge_gdb):
    pass
else:
    arcpy.CreateFileGDB_management(out_path, out_gdb_name)
for gdb in gdb_list:
    print(gdb)
    arcpy.env.workspace = gdb
    fcs = arcpy.ListFeatureClasses()
    for fc in fcs:
        print("        " + fc)
        if arcpy.Exists(merge_gdb + "\\" + fc):
            arcpy.Append_management(gdb + "\\" + fc,
                                    merge_gdb + "\\"  + fc,
                                    "NO_TEST")
        else:
            arcpy.FeatureClassToFeatureClass_conversion(fc, merge_gdb, fc)

1:25万全国基础地理数据库

import os
import arcpy
path = 'data'  # put the unzipped file in this folder
fileList = os.listdir(path)
a = []
b = []
c = []
d = []
e = []
f = []
g = []
h = []
k = []
for i in fileList:
    path_name = path + '\\' + i
    arcpy.env.workspace = path_name
    fsc = arcpy.ListFeatureClasses()
    for j in fsc:
        if j == "agnp.shp":
            a.append(path_name + "\\" + j)
        elif j == 'aanp.shp':
            b.append(path_name + "\\" + j)
        elif j == 'hyda.shp':
            c.append(path_name + "\\" + j)
        elif j == 'hydl.shp':
            d.append(path_name + "\\" + j)
        elif j == 'hydp.shp':
            e.append(path_name + "\\" + j)
        elif j == 'lrdl.shp':
            f.append(path_name + "\\" + j)
        elif j == 'lrrl.shp':
            g.append(path_name + "\\" + j)
        elif j == 'resa.shp':
            h.append(path_name + "\\" + j)
        elif j == 'resp.shp':
            k.append(path_name + "\\" + j)
arcpy.env.workspace = 'dataset_25w\\'
arcpy.env.overwriteoutput = True
arcpy.Merge_management(a,"agnp_25w.shp")
arcpy.Merge_management(b,"aanp_25w.shp")
arcpy.Merge_management(c,"hyda_25w.shp")
arcpy.Merge_management(d,"hydl_25w.shp")
arcpy.Merge_management(e,"hydp_25wn.shp")
arcpy.Merge_management(f,"lrdl_25w.shp")
arcpy.Merge_management(g,"lrrl_25w.shp")
arcpy.Merge_management(h,"resa_25w.shp")
arcpy.Merge_management(k,"resp_25w.shp")

最后

这两段代码经过了最新测试，没有出现问题，可以放心使用！其中，1：100万的结果合并成了gdb数据库格式，1：25万的结果合并成了shape格式，都可以用ArcGIS打开。下图是1：100万居民地地名在ArcGIS中的展示结果。