文章/答案/技术大牛

发布

社区首页 >专栏 >六、处理几何数据【ArcGIS Python系列】

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

renhai

发布于 2023-11-24 08:24:17

56200

代码可运行

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

运行总次数：0

代码可运行

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])
array2 = arcpy.Array([point2a, point2b])

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

会输出：

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文件，把人口数据在空间上呈现。通常，这是做研究的基础工作，方便了解我们数据在空间上是如何分布的，比如横向对比每个省份之间的总人口差异有哪些，每个省份年龄构成差异有哪些，年龄结构和经济的关系，你可以纵向对比多次人口普查在空间上的差异，这些都是可以进行深挖的方向。

为此我们准备的数据有：

中国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 分年龄、性别的人口(全部数据).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")

查看此要素，为无人口数据的空地图：

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)