四、探索空间数据【ArcGIS Python系列】

renhai

发布于 2023-11-24 16:20:36

2830

发布于 2023-11-24 16:20:36

文章被收录于专栏：renhailab数据分析

微信公众号无法嵌入超链接，推荐去博客阅读。点击底部阅读原文获得更好的阅读体验。代码练习notebook :4.2.4-探索空间数据(代码练习).ipynb

本篇介绍了探索空间数据的几种方法。

1.检查数据是否存在

在 Python 脚本中，您可能需要确定数据集是否确实存在。该任务可以使用 arcpy.Exists()函数完成。返回函数返回一个布尔值True或False，指示该元素是否存在。

import arcpy
arcpy.env.workspace = os.path.join(os.getcwd(), "resource/data1", "demo.gdb")
print(arcpy.Exists("streets")) # 可以写绝对路径
>>> True

arcpy.Exists()与os.path.exists()的区别有两点：

第一是 arcpy.Exists()可以识别ArcGIS工作空间所在目录，os.path.exists()只能是系统目录。第二是 arcpy.Exists()可以识别数据库中的文件，而os.path.exists()不可以。

注意：ArcGIS Pro 不再支持使用个人地理数据库（文件扩展名 .mdb）。如果您有较旧的 .mdb 文件，则必须使用 ArcGIS Desktop 10.x 应用程序（例如 ArcMap 或 ArcCatalog）将它们转换为文件地理数据库 (.gdb)，然后才能在 ArcGIS Pro 中使用它们。

2.描述数据的两种方法

ArcPy 中有两个函数用于描述数据集：arcpy.Describe() 和 arcpy.da.Describe()。两者执行相同的任务，但功能结构略有不同。这两个函数都很常用，因此了解这两个函数是有好处的。

Describe()函数是ArcPy的常规函数，返回一个对象。da.Describe() 函数是 arcpy.da 模块或数据访问模块的函数，返回的值是字典。除了描述数据之外，此 ArcPy 模块还用于编辑任务并支持数据库工作流程，例如版本控制、副本、域和子类型。

desc = arcpy.Describe("streets")
print(desc.shapeType)
>>> Polyline

desc = arcpy.da.Describe("streets")
print(desc["shapeType"])
>>>Polyline

举例：

打印空间参考

fc = "streets"
desc = arcpy.da.Describe(fc)
sr = desc["spatialReference"]
print("空间参考是：" + sr.name)
>>> 空间参考是：NAD_1983_HARN_Adj_MN_Clay_Feet

3.列出数据

在脚本中列出数据，方便后续操作。ArcPy列表函数包括ArcPy列表函数包括 ListFields() 、 ListIndexes() 、 ListDatasets() 、 ListFeatureClasses() 、 ListFiles() 、 ListRasters() 、 ListTables() 、 ListWorkspaces() 和 ListVersions() 。

最常用的是arcpy.ListFeatureClasses()他会返回当前工作空间的要素类列表，用于检索数据：

arcpy.ListFeatureClasses({wild_card}, {feature_type},                          {feature_dataset})

wild_card指通配符限制列表，类似于git的忽略文件中的匹配模式，""代表的是通配符，等同字符串""。"s"可以匹配shanghai也可以匹配shenzhen。
feature type将限制结果的要素类型。

举例：

# 查看data1文件夹下有什么文件
arcpy.env.workspace = os.path.join(os.getcwd(), "resource/data1")
fclist = arcpy.ListFeatureClasses()
print(fclist)
# 返回列表元素
> ['parcels.shp', 'streets.shp']

与之类似的是返回当前工作空间中的数据集的列表ListDatasets。他们俩通常可以结合来用：

arcpy.env.workspace = os.path.join(os.getcwd(), "resource/data1", "demo.gdb")
datasets = arcpy.ListDatasets(feature_type='feature')
datasets = [''] + datasets if datasets is not None else []

# datasets为数据集Polyline，我们通过遍历查询此数据集内部的要素
for ds in datasets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
        path = os.path.join(arcpy.env.workspace, ds, fc)
        print(path)
