六、处理几何数据【ArcGIS Python系列】

1.了解几何对象

要素类中的每个要素都由一个或多个顶点组成,这些顶点定义了点、多段线或多边形要素。在点要素类的情况下,每个点要素由单个顶点组成。多段线和多边形要素由多个顶点组成。每个顶点是由一对 x、y 坐标定义的位置。该图说明了点、多段线和多边形如何在笛卡尔坐标空间中由顶点定义。

使用几何体对象可以将要素写入要素类,我们可以从坐标值表创建要素。几何对象可用于地理处理操作,能在内存中直接创建并在地理处理工具中使用,无需创建临时要素类来保存。

image-20230810145329873

ArcGIS 中主要有五中几何对象:PointMultiPointPointGeometryPolygonPolyline 。ArcPy 通常用 arcpy.Geometry 类创建几何对象。 Geometry 类的一般语法如下:

arcpy.Geometry(geometry, inputs, {spatial_reference}, {has_z}, {has_m})
参数 说明 数据类型
geometry 几何类型:点、面、折线或多点。 String
inputs 用于创建对象的坐标。数据类型可以是 Point 或 Array 对象。 Object
spatial_reference 新几何的空间参考。(默认值为 None) SpatialReference
has_z Z 状态:如果启用 Z,则为几何的 True,如果未启用,则为 False。(默认值为 False) Boolean
has_m M 状态:如果启用 M,则为几何的 True,如果未启用,则为 False。(默认值为 False) Boolean

此外,ArcPy 还使用借助两个类来帮助构建几何图形: ArrayPoint

扩展:

对比Shapely包:Shapely中有PointLineStringPolygonMultiPointMultiLineStringMultiPolygonGeometryCollection,也支持从 numpy 的 array 对象创建几何对象。

Point 和 PointGeometry

Point 无空间参考信息,通常是一对点坐标,PointGeometry有空参考信息,是一个几何对象。下面的代码演示了 Point 对象如何使用 PointGeometry 类构造几何对象:

point = arcpy.Point(4.900160, 52.378424)
pointgeo = arcpy.PointGeometry(point, 4326) # 4326等同于GCS_WGS_1984

Polyline

多段线和多边形要素由多个顶点组成,并使用两个或多个 Point 对象构造。为便于处理多个 Point 对象,ArcPy 使用 Array 类。此类专门为构造多段线和多边形几何对象而创建。以下示例显示如何使用两个 Point 对象创建一个 Polyline 对象:

point1 = arcpy.Point(0, 0)
point2 = arcpy.Point(100, 100)
array = arcpy.Array([point1, point2])
polyline = arcpy.Polyline(array)

print(polyline.length)
>>> 141.4213562373095

可以创建多个几何体对象,并使用几何体对象的方法直接比较它们。例如,下面的代码创建了两个 Polyline 对象,并确定它们是否相互交叉:

import arcpy
point1a = arcpy.Point(0,0)
point1b = arcpy.Point(100, 100)
point2a = arcpy.Point(100, 0)
point2b = arcpy.Point(0, 100)
array1 = arcpy.Array([point1a, point1b])
arrayhttps://blog.renhai.onlineoint2b])

polyline1 = arcpy.Polyline(array1)
polyline2 = arcpy.Polyline(array2)

print(polyline1.crosses(polyline2))
>>> True # 两条折线在位置(50,50)相交。

Polygon

创建 Polygon 对象的语法和Polyline类似,但有意义的面要素至少需要三个 Point 对象。

2.读取几何对象属性

我们已经理解了几何对象,现在可以通过搜索游标来访问要素类的几何对象。再此之前我们要先了解**几何令牌**:

几何令牌可以作为快捷方式来替代访问完整几何对象。附加几何令牌可用于访问特定几何信息。访问完整几何往往更加耗时。如果只需要几何的某些特定属性,可使用令牌来提供快捷方式从而访问几何属性。例如,SHAPE@XY 会返回一组代表要素质心的 x,y 坐标。常用的几何令牌有:

令牌 说明
SHAPE@ 要素的几何对象。
SHAPE@XY 一组要素的质心 x,y 坐标。
SHAPE@Z 要素的双精度 z 坐标。
SHAPE@AREA 要素的双精度面积。
SHAPE@LENGTH 要素的双精度长度。

我们结合游标和几何令牌探索Point要素的坐标:

