背景
数据来源:全国地理信息资源目录服务系统
这个网站是由自然资源部主管, 国家基础地理信息中心进行运行维护的,因此数据是比较权威的。目前,网站上提供了以下三类数据:
-
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中的展示结果。