"""
你会看到类似的输出：
~\Documents\Python_\Urban-Spatial-Data-Analysis-For-Beginners\4-空间数据分析\4.2-arcpy\resource/data1\demo.gdb\Polyline\streets
~\Documents\Python_\Urban-Spatial-Data-Analysis-For-Beginners\4-空间数据分析\4.2-arcpy\resource/data1\demo.gdb\streets_Buffered_1
~\Documents\Python_\Python辅助城市研究\Urban-Spatial-Data-Analysis-For-Beginners\4-空间数据分析\4.2-arcpy\resource/data1\demo.gdb\streets_Buffered_2
.......
"""

其他常用函数：

列出栅格：要在当前工作空间中创建栅格数据集列表，使用 arcpy.ListRasters() 函数。

列出字段：还有一个常用的列出数据的函数是 arcpy.ListFields() 。此函数列出指定数据集的要素类或表中的字段。语法为arcpy.ListFields(dataset, {wild_card}, {field_type})。

列出文件：另一个有用的函数是 arcpy.ListFiles() 。此函数返回工作区中所有文件的列表，通常用于列出不是空间数据集的独立文件，包括CSV文件（.csv）、Microsoft Excel文件（.xlsx）和纯文本文件（.txt）。ListFiles() 函数的一般语法为arcpy.ListFiles({wild_card})。

需要注意的是 ListFiles() 不显示子文件夹内和数据库中的内容，但是子文件夹和数据库本身会显示。如果需要遍历子文件夹，可以用da.Walk()，ArcPy的 da.Walk() 函数类似于内置 os 模块的 os.walk() 函数。还有一个ListDatasets()返回当前工作空间中的数据集的列表。举个例：

对于以上"Project"文件夹，我们使用要使用da.Walk()]列出所有的要素类：

import arcpy
import os
walk = arcpy.da.Walk("C:/Project", datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
"""
da.Walk的返回值为三个元组，包括：工作空间、目录名称和文件名称。
dirpath 是字符串形式的工作空间路径。
dirnames 是子目录的名称列表和 dirpath 中的其他工作空间。
filenames 是 dirpath 中的非工作空间内容的名称列表。
"""
    for file in filenames:        
        print(os.path.join(dirpath, file))

结果打印如下：

C:/Project\boundaries.shp
C:\Project\City.gdb\bus_routes
C:\Project\City.gdb\bus_stops
C:\Project\Planning\Education.gdb\schools
C:\Project\Planning\Education.gdb\units
C:\Project\Shapefiles\facilities.shp
C:\Project\Shapefiles\parks.shp
C:\Project\Shapefiles\roads.shp
C:\Project\Water.gdb\streams
C:\Project\Water.gdb\monitoring
C:\Project\Water.gdb\watersheds

对于列出文件，Arcpy中针对某一类的数据类型进行分类编写函数，其实不必太纠结用哪个，甚至你熟悉glob()也可以用它来列出数据路径，同样支持文件匹配。ps：glob()不支持gdb数据库内的内容，但是高版本的geopandas支持。

在列出的文件中处理数据

在GIS工作流中，我们通过列出csv表格文件，通常只是第一步，还会后续操作，例如将每个csv表格中的数据复制到地理数据库表或使用每个csv表格中的数据创建新要素。

ArcPy 使用 Python 列表类型作为其全部列表函数结果的返回类型，因为列表支持简单数据访问所需的灵活性和多种数据类型。for 循环非常适用于处理列表，因为使用它可以一次一个项目的方式浏览列表。for 循环可遍历表中的每一个项目。下面的示例是使用 for 循环遍历前一个示例中生成的列表：

import arcpy
import os
arcpy.env.workspace = "C:/Transportation"
outgdb = "C:/Transportation/City.gdb"
fcs = arcpy.ListFeatureClasses()
for fc in fcs:    
    desc = arcpy.da.Describe(fc)    
    outfc = os.path.join(outgdb, desc["baseName"])    
    arcpy.CopyFeatures_management(fc, outfc)

如果存在子文件夹：

for dirpath, dirnames, filenames in walk:    
    for file in filenames:        
        infc = os.path.join(dirpath, file)        
    desc = arcpy.da.Describe(infc)        
    outfc = os.path.join(outgdb, desc["baseName"])        
    arcpy.CopyFeatures_management(infc, outfc)

使用游标（cursor）访问数据

游标是一个数据库用于访问表（table）中的一组记录或者操作此记录，表中的记录也称作行（rows）。

游标有三种类型：搜索，插入和更新。分别创建于arcpy.da模块的arcpy.da.SearchCursor, arcpy.da.InsertCursor, arcpy.da.UpdateCursor 三个类，其返回一个Python的对象，三种游标都可以在表和要素或要素图层上工作。

除了作为数据访问模块一部分的游标类之外，ArcPy还包含游标函数 arcpy.SearchCursor() 、 arcpy.UpdateCursor() 和 arcpy.InsertCursor() 。数据访问模块中的游标类提供了比这些游标函数更好的性能。

image-20230810120015202

1) 数据锁