fc = "./resource/第七次人口普查数据/中国各城市中心/中国各城市中心.shp"
with arcpy.da.SearchCursor(fc, ["SHAPE@XY", "city"]) as cursor: # 返回几何对象的特定值,此处为xy坐标,元祖格式
    for row in cursor:
        x, y = row[0] # 解开元祖
        city = row[1]
        print(f"{city}中心点坐标: {x}, {y}") # 打印点坐标

>>> 北京市中心点坐标: 116.39904202900004, 39.90358526800003
天津市中心点坐标: 117.18386143100008, 39.124546990000056
石家庄市中心点坐标: 114.49643611600004, 38.04481885200005
......
......

我们结合游标和几何令牌探索Polygon要素的坐标:

import arcpy
arcpy.env.workspace = "C:/Data/Demo.gdb"
fc = "pipes"
with arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"]) as cursor: # 使用字符串 "OID@" 引用唯一标识符字段。
	for row in cursor:
		print("Feature {0}: ".format(row[0]))
			for point in row[1].getPart(0):
				print("{0}, {1}".format(point.X, point.Y))

getPart()方法用于获取几何形状的各个部分,这里使用索引 0 表示获取第一部分。对于只有一个部分的要素类,第一个部分也是唯一的部分。图示就是一个包含多个多边形的多部分集合图形。

image-20230810162922634

多边形要素类

import arcpy
arcpy.env.workspace = "C:/Data/Demo.gdb"
fc = "parcels"
	with arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"]) as cursor:
		for row in cursor:
			print("Feature {0}: ".format(row[0]))
				for point in row[1].getPart(0):
					print(f"{point.X:.2f}, {point.Y:.2f}") # `{point.X:.2f}`表示将`point.X`的值格式化为浮点数,并保留两位小数。等同于`print("{0:.2f}, {1:.2f}".format(point.X, point.Y))`

image-20230810162738533

会输出:

Feature 1:
427837.21, 5450864.65
427837.80, 5450898.34
427842.31, 5450898.89
427846.84, 5450899.18
427851.39, 5450899.24
427859.98, 5450899.11
427859.42, 5450867.42
427837.21, 5450864.65

示例三:

多部分几何图形

import arcpy
arcpy.env.workspace = "C:/Data/Demo.gdb"
fc = "multipart_features"
with arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"]) as cursor:
	for row in cursor:
    print("Feature {0}: ".format(row[0]))
    partnum = 1
    for part in row[1]:
      print("Part {0}:".format(partnum))
      for point in part:
        print("{0}, {1}".format(point.X, point.Y))
      partnum += 1

会输出:

image-20230810163433024

3.写入几何数据

如果你熟悉对 geopandas 更加熟悉,推荐使用 geopandas 和 shapely。实际使用的时候更多的还是从 csv、json 构建几何对象,还是直接读取 shp、geojson 等文件,这些库处理起来都会比 arcpy 顺手很多。

有两种方法写入几何数据:使用arcpy.CopyFeatures_management()将几何对象复制到要素类和使用arcpy.da.InsertCursor()插入游标。

arcpy.CopyFeatures_management()可以让代码更加简洁,但是也有缺点:

  • 复制几何图元时不能创建或更新特征的属性。因此,如果必须同时创建新要素和属性值,则必须使用游标。

  • 创建许多要素(尤其是由许多顶点组成的要素)可能会降低性能,因为必须同时将所有几何图形对象加载到内存中,才能将它们复制到要素类。使用游标时,可以在游标对象上的每次迭代中创建每个新特征,这样可以在处理许多特征时获得更好的性能。

下面我们从 x,y 坐标对列表创建新多边形要素的完整实例来看两种方法的区别,首先演示使用arcpy.CopyFeatures_management()

import arcpy
point = arcpy.Point()
array = arcpy.Array()
coordinates = [[3116036.11, 10071403.50], [3115768.36, 10071482.07], [3115847.82, 10071747.21], [3116114.23, 10071667.17]] # 列表形式的坐标

# 接下来,代码遍历坐标对列表,并在每次迭代中创建一个新的 Point 对象。
for coord in coordinates:
  point.X = coord[0]
  point.Y = coord[1]
  array.add(point)

# 将array对象构造成多边形
polygon = arcpy.Polygon(array, 2277)

# 将内存中的多边形创建为数据库的新要素
arcpy.CopyFeatures_management(polygon, fc)

使用arcpy.da.InsertCursor()插入游标

import arcpy
point = arcpy.Point()
array = arcpy.Array()
coordinates = [[3116036.11, 10071403.50], [3115768.36, 10071482.07], [3115847.82, 10071747.21], [3116114.23, 10071667.17]] # 列表形式的坐标

