三、Arcpy基础【ArcGIS Python系列】

修改时间：2023 年 8 月 24 日

阅读时间：15~20 分钟

在阅读本篇和后续篇章时，推荐要有一点点 Python 基础。

代码练习 notebook 为4.2.3-arcpy 基础(代码练习).ipynb

1.地理数据处理基础知识

ArcGIS 中的地理处理允许您执行空间分析和建模以及自动执行 GIS 任务。典型的地理处理工具获取输入数据（要素类、栅格或表），执行地理处理任务，然后生成输出数据作为结果。ArcGIS 包含数百种地理处理工具。地理处理工具的示例包括用于创建缓冲区、用于向表添加字段以及用于对地址表进行地理编码的工具。

地理处理通过创建组合了一系列工具的序列来支持工作流的自动化。一个工具的输出实际上成为下一个工具的输入。通过使用模型（model builder）和脚本，可以在 ArcGIS 中地理处理工具的自动化工作流。

2.从导入 ArcPy 开始

ArcPy 包含许多模块、类和函数，这使得可以在 Python 脚本中使用 ArcGIS Pro 中的所有地理处理工具。

import arcpy

执行上述语句后，就可以运行随 ArcGIS Pro 安装的工具箱中的所有地理处理工具。包括用于处理数据的模块 (arcpy.da)、地图脚本模块 (arcpy.mp)、用于图像分析和解释的模块 (arcpy.ia) 以及用于地图代数和栅格分析的模块（arcpy.sa）。导入 ArcPy 后，您就可以开始使用其模块、函数和类。

在脚本中导入 ArcPy 不仅会导入 ArcPy 的功能，还会执行两项重要检查：ArcPy 的可用性和许可证的可用性。如果输出RuntimeError: NotInitialized错误消息，请参照4.2.1-arcpy 介绍和安装.md安装 ArcGIS Pro。

3.设置工作区

首先得理解 Python 中绝对路径和相对路径的区别，简单提示一下：

绝对路径（Absolute Path）是从文件系统的根目录开始的完整路径。 它包含了从根目录到目标文件或目录的所有目录层级。在不同的操作系统中，根目录的表示方式可能不同。例如，在 Windows 系统中，绝对路径可以以盘符（如 C:\）开始，而在 Linux 或 Mac 系统中，绝对路径以斜杠（/）开始。在代码中如果是反斜杠 “"，应该改为 “/”(正斜杠)或’'\'(两个反斜杠)。或者写成 r"C:/data1”。
相对路径（Relative Path）是相对于当前工作目录的路径。 当前工作目录是指运行 Python 程序时所在的目录。相对路径指定了从当前工作目录到目标文件或目录的路径。相对路径可以是简单的文件名或目录名，也可以是包含目录层级关系的路径。

ArcPy 中的工作空间指定的就是工作目录，对应的可以使用相对路径引用。独立的 Python 脚本默认情况有一个当前工作目录，默认情况下该目录是脚本的位置。当设置 arcpy.env.workspace 时，ArcGIS Pro 将会在该路径下查找和操作数据。

tip: 您可以使用 os.getcwd() 获取当前工作目录，并且可以使用 os.chdir("/path") 更改当前工作目录。这样我们就能够在工作目录中使用相对路径指定路径了，保证了代码的可移植性。

arcpy.env.workspace 本质是一个 Python 类。

注意理解env 是一个 Python 类（class），workspace 是该类的一个属性（property）。arcpy.env.workspace对应arcpy.<class>.<property>，所以arcpy.<classname>.<property> = <value>就是工作空间的属性值。

例如，你有一个名为 “C:\Data” 的文件夹，其中包含了你要使用的地理数据，你可以通过以下方式将它设置为工作空间：

import arcpy

arcpy.env.workspace = r"C:\Data"

# 创建地理数据库
arcpy.CreateFileGDB_management(arcpy.env.workspace, "myGDB.gdb") # 在工作空间下创建名为myGDB.gdb的地理数据库

在这个例子中，arcpy.env.workspace 被设置为 “C:\Data”，这意味着在执行地理处理脚本时，ArcGIS Pro 将会在该文件夹下查找和操作数据。

使用 arcpy.env.workspace 的好处是，它可以确保地理处理脚本在不同的环境中都能正常工作，无论是在 Windows 还是其他操作系统上。它提供了一种统一的方式来设置工作空间，使得脚本可以在不同的计算机上或不同的工作目录中运行，而不需要手动更改路径。（此方法和 python 的相对路径的作用相同）例如你可以这样指定工作空间：