使用for来遍历游标对象之后需要关闭或重置游标，否则查询的数据会被锁定，同时在arcgis中操作表中的数据也会被锁定。

❗️提醒：使用游标的时候建议通过with语句可以保证运行代码之后关闭锁，释放数据。

import arcpy
fc = "C:/Data/study.gdb/roads"
with arcpy.da.SearchCursor(fc, "STREET_NAME") as cursor:    
   for row in cursor:        
       print("Street name = {0}".format(row[0]))

关于数据锁的扩展：

插入和更新游标支持 ArcGIS 应用程序设置的表格锁。锁可以防止多个进程同时更改同一个表。锁两种类型：共享和排他，如下所示：

只要访问表或数据集就会应用共享锁。同一表中可以存在多个共享锁，但存在共享锁时，将不允许存在排他锁。显示要素类和预览表是应用共享锁的示例。
对表或要素类进行更改时，将应用排他锁。 在 ArcGIS 中应用排他锁的示例包括：在地图中编辑和保存要素类；更改表的方案；或者在 Python IDE 中在要素类上使用插入游标。

如果数据集上存在排他锁，则无法为表或要素类创建更新和插入游标。UpdateCursor 或 InsertCursor 函数会因数据集上存在排他锁而失败。如果这些函数成功地创建了游标，它们将在数据集上应用排他锁，从而使两个脚本无法在同一数据集上创建更新和插入游标。

释放排他锁的方法：

游标支持 with 语句以重置迭代并帮助移除锁。
但是，为了防止锁定所有内容，应考虑使用 del 语句：del cursor来删除对象或将游标包含在函数中以使游标对象位于作用范围之外。
调用游标上的 reset() 方法也可以释放：cursor.reset()

编辑会话将在其会话期间对数据应用共享锁。保存编辑内容时将应用排他锁。已经存在排他锁时，数据集是不可编辑的。

2）insertRow

插入游标用于创建行并插入它们。创建游标后，insertRow 方法用于插入一组值，这些值会组成新行。表中任何不包含在游标中的字段都将被分配字段的默认值。

import arcpy

# Create insert cursor for table
cursor = arcpy.da.InsertCursor("c:/base/data.gdb/roads_lut", 
                               ["roadID", "distance"])

# Create 25 new rows. Set the initial row ID and distance values
for i in range(0,25):
    cursor.insertRow([i, 100])

3）updateRow

updateRow 方法用于对更新游标当前所在位置的行进行更新。从游标对象返回行后，可以根据需要对行进行修改，然后调用 updateRow 传入修改后的行。

import arcpy

# Create update cursor for feature class
with arcpy.da.UpdateCursor("c:/base/data.gdb/roads",
                           ["roadtype", "distance"]) as cursor:
    for row in cursor:
        # Update the values in the distance field by multiplying 
        #   the roadtype by 100. Road type is either 1, 2, 3 or 4.
        row[1] = row[0] * 100
        cursor.updateRow(row)

4）deleteRow

deleteRow 方法用于对更新游标当前所在位置的行进行删除。提取行后，可在游标上调用 deleteRow 来删除行。

import arcpy

# Create update cursor for feature class
with arcpy.da.UpdateCursor("c:/base/data.gdb/roads", 
                          ["roadtype"]) as cursor:
    # Delete all rows that have a roads type of 4
    for row in cursor:
        if row[0] == 4:
            cursor.deleteRow()

5）几何字段

在 ArcGIS 中，几何数据类型用于指示表中所存储几何的类型（点、线、面、多点或多面体）。通常（但不总是）为存储的名为 Shape 的几何字段。

令牌也可以替代几何字段名称以作为快捷键。返回几何对象的 SHAPE@ 令牌可用于访问要素类几何字段，而无需提前了解字段名称。

使用搜索光标打印点要素类的 x,y 坐标。

import arcpy

infc = "c:/data/fgdb.gdb/fc"