# 接下来,代码遍历坐标对列表,并在每次迭代中创建一个新的 Point 对象。
for coord in coordinates:
  point.X = coord[0]
  point.Y = coord[1]
  array.add(point)

# 将array对象构造成多边形
polygon = arcpy.Polygon(array) # , 2277可省略

# 将内存中的多边形创建为数据库的新要素
fgdb = "C:/Data/Demo.gdb"
fc = "newpoly"
arcpy.env.workspace = fgdb
# 先创建一个空要素
arcpy.CreateFeatureclass_management(fgdb, fc, "Polygon",                                     "", "", "", 2277)
# 将多边形插入
with arcpy.da.InsertCursor(fc, "SHAPE@") as cursor:
  cursor.insertRow([polygon])

结果如图(显示的顶点是为了强调):

image-20230810170736165

你也可以从硬盘中读取坐标点:

filename = "coordinates.txt"

with open(filename, "r") as file:
    coordinates = [line.strip().split(",") for line in file]

# 打印坐标列表
print(coordinates)

总的来说,Arcpy 中的几何对象可以提高代码的效率,大部分几何对象函数创建返回的对象也是几何对象,避免了创建临时要素类和使用光标读取所有要素的步骤。


示例:从 excel 表格制作分年龄的人口普查要素文件

代码文件在4.2.7-处理几何数据代码练习和示例 2.ipynb

此示例演示了如何通过表格数据制作分年龄、性别的人口_省份等级.shp文件,把人口数据在空间上呈现。通常,这是做研究的基础工作,方便了解我们数据在空间上是如何分布的,比如横向对比每个省份之间的总人口差异有哪些,每个省份年龄构成差异有哪些,年龄结构和经济的关系,你可以纵向对比多次人口普查在空间上的差异,这些都是可以进行深挖的方向。

为此我们准备的数据有:

  1. 中国 34 个省市区的空地图:中国各省份地图.shp
  2. 分年龄的人口统计数据:中国第七次人口普查-分年龄_性别的人口数据.xlsx

方法一:通过 Python 的 pandas 和 geopandas 处理

如果你需要在 arcpy 的环境下安装库,推荐用 conda 安装环境,避免库之间的冲突,出现错误了也能够回滚环境。
geopandas 通过运行 conda install geopandas -c esri 来安装。

1.读取 excel 文件
import pandas as pd
df = pd.read_excel("./resource/data2/中国第七次人口普查-分年龄_性别的人口数据.xlsx")
df.head()


通过对比上述df对象和原始表格,首先发现,表头需要处理,要将合并的单元格拆散,比如年龄0岁拆分成0岁男和0岁女。然后,表格中包含有省级的也有市县一级的数据,我们只需要省级信息,只是表格没有可以供筛选的字段,我们可以下一步通过 pandas 合并表格的时候直接扔掉不匹配的行。最后需要注意的是,表格内无港澳台人口信息,因为第七次人口普查就没有统计,但是地图必须有港澳台!!!:

表2 分年龄、性别的人口(全部数据).xlsx

2.处理dataframe

我们初步处理此dataframe,首先删除空值,去除第一行:

df = df.dropna(axis=0, how='all')
df = df.drop(1, axis=0).reset_index(drop=True)
df.head()

重命名列名:

# 此处直接指定列名
fields = [
    "Province",
    "M0",
    "F0",
    "M1_4",
    "F1_4",
    "M5_9",
    "F5_9",
    "M10_14",
    "F10_14",
    "M15_19",
    "F15_19",
    "M20_24",
    "F20_24",
    "M25_29",
    "F25_29",
    "M30_34",
    "F30_34",
    "M35_39",
    "F35_39",
    "M40_44",
    "F40_44",
    "M45_49",
    "F45_49",
    "M50_54",
    "F50_54",
    "M55_59",
    "F55_59",
    "M60_64",
    "F60_64",
    "M65_69",
    "F65_69",
    "M70_74",
    "F70_74",
    "M75_79",
    "F75_79",
    "M80_84",
    "F80_84",
    "M85+",
    "F85+"
]
df.columns = fields

命名之后的数据:

image-20230813115106174

删除不需要的第一行数据:

df = df.drop(0, axis=0).reset_index(drop=True)
df.head()

image-20230813115133806


3.读取省份地图

我们用geopandas读取地图数据,然后用pandas读取人口数据,然后通过merge方法进行匹配,最后用geopandas导出为shp文件。