import os

# 在整个脚本前指定一次绝对路径
data_dir = r'C:\Users\<用户名>\Documents\Python_\Github\arcgis-notebooks-tutorial\hurricane_analysis\data'
# 以后路径都是用相对路径 利用os.path.join处理路径能避免许多问题
hurricanes_raw_dir = os.path.join(data_dir,'hurricanes_raw')

# 利用mkdir创建检查和创建目录
if not os.path.exists(hurricanes_raw_dir):
    os.mkdir(os.path.join(data_dir,'hurricanes_raw'))

总而言之，arcpy.env.workspace 是一个用于设置地理处理脚本工作空间的变量，它确保脚本能够在不同的环境中正确访问和操作数据。

4.使用地理处理工具

ArcPy 使您可以访问 ArcGIS Pro 中的所有地理处理工具。打开软件能看到有很多地理处理工具：

（1）调用工具的方法两种方法

作为 python函数调用：arcpy.<toolname_toolboxalias>(<parameters>)

例如调用裁剪工具：

import arcpy
arcpy.env.workspace = "C:/Data"
arcpy.Clip_analysis("streams.shp", "study.shp", "result.shp")

还可以通过使用与工具箱别名匹配的模块来使用工具。首先将工具箱作为模块调用，然后将工具作为该模块的函数调用，最后是工具的参数。

arcpy.<toolboxalias>.<toolname>(<parameters>)

import arcpy
arcpy.env.workspace = "C:/Data"
arcpy.analysis.Clip("streams.shp", "study.shp", "result.shp")

两种方法都是正确的。

小 tips：

Python 区分大小写，Clip 不等于 clip
在代码行中空格对执行没有影响，但是对可读性有影响 👀。
函数和变量之间不要有空格，() 正确， <toolname> (<parameters>)不正确。
arcpy .analysis. Clip()不产生错误但是不是正确格式。

（2）示例：使用缓冲区 buffer

而具体函数的使用可以参照帮助文档。我们从缓冲区 buffer的帮助文档工具中举例说明：

缓冲工具图示

程序中通过搜索找到buffer工具，可以看到通过此地图处理工具的可视化操作的参数：带星号的是必填此参数，分别是输入要素、输出要素类和距离。

点击右上角的问号 ❔，进入帮助页面，找到 Python 代码:

可以看到缓冲工具有四个参数，"{}“中的参数是可选的，可以不填或者用 None 和”"表示, 其余的是必选参数。缓冲工具的语法是:

arcpy.analysis.Buffer(in_features, out_feature_class, buffer_distance_or_field, {line_side}, {line_end_type}, {dissolve_option}, {dissolve_field}, {method})

Python 函数中的参数和软件中的参数基本是一一对应的，前三个参数为必选参数。我们简单浏览整个表格后然后一一说明：

1）必选参数

名称	说明	数据类型
in_features	要进行缓冲的输入点、线或面要素。	Feature Layer
out_feature_class	包含输出缓冲区的要素类。	Feature Class
buffer_distance_or_field	与要缓冲的输入要素之间的距离。该距离可以用表示线性距离的某个值来指定，也可以用输入要素中的某个字段（包含用来对每个要素进行缓冲的距离）来指定。如果未指定线性单位或输入了“未知”，则将使用输入要素空间参考的线性单位。指定距离时，如果所需线性单位含有两个单词，如 Decimal Degrees，请将两个单词合并成一个词（例如，20 DecimalDegrees）。	Linear Unit; Field

2）可选参数：