# Enter a for loop for each feature
with arcpy.da.SearchCursor(infc, ['OID@', 'SHAPE@']) as cursor:
    for row in cursor:
        # Print the current point's object ID and x,y coordinates
        print(f'Feature {row[0]}: {row[1][0].X}, {row[1][0].Y}')

可以使用其他几何令牌访问特定几何信息。访问完整几何往往更加耗时。如果仅需要几何的特定属性，可以使用令牌以提供访问几何属性的快捷方式。例如，SHAPE@XY 将返回一组代表要素质心的 x,y 坐标。

还可以进行游标操作的字段类型：全局标识符字段（uuid），BLOB字段（储存图片的二进制信息）、栅格字段。

在Python中使用SQL表达式

在 ArcGIS 中使用的查询表达式的 SQL 参考

地理处理中最常见的处理步骤之一是使用结构化查询语言（SQL）应用查询。SQL基于属性、运算符和计算定义了一个或多个条件。例如，SQL用于ArcGIS Pro地理处理工具，如选择工具和按属性选择图层工具。

SQL的基本教程：

菜鸟教程-SQL。
在 ArcGIS 中使用的查询表达式的 SQL 参考

在搜索游标中使用SQL查询：

arcpy.da.SearchCursor(in_table, field_names {where_clause}, {spatial_reference}, {fields}, {explode_to_points}, {sql_clause})

实用使用技巧：确定表名是否存在的另一种方法是使用 arcpy.CreateUniqueName() 函数。此函数通过在输入名称后附加数字，在指定的工作区中创建唯一名称。此数字会增加，直到名称唯一为止。例如，如果名称“Clip”已经存在，则 CreateUniqueName() 函数将其更改为“Clip 0”;如果这个名称也存在，函数将名称更改为“Clip 1”，依此类推。此函数只能在工作区内创建唯一的表名。它不适用于字段名称。

举例来说：

import arcpy
arcpy.env.workspace = "C:/Data"
unique_name = arcpy.CreateUniqueName("buffer.shp")
arcpy.Buffer_analysis("roads.shp", unique_name, "100 Meters")

计算字段

可以使用Python进行表的字段操作，通过ArcGIS Pro软件中的计算字段或者ArcPy函数CalculateField()实现。等同于更新游标。

示例1：使用ArcPy进行GIS人口空间分布数据探索

本示例简单演示了通过使用arcpy的几种列出数据的方法查看中国人口数据shp文件的信息，通过游标查询单个shp文件的属性表，探索其中的字段，并进行总人口的计算。本示例由jupyter notebook转换而来，可以点击访问原始notebook：4.2.5-示例1：使用Arcpy进行GIS人口空间分布数据探索.ipynb。本示例的数据文件在第七次人口普查数据文件夹中。

数据来源

本次数据为已处理好的分年龄、分性别的人口普查数据，来源于公众号"立方数据学社"。在文件夹的目录结构如下：

'resource\第七次人口普查数据'文件夹的目录

用代码遍历数据

可以用arcpy.da.Walk查看此数据目录：

import arcpy
import os
arcpy.env.workspace = os.path.join(os.getcwd(), "resource/第七次人口普查数据")