读取地图用geopandas

import geopandas as gpd
gdf = gpd.read_file("./resource/data2/中国分省份空地图/中国各省份地图.shp", encoding="utf-8")
gdf.head()

image-20230813120628597


4.进行数据合并
gdf_new = pd.merge(gdf, df, left_on="省", right_on="Province", how="left").drop("Province", axis=1)
gdf_new.head()

合并数据如下:

image-20230813120720612

查看尾部数据,可以看到,我们的数据集中包含了中国的 34 个省级行政区,以及港澳台地区。只不过港澳台地区的数据是空的,因为我们的数据集中没有这些地区的数据。

5.处理数据类型

人口数量字段肯定是数字类型,我们通过astype将字段转化为整数型:

可以先查看数据类型,将人口字段转换为整数型 int:

gdf_new.dtypes
>>>
省             object
省级码            int64
geometry    geometry
M0            object
F0            object
M1_4          object
F1_4          object
M5_9          object
F5_9          object
.....................
M85+          object
F85+          object
dtype: object

将数据列转化为整数

# 选取需要转换的列
cols_to_convert = gdf_new.columns[3:]
# 应用匿名函数进行转换
gdf_new[cols_to_convert] = gdf_new[cols_to_convert].apply(lambda x : x.astype(int) if not pd.isnull else x)
# 再次查看类型
gdf_new.dtypes
>>>
省             object
省级码            int64
geometry    geometry
M0            object
F0            object
M1_4          object
F1_4          object
M5_9          object
F5_9          object
....................
F75_79        object
M80_84        object
F80_84        object
M85+          object
F85+          object
dtype: object

6.保存为 shp 文件:
import os
if not os.path.exists("./resource/data2/output"):
    os.mkdir("./resource/data2/output")
gdf_new.to_file("./resource/data2/output/分年龄、性别的人口_省份等级_方法1.shp", encoding="utf-8", driver="ESRI Shapefile",  engine="pyogrio") # pip安装pyogrio

我本来用的中文列名。然后遇到了编码问题。最后还是老老实实把字段写成英文,就没问题了。


方法二:使用 Arcpy 的游标来管理数据

此方法如果只用 arcpy 的游标更新数据,相对来说没有merge方便。

import arcpy
import pandas as pd
import os

# 继续用这个表格
df.head()
1.创建数据库和要素

我们先创建数据库,然后将数据导入到数据库中,这样就可以避免覆盖原有的数据了。

arcpy.env.workspace = "./resource/data2"
arcpy.env.overwriteOutput = True

# 创建一个数据库来操作 避免覆盖
if os.path.exists("./resource/data2/output.gdb") == False:
    arcpy.management.CreateFileGDB("./resource/data2", "output.gdb")

# 1.选择空的分省份地图要素
fc = "./resource/data2/中国分省份空地图/中国各省份地图.shp"

# 复制到数据库
arcpy.conversion.FeatureClassToFeatureClass(fc, "./resource/data2/output.gdb", "分年龄、性别的人口_省份等级_方法2")

查看此要素,为无人口数据的空地图:

fc

2.进行 arcpy 的字段操作:
# 3.添加人口字段
new_field = df.columns[1:].tolist()

# 使用插入游标添加字段
for field in new_field:
    arcpy.management.AddField("output.gdb/分年龄、性别的人口_省份等级_方法2", field, "LONG")
# 列出字段 发现字段名被arcgis规整化了
gis_field = [field.name for field in arcpy.ListFields("output.gdb/分年龄、性别的人口_省份等级_方法2")]
gis_field = gis_field[2:3] + gis_field[6:]
gis_field
>>>
['省',
 'M0',
 'F0',
 'M1_4',
 'F1_4',
 'M5_9',
 'F5_9',
 'M10_14',
 'F10_14',
 'M15_19',
 'F15_19',
 'M20_24',
 'F20_24',
 'M25_29',
 'F25_29',
 'M30_34',
 'F30_34',
 'M35_39',
 'F35_39',
 'M40_44',
 'F40_44',
 'M45_49',
 'F45_49',
 'M50_54',
 'F50_54',
 'M55_59',
 'F55_59',
 'M60_64',
 'F60_64',
 'M65_69',
 'F65_69',
 'M70_74',
 'F70_74',
 'M75_79',
 'F75_79',
 'M80_84',
 'F80_84',
 'M85_',
 'F85_']