名称	说明	数据类型
line_side(可选)	指定将在输入要素的哪一侧进行缓冲。 *该参数仅支持面和线要素。1. FULL—对于线，将在线两侧生成缓冲区。对于面，将在面周围生成缓冲区，并且这些缓冲区将包含并叠加输入要素的区域。这是默认设置。2. LEFT—对于线，将在线的拓扑左侧生成缓冲区。此选项不支持面输入要素。3. RIGHT—对于线，将在线的拓扑右侧生成缓冲区。此选项不支持面输入要素。4. OUTSIDE_ONLY—对于面，仅在输入面的外部生成缓冲区（输入面内部的区域将在输出缓冲区中被擦除)	String
line_end_type(可选)	指定线输入要素末端的缓冲区形状。此参数对于面输入要素无效。ROUND—缓冲区的末端为圆形，即半圆形。这是默认设置。FLAT—缓冲区的末端很平整或者为方形，并且在输入线要素的端点处终止。	String
dissolve_option(可选)	指定移除缓冲区重叠要执行的融合类型。NONE—不考虑重叠，将保持每个要素的独立缓冲区。这是默认设置。ALL—将所有缓冲区融合为单个要素，从而移除所有重叠。LIST—将融合共享所列字段（传递自输入要素）属性值的所有缓冲区。	String
dissolve_field[dissolve_field,…]（可选）	融合输出缓冲区所依据的输入要素的字段列表。将融合共享所列字段（传递自输入要素）属性值的所有缓冲区。	Field
method(可选)	指定是使用平面方法还是测地线方法来创建缓冲区。PLANAR—如果输入要素位于投影坐标系中，则将创建欧氏缓冲区。如果输入要素位于地理坐标系中且缓冲距离的单位为线性单位（米、英尺等，而非诸如度之类的角度单位），则会创建测地线缓冲区。这是默认设置。您可以使用输出坐标系环境设置指定要使用的坐标系。例如，如果输入要素位于投影坐标系中，您可以将环境设置为地理坐标系，以便创建测地线缓冲区。GEODESIC—无论使用哪种输入坐标系，均使用形状不变的测地线缓冲区方法创建所有缓冲区。	St

看起来很复杂，如果你了解 ArcGIS Pro 中使用缓冲区的方法，其实只需要将相应参数填入函数内就可以了，例如：

import arcpy, os

# 用os.getcwd()设置文件夹
data1_dir = os.path.join(os.getcwd(), "resource/data1")

arcpy.env.workspace = data1_dir

# 将shp导入数据库
arcpy.CopyFeatures_management("streets.shp", os.path.join("demo.gdb", "streets"))

#  可以更改工作空间
arcpy.env.workspace = os.path.join(data1_dir, "demo.gdb") # 改成你的data1文件夹位置

# 缓冲区分析
arcpy.analysis.Buffer("streets", "streets_Buffered", "20 Meters", "FULL", "ROUND", "LIST", "LABEL_CLAS")

我们将输出的文件 streets_Buffered 拖入地图中：

以上也可以改写成以下形式方便阅读：

arcpy.analysis.Buffer(in_features="streets",
                      out_feature_class= "streets_Buffered_1",
                      buffer_distance_or_field="20 Meters",
                      line_side="FULL",
                      line_end_type="ROUND",
                      dissolve_option="LIST",
                      dissolve_field="LABEL_CLAS")

你也可以单独定义变量，方便代码复用和制作脚本：

in_features="streets"
out_feature_class="streets_Buffered_2"
buffer_distance_or_field="20 Meters"
line_side="FULL"
line_end_type="ROUND"
dissolve_option="LIST"
dissolve_field="LABEL_CLAS"

arcpy.analysis.Buffer(in_features,out_feature_class, buffer_distance_or_field, line_side, line_end_type, dissolve_option, dissolve_field)

（3）调用工具的小技巧

在软件中打开某个地理处理工具的界面，你可以点击右上角的问号进入帮助页面，在填好参数之后就可以将此操作复制为python命令，甚至你可以在运行完地理处理工具确保没有错误之后再复制为python命令。如图所示：

你也可以打开历史记录，复制你运行过的命令：

5.空间参考

在 ArcGIS 中，为什么统一地图标准，通常使用两种坐标系：

地理坐标系（Geographic Coordinate System）：
地理坐标系使用经度（Longitude）和纬度（Latitude）来表示地球上的位置。经度表示东西方向上的位置，纬度表示南北方向上的位置。常用的地理坐标系包括经度-纬度坐标系，如 WGS84（世界大地测量系统 1984）， GCJ02（火星坐标系），BD09（百度坐标系）等。前者是目前 GPS 使用的坐标系，后两者是国内使用常使用的坐标系，被加密，WGS84 转后者可以使用百度或高德提供的地图转换服务，反过来转为 WGS84 需要用单独的方法。此处有吐槽。
投影坐标系（Projected Coordinate System）：
投影坐标系是将地理坐标系映射到平面上的二维坐标系。它使用笛卡尔坐标系，其中位置由 X 和 Y 坐标表示。投影坐标系通常用于地图制作和空间分析。常见的投影坐标系包括通用横轴墨卡托投影（Universal Transverse Mercator，UTM）。

（1）理解空间参考类

我们通过空间参考类（SpatialReference）来指定和引用空间参考。一般在创建空白要素类的时候以及投影转换的时候使用。

