利用ArcGIS Python批量处理地理数据的坐标系
利用ArcGIS Python批量处理地理数据的坐标系
本文整理自使用 Python 自动化地理处理工作流，原教程为网页，我整理为了
jupyter notebook，加上了代码注释，方便学习和实操，难度1颗星，适合初学者，大佬请跳过。jupyter notebook地址：4.2.12-实操1-利用ArcGIS_Python批量处理地理数据的坐标系.ipynb 项目资源文件：PythonWorkflow.zip 由于微信的限制，可以点击最下方的阅读原文去我的博客阅读。
试想一下，你是一名GIS工作新人，你的领导总是让你做一些基础的工作，这一次他交给你政府和甲方提供的shp格式用地数据、兴趣点数据、街道数据等，你需要将分散在各个文件夹的一些数据集转换为统一的坐标系，然后将其导入到地理数据库中。你可能会在ArcGIS Pro中手动完成这些工作，但是如果你需要重复这些工作，那么手动完成这些工作就会变得很繁琐。在这种情况下，你可以使用30行Python代码自动化完成这些工作流程，然后你就可以摸鱼了。。。
1.数据准备
本次演示文件在"Workflow"文件夹中，为了便于理解，只包含一个地理数据库，数据库中的要素类和交通运输有关，我们需要将要素类需要存储在相同的同一个要素数据集中，并且具有相同的坐标系。
在处理 Python 代码之前，用ArcGIS打开名为Workflow.aprx的工程文件，检查一下数据库结构，如何所示：
数据库结构
2.手动流程
试想以下手动执行检查并统一坐标系的流程：检查要素类的坐标系，然后手动运行投影工具，借此能了解使用Arcpy时的工作流程。
虽然软件提供了
批量投影工具，但是也不能同时批量处理所有文件，有时候还会遇见有的要素类没有投影，还需要先定义投影，然后才能进行投影工作。
- 检查坐标系，并决定是否需要运行投影工具。投影必须进行投影的要素类。在不进行投影的情况下复制其他要素类。
- 创建一个要素类。
- 检查图层和地图的坐标系：右键选择属性-选择源-空间参考-查看坐标系：有的为NAD 1983，有的为WGS 1984。
例如，bike_racks 要素类位于名为 WGS 1984 的地理坐标系中，而 roads 和 boundary 要素类位于经过投影的美国国家平面坐标系（马里兰州）NAD 1983 StatePlane Maryland FIPS 1900 (US Feet) 中。
- 选择投影工具-选择输入要素类-选择投影坐标系-确定-等待完成。
投影工具
3.自动化流程
为了不重复这些步骤，我们用 Python 代码自动化完成此过程。
为了统一坐标系，我使用 Python 代码检查要素类的坐标系，并使用投影工具对所有当前不在正确坐标系中的数据集进行转换，从而将其复制到新地理数据库和要素数据集。最后可以通过此操作创建一个网络，但是本文不做演示。
（1）检查坐标系
我们在不打开软件的情况下检查数据库中所有要素类的坐标：
导入和环境设置：
import os
import arcpy
# 设置工作空间
mypath = os.path.join(os.getcwd(), r"resource\PythonWorkflow") # 修改为你的工作目录
gdb = "Transportation.gdb"
arcpy.env.workspace = os.path.join(mypath, gdb)
# print(arcpy.env.workspace)查看数据库中的要素类：
fcs = arcpy.ListFeatureClasses() # 获取所有要素类
len(fcs), fcs
10 ['bus_lines', 'boundary', 'street_lights', 'roads', 'traffic_analysis_zones', 'bike_routes', 'bike_racks', 'parking_zones', 'intersections', 'station_entrances']
fcs是要素类名称的 Python 列表。列表使用方括号括起来，而要素类名称为 Python 字符串，使用逗号分隔。
我们先检查一个要素类的投影信息：
desc = arcpy.da.Describe(fcs[0]) # 获取第一个要素类的描述信息
sr = desc["spatialReference"] # 获取要素类的空间参考
sr
sr是空间参考的 Python 字典。字典使用花括号括起来，而键(key)和值(value)之间使用冒号分隔。例如，键name对应于空间参考的名称，而键factoryCode对应于空间参考的WKID代码。
接下来可以通过for循环查看数据库中所有要素类的投影信息：
# 我们检查所有的要素类的坐标系
fcs = arcpy.ListFeatureClasses() # 获取所有要素类
# 遍历所有要素类
for fc in fcs:
desc = arcpy.da.Describe(fc) # 获取要素类的描述信息
sr = desc["spatialReference"] # 获取要素类的空间参考
print(fc + " : " + sr.name ) # 打印要素类名称和空间参考名称
bus_lines : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet
boundary : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet
street_lights : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet
roads : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet
traffic_analysis_zones : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet
bike_routes : GCS_WGS_1984
bike_racks : GCS_WGS_1984
parking_zones : GCS_WGS_1984
intersections : GCS_WGS_1984
station_entrances : GCS_WGS_1984
这些要素类中有一半是"NAD_1983_StatePlane_Maryland_FIPS_1900_Feet"，另一半是"GCS_WGS_1984"。
下一步就是将GCS_WGS_1984的所有要素类投影到一个坐标系中。但是，在此之前我们创建一个新的地理数据库用于储存投影后的要素。
（2）创建新的地理数据库
new_gdb = "Metro_Transport.gdb"
arcpy.env.overwriteOutput = True # 允许覆盖输出
new_gdb_path = arcpy.CreateFileGDB_management(mypath, new_gdb) # 创建新的地理数据库 返回值为新的地理数据库的对象, 可以作为工作空间引用
# 可以加一个数据库判断 避免重复操作 也可以用于检查数据库是否存在
# 判断数据库是否存在，如果不存在则创建
if not arcpy.Exists(os.path.join(mypath, new_gdb)):
arcpy.CreateFileGDB_management(mypath, new_gdb) # 创建新的地理数据库
else:
print("数据库已存在")
>>> 数据库已存在
（3）投影
接下来我们进行投影操作，投影后的要素类存储在新的地理数据库中。
首先确定我们想要投影到的坐标系，本次会使用roads要素类的投影坐标系NAD_1983_StatePlane_Maryland_FIPS_1900_Feet。使用坐标系的名字会很冗长，我们会使用WKID代码来代替坐标系的名字，WKID代码是唯一的，可以代表坐标系。
国内一般使用CGCS2000坐标系，也有使用百度和高德的坐标系的，具体使用哪个坐标系需要根据实际情况而定。
除了通过WKID代码设置坐标系，另一种方法是提供现有要素类的路径，例如，提供 Transportation.gdb 地理数据库中 roads 要素类的路径。
# 查询roads坐标系NAD_1983_StatePlane_Maryland_FIPS_1900_Feet的WKID代码
desc = arcpy.da.Describe("roads")
sr = desc["spatialReference"] # 获取要素类的空间参考
sr.factoryCode # 获取坐标系的WKID代码
>>> 2248
# 将WKID代码储存到变量
out_wkid = sr.factoryCode
# 同样我们创建一个空的要素数据集
fds = "Metro_Network" # 要素数据集名称
arcpy.CreateFeatureDataset_management(new_gdb_path, fds, out_wkid) # 创建要素数据集
如下"Messages"则为运行成功：
Messages
根据条件复制或投影要素类到新要素集复制要素工具arcpy.CopyFeatures_management和投影工具arcpy.Project_management都会使用一个输入要素类并生成一个输出要素类。虽然要素类的名称可以保持相同，但输出的路径将有所不同，因为新的要素类将位于新的地理数据库中。
fcs = arcpy.ListFeatureClasses() # 获取所有要素类
# 遍历所有要素类
for fc in fcs:
desc = arcpy.da.Describe(fc) # 获取要素类的描述信息
sr = desc["spatialReference"] # 获取要素类的空间参考
new_fc = os.path.join(mypath, new_gdb, fds, fc) # 新要素类的路径
if sr.factoryCode == out_wkid: # 如果要素类的WKID代码为2248 则进行复制 否则投影
arcpy.CopyFeatures_management(fc, new_fc) # 复制
print(fc + " : " + sr.factoryCode + " : " + "复制成功")
else:
arcpy.Project_management(fc, new_fc, out_wkid) # 投影
print(fc + " : " + sr.factoryCode + " : " + "投影成功")
bus_lines : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet : 复制成功
boundary : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet : 复制成功
street_lights : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet : 复制成功
roads : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet : 复制成功
traffic_analysis_zones : NAD_1983_StatePlane_Maryland_FIPS_1900_Feet : 复制成功
bike_routes : GCS_WGS_1984 : 投影成功
bike_racks : GCS_WGS_1984 : 投影成功
parking_zones : GCS_WGS_1984 : 投影成功
intersections : GCS_WGS_1984 : 投影成功
station_entrances : GCS_WGS_1984 : 投影成功