3.使用游标添加数据,通过省份字段匹配:
with arcpy.da.UpdateCursor("output.gdb/分年龄、性别的人口_省份等级_方法2", gis_field
) as cursor:
    for i, row in enumerate(cursor):
        # 通过游标所在行的"省"字段匹配df中"Province"的字段
        filtered_df = df[df['Province'] == row[0]]
        # 判断是否匹配
        if filtered_df.empty:
            # 如果不匹配此dataframe为空
            print("不匹配")
        else:
            print("匹配")

            # 赋值
            for j, field in enumerate(gis_field[1:]):
                # 因为gis把字段规整化了 有两个例外我们替换一下
                if field == "F85_":
                    field = "F85+"
                elif field == "M85_":
                    field = "M85+"

                row[j+1] = filtered_df.loc[:, field].values[0] # 这里的j+1是因为gis_field中去掉了"省"字段

            # 最后执行更新游标
            cursor.updateRow(row)

            # break # 测试用

运行完整之后和原始 df 对比一下:

df.iloc[0, :]
Province        北京市
M0            79388
F0            73342
M1_4         447165
F1_4         416355
M5_9         485161
F5_9         447792
M10_14       335714
F10_14       306590
M15_19       345715
F15_19       287842
M20_24       715999
F20_24       634503
M25_29       995896
F25_29       908792
M30_34      1309189
F30_34      1193840
M35_39      1108136
F35_39      1035049
M40_44       838444
F40_44       763689
M45_49       842337
F45_49       776237
M50_54       842484
F50_54       777307
M55_59       828013
F55_59       799526
M60_64       677025
F60_64       709505
M65_69       568385
F65_69       626286
M70_74       314838
F70_74       353703
M75_79       184665
F75_79       230508
M80_84       152087
F80_84       196699
M85+         124749
F85+         160140
Name: 0, dtype: object
4.更改别名

更改 arcgis 属性表中的别名,方便阅读

alias = [
    "男性0岁人口",
    "女性0岁人口",
    "男性1-4岁人口",
    "女性1-4岁人口",
    "男性5-9岁人口",
    "女性5-9岁人口",
    "男性10-14岁人口",
    "女性10-14岁人口",
    "男性15-19岁人口",
    "女性15-19岁人口",
    "男性20-24岁人口",
    "女性20-24岁人口",
    "男性25-29岁人口",
    "女性25-29岁人口",
    "男性30-34岁人口",
    "女性30-34岁人口",
    "男性35-39岁人口",
    "女性35-39岁人口",
    "男性40-44岁人口",
    "女性40-44岁人口",
    "男性45-49岁人口",
    "女性45-49岁人口",
    "男性50-54岁人口",
    "女性50-54岁人口",
    "男性55-59岁人口",
    "女性55-59岁人口",
    "男性60-64岁人口",
    "女性60-64岁人口",
    "男性65-69岁人口",
    "女性65-69岁人口",
    "男性70-74岁人口",
    "女性70-74岁人口",
    "男性75-79岁人口",
    "女性75-79岁人口",
    "男性80-84岁人口",
    "女性80-84岁人口",
    "男性85岁以上人口",
    "女性85岁以上人口"
]

for i, field in enumerate(gis_field[1:]):
    arcpy.management.AlterField("output.gdb/分年龄、性别的人口_省份等级_方法2", field, "" , alias[i])

更改别名后的属性表:更改别名后的属性表

5.导出为 shp 文件
in_features = "output.gdb/分年龄、性别的人口_省份等级_方法2"
out_features = "output/分年龄、性别的人口_省份等级_方法2.shp"

arcpy.conversion.ExportFeatures(in_features, out_features)

image-20230813122447208

最终我们通过两种方法生成了 shp 文件。


文章索引

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

【ArcGIS Python 系列】教程部分:

【ArcGIS Python 系列】jupyter notebook:

托管在 Github:Urban-Spatial-Data-Analysis-Notebook/4-空间数据分析/4.2-ArcGIS Python 系列 at main · renhai-lab/Urban-Spatial-Data-Analysis-Notebook · GitHub


如果你觉得本系列文章有用,欢迎关注博客,点赞 👍 和收藏,也欢迎在评论区讨论,也欢迎访问我的爱发电支持我,或者对此文章进行赞赏。

donate

其他平台账号:
donate


六、处理几何数据【ArcGIS Python系列】
https://blog.renhai.online/archives/4.2.7-%E5%A4%84%E7%90%86%E5%87%A0%E4%BD%95%E6%95%B0%E6%8D%AE
作者
Renhai
发布于
2023年08月29日
更新于
2024年06月19日
许可协议