此类具有多个属性，包括坐标系参数。但是，若要使用这些属性，必须实例化类（instantiated），需要为此类创建一个对象。类就像一个此对象的蓝图，你可以通过实例化类在此蓝图的基础上创建一个对象。例如我们查看 shpfile 系列文件中.prj文件来看看其空间参考属性：

import arcpy, os
prjfile = os.path.join(data1_dir, "streets.prj")
spatialref = arcpy.SpatialReference(prjfile)
spatialref

spatialref对象

上图可以spatialref的全部属性。我们也可以单独获取其属性，比如

spatialref.name # 获取空间参考文件的名称
# >>> 'NAD_1983_HARN_Adj_MN_Clay_Feet'

以下演示定义一个地理坐标系空间参考对象：

sr1 = arcpy.SpatialReference("GCS_WGS_1984")
sr2 = arcpy.SpatialReference(4326) # ！！！数字类型不是字符串

# 判断两个参考系是否相等
sr1 == sr2
>>> True # 证明相等

可以同时对空间参考对象定义地理坐标系和投影坐标系。

（2）投影的概念

投影是一种将地球表面上的三维地理坐标（经度、纬度和高程）映射到二维平面上的方法。由于地球是一个三维椭球体，将其映射到平面上会引入形状、距离和方向的变形。因此，选择适当的投影方法对于特定的地理区域和任务非常重要。ArcGIS 提供了各种投影方法，包括等距圆柱投影、等距圆锥投影、等面积投影、等角投影等。每种投影方法都有其特定的优势和适用范围。ArcGIS 中使用投影和投影栅格工具进行投影变换，对应的 Arcpy 方法是arcpy.management.Project和arcpy.management.ProjectRaster，如果还未定义投影需要用定义投影工具。

（3）常用的投影坐标系

高斯-克吕格投影

高斯-克吕格投影又称等角横轴切圆柱投影，即横轴墨卡托投影（Transerse Mercator)。

高斯-克吕格投影有着中央经线上无变形，满足投影后长度比不变的条件的优点，是国内坐标（北京 54 坐标、西安 80 坐标、CGCS2000 坐标）最常采用的投影坐标系。我国各种大、中比例尺地形图采用了不同的高斯-克吕格投影带。其中大于 1:1 万的地形图采用 3° 带；1:2.5 万至 1:50 万的地形图采用 6° 带。
高斯-克吕格投影分为 3° 分带和 6° 分带两种分带方法。如下图所示：

在这里插入图片描述

3° 分带法：
从东经 1°30′起，每 3° 为一带，将全球划分为 120 个投影带；
3° 分带常应用于大比例尺地形图，大于 1：1 万的地形图均采用 3° 分带，城建坐标多采用 3° 分带。

6° 分带法：
从 0° 经线（格林威治）起，每 6° 分为一个投影带，全球共分为 60 个投影带；
6° 分带常应用于小比例尺地形图，包括 1：100 万、1：50 万、1：25 万、1：10 万、1：5 万、1：2.5 万地形图。

在 Arcgis 中，例如CGS2000 3 Degree GK Zone 38用来表示高斯-克吕格投影坐标：

CGCS2000代表该投影坐标系基于的地理坐标系为 CGCS2000

3 Degree代表 3° 分带，当没有标注 3 Degree 时就代表是 6° 分带

GK 代表高斯-克吕格投影方

Zone 代表投影带

UTM 投影

通用墨卡尔投影是英、美、日、加拿大等国地形图最通用的投影。简称“UTM 投影”。

（4）哪些情况需要采用投影坐标系

选择投影坐标系解决地球表面的曲面到平面的映射问题，避免地球曲面产生的误差。以下情况需要使用投影坐标系：

地图制作：当需要制作地图时，通常需要将地球表面的曲面映射到平面上。由于地球是一个三维椭球体，直接在平面上表示地球上的地理坐标会引入形状、距离和方向的变形。通过采用适当的投影坐标系，可以将地理坐标转换为平面坐标，以在地图上准确地表示地理特征、距离和方向。
空间分析：在进行空间分析时，需要进行地理数据的测量、叠加和分析。在地理坐标系下，直接进行距离、面积和方向的计算可能不准确，因为地球的曲面会引入误差。通过将数据转换到适当的投影坐标系，可以进行准确的空间分析，确保测量和计算的精度。
数据叠加：当需要将来自不同数据源的地理数据进行叠加时，这些数据可能使用不同的地理坐标系。为了进行准确的叠加，您需要将数据转换到相同的投影坐标系，以确保它们在平面上的位置和几何关系正确匹配。
可视化和展示：在将地理数据可视化和展示时，使用投影坐标系可以确保地图的形状和比例符合实际。通过选择适当的投影坐标系，可以在地图上准确地显示地理特征和空间分布，使观众能够更好地理解和解读地理信息。

