六、处理几何数据【ArcGIS Python系列】
1.了解几何对象
要素类中的每个要素都由一个或多个顶点组成，这些顶点定义了点、多段线或多边形要素。在点要素类的情况下，每个点要素由单个顶点组成。多段线和多边形要素由多个顶点组成。每个顶点是由一对 x、y 坐标定义的位置。该图说明了点、多段线和多边形如何在笛卡尔坐标空间中由顶点定义。
使用几何体对象可以将要素写入要素类，我们可以从坐标值表创建要素。几何对象可用于地理处理操作，能在内存中直接创建并在地理处理工具中使用，无需创建临时要素类来保存。
ArcGIS 中主要有五中几何对象：Point、MultiPoint 、 PointGeometry 、 Polygon 和 Polyline 。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 还使用借助两个类来帮助构建几何图形： Array 和 Point 。
扩展：
对比Shapely包：Shapely中有
Point、LineString、Polygon、MultiPoint、MultiLineString、MultiPolygon、GeometryCollection，也支持从 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 表示获取第一部分。对于只有一个部分的要素类，第一个部分也是唯一的部分。图示就是一个包含多个多边形的多部分集合图形。
多边形要素类
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))`
会输出：
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
会输出：
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])
结果如图（显示的顶点是为了强调）：
你也可以从硬盘中读取坐标点：
filename = "coordinates.txt"
with open(filename, "r") as file:
coordinates = [line.strip().split(",") for line in file]
# 打印坐标列表
print(coordinates)
总的来说，Arcpy 中的几何对象可以提高代码的效率，大部分几何对象函数创建返回的对象也是几何对象，避免了创建临时要素类和使用光标读取所有要素的步骤。
示例：从 excel 表格制作分年龄的人口普查要素文件
此示例演示了如何通过表格数据制作分年龄、性别的人口_省份等级.shp文件，把人口数据在空间上呈现。通常，这是做研究的基础工作，方便了解我们数据在空间上是如何分布的，比如横向对比每个省份之间的总人口差异有哪些，每个省份年龄构成差异有哪些，年龄结构和经济的关系，你可以纵向对比多次人口普查在空间上的差异，这些都是可以进行深挖的方向。
为此我们准备的数据有：
- 中国 34 个省市区的空地图：中国各省份地图.shp
- 分年龄的人口统计数据：中国第七次人口普查-分年龄_性别的人口数据.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.处理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
命名之后的数据：
删除不需要的第一行数据：
df = df.drop(0, axis=0).reset_index(drop=True)
df.head()
3.读取省份地图
我们用geopandas读取地图数据，然后用pandas读取人口数据，然后通过merge方法进行匹配，最后用geopandas导出为shp文件。
读取地图用geopandas：
import geopandas as gpd
gdf = gpd.read_file("./resource/data2/中国分省份空地图/中国各省份地图.shp", encoding="utf-8")
gdf.head()
4.进行数据合并
gdf_new = pd.merge(gdf, df, left_on="省", right_on="Province", how="left").drop("Province", axis=1)
gdf_new.head()
合并数据如下：
查看尾部数据，可以看到，我们的数据集中包含了中国的 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")
查看此要素，为无人口数据的空地图：
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)
最终我们通过两种方法生成了 shp 文件。
文章索引
【ArcGIS Python 系列】系列笔记为学习 ArcGIS Pro 和
Arcpy过程中的总结，记下来方便回看，最新版本会优先发布在我的博客和GITHUB。
【ArcGIS Python 系列】教程部分：
- 一、Arcpy 介绍和安装【ArcGIS Python 系列】
- 二、ArcGIS Pro 和 ArcMap 的区别【ArcGIS Python 系列】
- 三、Arcpy 基础【ArcGIS Python 系列】
- 四、探索空间数据【ArcGIS Python 系列】
- 五、处理地理数据异常【ArcGIS Python 系列】
- 六、处理几何数据【ArcGIS Python 系列】
- 七、处理栅格数据【ArcGIS Python 系列】
- 八、制图模块【ArcGIS Python 系列】
- 九、自定义工具箱【ArcGIS Python 系列】
- 十、ArcGIS_Pro 常见问题【ArcGIS Python 系列】
- 利用 ArcGIS Python 批量处理地理数据的坐标系
- 使用 ArcGIS Python 检测洪水影像的区域
- 利用 ArcGIS_Python 制作考虑路况的交通等时圈
- 利用 ArcGIS Pro 制作弧线 OD 图【ArcGIS Python 系列】
- 使用 ArcGIS Pro 对卫星图进行建筑轮廓识别和车辆检测
- ArcGIS_Pro 官方课程整理
- 持续更新…
【ArcGIS Python 系列】jupyter notebook：
- 4.2.3-arcpy 基础(代码练习).ipynb
- 4.2.4-探索空间数据(代码练习).ipynb
- 4.2.5-示例 1：使用 Arcpy 进行 GIS 人口空间分布数据探索.ipynb
- 4.2.7-处理几何数据代码练习和示例 2.ipynb
- 4.2.8-栅格数据(代码练习).ipynb
- 4.2.9-制图模块.ipynb
- 4.2.12-实操 1-利用 Python 批量处理地理数据的坐标系.ipynb
- 4.2.13-实操 2-使用 Python 对图像中的洪水进行分类.ipynb
- 4.2.14-实操 3-制作考虑路况的交通等时圈.ipynb
如果你觉得本系列文章有用，欢迎关注博客，点赞 👍 和收藏，也欢迎在评论区讨论，也欢迎访问我的爱发电支持我，或者对此文章进行赞赏。
其他平台账号：
