修改时间: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 中的所有地理处理工具。打开软件能看到有很多地理处理工具:

image-20230824120611163

(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:

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

(2)示例:使用缓冲区buffer

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

缓冲工具图示

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

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

  3. 可以看到缓冲工具有四个参数,"{}"中的参数是可选的,可以不填或者用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 中,为什么统一地图标准,通常使用两种坐标系:

  1. 地理坐标系(Geographic Coordinate System):
    地理坐标系使用经度(Longitude)和纬度(Latitude)来表示地球上的位置。经度表示东西方向上的位置,纬度表示南北方向上的位置。常用的地理坐标系包括经度-纬度坐标系,如 WGS84(世界大地测量系统1984), GCJ02(火星坐标系),BD09(百度坐标系)等。前者是目前GPS使用的坐标系,后两者是国内使用常使用的坐标系,被加密,WGS84转后者可以使用百度或高德提供的地图转换服务,反过来转为WGS84需要用单独的方法。此处有吐槽。
  2. 投影坐标系(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.Projectarcpy.management.ProjectRaster,如果还未定义投影需要用定义投影工具。

img

(3)常用的投影坐标系

高斯-克吕格投影

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

img高斯-克吕格投影有着中央经线上无变形,满足投影后长度比不变的条件的优点,是国内坐标(北京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)哪些情况需要采用投影坐标系

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

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

(6)处理CAD投影

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

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

(5)WKID代号的查询

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

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

  1. geographic_coordinate_systems.pdf
  2. projected_coordinate_systems.pdf

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

image-20230812074808315


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

image-20230824132330857

image-20230824132354731

(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文件创建成功

以上示例通过遍历列表的方式可以将重复性的导入工作自动化,省时省力~


更多:

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

【ArcGIS Python系列】教程部分:

【ArcGIS Python系列】jupyter notebook:


如果你觉得本系列文章有用,欢迎关注博客,点赞和收藏,也欢迎在评论区讨论:

更多账号