（6）处理 CAD 投影

通常，城市设计和建筑设计一类的数据都是从 CAD 导入坐标点，当你直接导入 Arcgis 中，此时会提示该文件无空间参考信息，此时需要先知晓 CAD 的空间数据之后，通过定义投影来确定导入 CAD 的空间参考信息。

如果连 CAD 的空间参考都不知道，就只能一个个空间参考去和卫星图对比了，可以先按照 CGCS2000、西安 80、北京 54 坐标系的顺序来尝试，直到选出合适的坐标系，祝你好运。如何实在是选不出坐标系，可以尝试地理配准工具。

（5）WKID 代号的查询

WKID 空间参考代号查询方式：

通过查询空间参考代号，可以找到对应的空间参考文件，从而确定空间参考。空间参考代号可以通过以下两个文件查询：

比如我们要查询 WGS1984 的 WKID 和标准命名，打开文件 1 搜索：

也可以在 arcgis pro 中点开任何一个图层或者使用打开投影工具，进行可视化操作：

（6）投影相关函数的使用

简单来说，对于矢量数据采用投影arcpy.management.Project，对于栅格数据采用投影栅格arcpy.management.ProjectRaster，如果没有数据空间参考采用定义投影arcpy.management.DefineProjection。他们都可以传入空间参考类的实例化对象作为参数传入，拿定义投影举例：

import arcpy
infc = r"C:\data\demo.shp"
sr = arcpy.SpatialReference(4326) # 建议输入WKID
arcpy.DefineProjection_management(infc, sr)

6.示例

我们看个简单的 Python 例子，如果要将文件夹中所有的 shapefile 复制到地理数据库，我们应该怎么办：

# 1导入包
import arcpy, os

# 2定义相关参数
gdb = "demo.gdb"
workspace = r"Z:\Sync\Urban-Spatial-Data-Analysis-For-Beginners\4-空间数据分析\4.2-arcpy\resource\data1"
arcpy.env.workspace = workspace

# 3创建地理数据库
arcpy.CreateFileGDB_management(workspace, gdb)

# 列出文件夹中的要素类 shp文件
fc_list = arcpy.ListFeatureClasses()

for fc in fc_list:
    fc_desc = arcpy.Describe(fc)
    if fc_desc.shapeType == "Polyline":
        newfc = os.path.join(gdb, "Polyline",fc_desc.basename)

prj = arcpy.SpatialReference(os.path.join(workspace, "streets.prj")) # 读取shp文件的投影信息
arcpy.CreateFeatureDataset_management(gdb, "Polyline", prj) # 在数据库中创建名叫Polyline的空白要素类 指定其空间参考为 "streets.prj"的空间参考

arcpy.CopyFeatures_management(fc, newfc)

你可以将以下代码按照我的分块形式复制到 jupyter notebook 中执行，以便更清楚的了解每行代码发生了什么：

在第 3 步代码运行之后，你会发现 data1 文件夹下多了一个空的 gdb 数据库：
生成结果

第 4 步我们想把 data1 文件夹里所有(其实只有一个)多段线的要素导入到此数据库，首先列出当前工作空间的要素类：

fc_list = arcpy.ListFeatureClasses()
fc_list
"""
>>> (">>>"我用来表示notebook单元格的输出，下文不再另做说明)
>>> ['parcels.shp', 'streets.shp']
"""

然后把这两个要素类通过数据arcpy.Describe返回的对象中的数据类型shapeType进行判定，如果是多段线则构建输出文件名。

for fc in fc_list:
    fc_desc = arcpy.Describe(fc)
    if fc_desc.shapeType == "Polyline":
        newfc = os.path.join(gdb, "Polyline",fc_desc.basename)
newfc
"""
>>> 'demo.gdb\\Polyline\\streets'
"""

我们复制此要素到数据库：

prj = arcpy.SpatialReference(os.path.join(workspace, "streets.prj")) # 读取shp文件的投影信息
arcpy.CreateFeatureDataset_management(gdb, "Polyline", prj) # 在数据库中创建名叫Polyline的空白要素类

arcpy.CopyFeatures_management(fc, newfc)

创建成功后此要素会被自动加载到 arcgis 的地图中：

streets文件创建成功