# 此目录为多层结构，需要递归遍历
walk = arcpy.da.Walk(arcpy.env.workspace, datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
    for file in filenames:
        # print(os.path.join(dirpath, file))
        fc = os.path.join(dirpath, file)
        desc = arcpy.da.Describe(fc)
        print(desc["baseName"], desc["shapeType"])

输出 >>> 
分年龄、性别的人口_区县等级 Polygon
分年龄、性别的人口_地级市等级 Polygon
分年龄、性别的人口_省份等级 Polygon

将shp导入数据库

先创建一个文件数据库：

outgdb = "census.gdb"
arcpy.CreateFileGDB_management(arcpy.env.workspace, outgdb)

创建成功：

将所有的多边形(Polygon)类型的要素集都复制到这个文件数据库中：

walk = arcpy.da.Walk(arcpy.env.workspace, datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
    for file in filenames:
        fc = os.path.join(dirpath, file)
        desc = arcpy.da.Describe(fc)
        if desc["shapeType"] == "Polygon":
            outfc = os.path.join(outgdb, desc["baseName"].replace("、", "_")) # 不能有特殊字符："、"
            arcpy.CopyFeatures_management(fc, outfc)

查看结果，我们已经将shp文件导入到了数据库中：

另一种数据导入方式： ArcGIS Pro中也提供了一种只能得导入方式叫批量导入数据 (智能)，可以将 KML、KMZ、shapefile、Excel 工作表、表格文本文件、GeoJSON等文件导入地理数据库，支持文件夹和子文件夹导入，也支持文件过滤。

查看数据库中的要素

使用ListFeatureClasses()查看数据库中的要素类信息：

arcpy.env.workspace = os.path.join(os.getcwd(), "resource/第七次人口普查数据", "census.gdb")
fc = arcpy.ListFeatureClasses()
print("共有{}个shp文件".format(len(fc)) )
>>> 共有3个shp文件

print(f"第一个要素文件的名字是: {fc[0]}")
>>> 第一个要素文件的名字是: 分年龄_性别的人口_区县等级

# 查看第一个要素类的空间参考
sr = arcpy.Describe(fc[0]).spatialReference.name
sr
>>> 'GCS_WGS_1984'

# 查看第一个要素文件的字段
fields = arcpy.ListFields(fc[0])
col = [field.name for field in fields]
col
>>>
['OBJECTID',
 'Shape',
 '地名',
 '区划码',
 '县级',
 '县级码',
 '县级类',
 '地级',
 '地级码',
 '地级类',
 '省级',
 '省级码',
 '省级类',
 'F0_男',
  ····省略·····
 'F85及以上女',
 'Shape_Length',
 'Shape_Area']

下面我们通过游标查询一下分省份的人口数据

# 重复一下上面的步骤，过滤出分省份的要素
fc = arcpy.ListFeatureClasses("*省份*")
fc
>>>
['分年龄_性别的人口_省份等级']

查看其字段

fields = arcpy.ListFields(fc[0])
# 我们储存一下字段名
field_names = [field.name for field in fields]

我们通过查询游标来查看一下是不是34个省份的数据都在

num = 0
with arcpy.da.SearchCursor(fc[0], "省") as cursor:
    for row in cursor:
        print(row[0])
        num += 1
print(f"一共有{num}个省份的数据")
>>>
北京市
天津市
河北省
山西省
内蒙古自治区
辽宁省
吉林省
黑龙江省
上海市
浙江省
安徽省
福建省
江西省
山东省
河南省
湖北省
湖南省
广东省
广西壮族自治区
海南省
重庆市
四川省
贵州省
云南省
西藏自治区
陕西省
甘肃省
青海省
宁夏回族自治区
新疆维吾尔自治区
江苏省
台湾省
香港特别行政区
澳门特别行政区
一共有34个省份的数据

我们再通过总人口数来验证一下数据是否正确

# 首先要素表中并未统计总人口数，我们需要自己添加一个字段
arcpy.AddField_management(fc[0], "total", "LONG")

# 重新获取一下字段 查看是否添加成功
fields = arcpy.ListFields(fc[0])
field_names = [field.name for field in fields]
field_names
>>>
['OBJECTID',
 'Shape',
 '省',
 '省级码',
 '省类型',
 'F0_男',
 'F0_女',
 'F1_4_男',
 111111111
 'F80_84_女',
 'F85及以上男',
 'F85及以上女',
 'Shape_Length',
 'Shape_Area',
 'total']

进行各省总人口的计算

# 我们需要传入以下分年龄的人口字段进行游标查询
sum_fields = field_names[5:-3] + field_names[-1:]
sum_fields
>>>
['F0_男',
 'F0_女',
 'F1_4_男',
 ··········
 'F80_84_男',
 'F80_84_女',
 'F85及以上男',
 'F85及以上女',
 'total']

通过更新游标计算总人口数

total_census = 0
with arcpy.da.UpdateCursor(fc[0], sum_fields) as cursor:
    for row in cursor:
        row[-1] = sum(row[:-1])
        total_census += row[-1]
        # cursor.updateRow(row) # 先注释掉，测试代码成功后再运行
print(total_census)
>>>
1409778724

可以看到14亿人口的数据是正确的，我们可以将其写入到一个"total"的字段中了。

# 其次，通过更新游标计算总人口数
total_census = 0
with arcpy.da.UpdateCursor(fc[0], sum_fields) as cursor:
    for row in cursor:
        row[-1] = sum(row[:-1])
        total_census += row[-1]
        cursor.updateRow(row)