分享

利用ArcGIS Python批量处理地理数据的坐标系

 晴耕雨读天 2023-09-09

本文整理自使用 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时的工作流程。

虽然软件提供了批量投影 工具,但是也不能同时批量处理所有文件,有时候还会遇见有的要素类没有投影,还需要先定义投影,然后才能进行投影工作。

  1. 检查坐标系,并决定是否需要运行投影工具。投影必须进行投影的要素类。在不进行投影的情况下复制其他要素类。

  2. 创建一个要素类。

  3. 检查图层和地图的坐标系:右键选择属性-选择源-空间参考-查看坐标系:有的为NAD 1983,有的为WGS 1984。

    图片例如,bike_racks 要素类位于名为 WGS 1984 的地理坐标系中,而 roads 和 boundary 要素类位于经过投影的美国国家平面坐标系(马里兰州)NAD 1983 StatePlane Maryland FIPS 1900 (US Feet) 中。

  4. 选择投影工具-选择输入要素类-选择投影坐标系-确定-等待完成。

图片
投影工具

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 : 投影成功

更多:

【ArcGIS Python系列】系列笔记为学习ArcGIS Pro和ArcPy过程中的总结,记下来方便回看,最新版本会优先发布在我的博客和GITHUB。

【ArcGIS Python系列】教程部分:

  • [x] 一、Arcpy介绍和安装【ArcGIS Python系列】

  • [x] 二、ArcGIS Pro和ArcMap的区别【ArcGIS Python系列】

  • [x] 三、Arcpy基础【ArcGIS Python系列】

  • [x] 四、探索空间数据【ArcGIS Python系列】

  • [x] 五、处理地理数据异常【ArcGIS Python系列】

  • [x] 六、处理几何数据【ArcGIS Python系列】

  • [x] 七、处理栅格数据【ArcGIS Python系列】

  • [x] 八、制图模块【ArcGIS Python系列】

  • [x] 九、自定义工具箱【ArcGIS Python系列】

  • [x] 十、ArcGIS_Pro常见问题【ArcGIS Python系列】

【ArcGIS Python系列】jupyter notebook:

  • [x] 4.2.3-arcpy基础(代码练习).ipynb
  • [x] 4.2.4-探索空间数据(代码练习).ipynb
  • [x] 4.2.5-示例1:使用Arcpy进行GIS人口空间分布数据探索.ipynb
  • [x] 4.2.7-处理几何数据代码练习和示例2.ipynb
  • [x] 4.2.8-栅格数据(代码练习).ipynb
  • [x] 4.2.12-实操1-利用ArcGIS_Python批量处理地理数据的坐标系.ipynb
  • [x] 4.2.13-实操2-使用 Python 对图像中的洪水进行分类.ipynb
  • [x] 4.2.14-实操3-制作考虑路况的交通等时圈.ipynb

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多