ArcPy高级开发教程—要素操作

今天，介绍使用ArcPy对矢量数据的要素feature进行操作方法。

目录

1. 空间分析相关

建立缓冲区

计算要素面积和周长

按属性选择

2. 要素遍历

3. 导出每个要素

4. 生成shpfile文件

---

1. 空间分析相关

• 建立缓冲区

# _*_ coding: UTF-8 _*_
import arcpy, os
def buildBufffer(shpPath,workspace):
    buffer_name = os.path.join(workspace, 'buffer.shp')
    distance = '200 meters'
    if not os.path.exists(buffer_name):
        arcpy.Buffer_analysis(shpPath, buffer_name, distance, "FULL", "ROUND", "NONE")
    return buffer_name

上述函数能够在workspace目录下生成名为buffer的缓冲区文件。shpPath是shpfile的路径，缓冲半径为200米，在哪一侧缓冲、末端类型和是否融合都可以在Buffer_analysis函数中更改。

• 计算要素面积和周长

def caculateAreaLength(inputShp):
    # 计算面积
    arcpy.AddGeometryAttributes_management(inputShp,
                                           ['AREA_GEODESIC'], '', 'SQUARE_METERS', '')
    # 计算周长
    arcpy.AddGeometryAttributes_management(inputShp,
                                           ['PERIMETER_LENGTH_GEODESIC'], 'METERS','', '')

上述函数为输入的shp的元素自动计算面积和周长，注意输入的shpfile必须为投影坐标系，如果原始shpfile为地理坐标系（例如WGS84），那么可以在AddGeometryAttributes_management函数前添加下面代码，帮助输出更改为投影坐标。

arcpy.env.outputCoordinateSystem=arcpy.SpatialReference('WGS 1984 UTM Zone 49N')

此外，计算图斑面积和周长，据笔者所知，还有以下方法：

# 通过字段计算器
arcpy.CalculateField_management(shpPath, 'area', "!shape.geodesicArea@SQUAREMETERS!", 'PYTHON_9.3')

# 通过空间统计工具
arcpy.CalculateAreas_stats(input,output)

• 按属性选择

这是一个十分常用的功能，比如我已经计算出了每个图斑的面积，现在我要删除面积小于100的，保留面积大于等于100的图斑，参考代码如下：

def SelectArea(inputShp):
    input_lyr = "lyr1"
    arcpy.MakeFeatureLayer_management(inputShp, input_lyr)
    arcpy.SelectLayerByAttribute_management(input_lyr, "NEW_SELECTION", '"AREA_GEO" < 100')
    if int(arcpy.GetCount_management(input_lyr).getOutput(0)) > 0:
        arcpy.DeleteFeatures_management(input_lyr)
    arcpy.Delete_management(input_lyr)

注意，这个操作是在原始shp上进行的，会直接删除不符合条件的要素。SelectLayerByAttribute_management函数中AREA_GEO就是原始shp中代表面积的字段。这个操作就类似于你在ArcGIS里面按属性选择，说完这个，好奇宝宝会问，那么按位置选择如何实现，简要代码如下

arcpy.SelectLayerByLocation_management(layer1, 'intersect',layer2)

这句代码的作用就是选择layer1中与layer2相交的要素，此外，我们还可以统计选择的数量，例如在上句代码后加入

matchcount = int(arcpy.GetCount_management(layer1)[0])

这里的matchcount就是layer1中被选中的要素个数

最后，如果你要对按属性选择或者按位置选择后的结果导出，命令如下

arcpy.CopyFeatures_management(layer1,outputshp)

先选择，选择后就可以导出结果了。

2. 要素遍历

遍历每一个要素，用到的主要是SearchCursor和UpdateCursor，顾名思义，SearchCursor主要用于你要查询某个要素的属性值，UpdateCursor用于更新某要素的属性值。

fields = ['FID', 'field1','field2']
with arcpy.da.SearchCursor(inputShp, fields) as cursor:
    for r in cursor:
        id = r[0]

上面代码中通过一个列表fields传入你要查询的字段，inputShp是输入的shp，通过游标cursor进行遍历，遍历的结果仍然以列表的形式返回，返回结果的顺序与fields中查询字段的顺序是一致的。

下面举个复杂点的例子，inputShp包含FID, AREA和Tag三个字段，AREA存储了每个要素的面积，现在我要把AREA>100的要素的Tag的值设为1，代码如下：

id_list = []
fields = ['FID','AREA']
with arcpy.da.SearchCursor(inputShp, fields) as cursor:
    for r in cursor:
        if int(r[1]) > 100:   
            id_list.append(r[0])
fields = ['FID', 'Tag']
with arcpy.da.UpdateCursor(inputShp, fields) as cursor:
    for r in cursor:
        if r[0] in id_list:
            r[1] = 1
        cursor.updateRow(r)

另外，一个很有用且相对复杂的，是ArcPy提供了一个叫做令牌'SHAPE@'的东西，帮助在遍历要素时直接获取要素的几何属性、空间参考等信息，使用时直接放在fields中就行。示例代码如下

# 面要素，获取每个要素的几何范围
with arcpy.da.UpdateCursor(inputShp, ['SHAPE@']) as cursor:
    for row in cursor:
        extent = row[0].extent
​
# 点要素，获取每个要素的XY坐标
with arcpy.da.SearchCursor(inputShp, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        x, y = row[0]

'SHAPE@'获取得到的是Geometry几何对象，"SHAPE@XY"得到的是要素的质心坐标，除了这两个，还有一些令牌可以直接使用，如下图：

3. 导出每个要素

这是一个经常有很多人问的问题，如何把一个shpfile按照ID把每个要素分别导出成shp,了解了要素遍历，这个就很简单了

• 方式1：

for row in arcpy.da.SearchCursor(inputShp,["FID"]):
    tempShpName=str(row[0])+'.shp'
    where_clause = '{fieldName}={data}' \
        .format(fieldName='FID', data=row[0])
    arcpy.Select_analysis(inputShp, tempShpName, where_clause)

• 方式2

with arcpy.da.SearchCursor(shp, ["SHAPE@","FID"]) as cursor:
    for row in cursor:
        out_name = str(row[1]) + '.shp'  
        arcpy.FeatureClassToFeatureClass_conversion(row[0], out_path, out_name)

4. 生成shpfile文件

经常我们会从txt或者excel等文件中读取点坐标，生成矢量文件。这里给定两个坐标点(x1,y1),(x2,y2)

• 生成包含这两个点的点shp文件：

pointGeoList = []
spatial_reference = arcpy.SpatialReference(4326)
pointList=[arcpy.Point(x1, y1),arcpy.Point(x2, y2)]
for pt in pointList:
    point_geo = arcpy.PointGeometry(pt, spatial_reference)
    pointGeoList.append(point_geo)
arcpy.CopyFeatures_management(pointGeoList, 'points.shp')

上述代码将两个点保存在points.shp中，且空间参考为WGS84

• 生成面文件

如果我们想以上面的两个点作为矩形的左上角和右下角坐标，生成一个面矢量，则

polygonList = []
spatial_reference = arcpy.SpatialReference(4326)
pointArray = []
point1 = arcpy.Point(x1, y1)
pointArray.append(point1)
point2 = arcpy.Point(x2, y1)
pointArray.append(point2)
point3 = arcpy.Point(x2, y2)
pointArray.append(point3)
point4 = arcpy.Point(x1, y2)
pointArray.append(point4)
polygonList.append(arcpy.Polygon(arcpy.Array([pointArray]), spatial_reference))
arcpy.CopyFeatures_management(polygonList, 'polygon.shp')

创建面的关键就是面也是由点组成的，将一个面中的点按顺序存放在list中，组成一个Polygon，再将多个Polygon放在list中组成包含多个面的矢量。

更多信息可以参考ArcPy的官方文档

原创不易，如果有问题我没有在博客上及时回复，可以关注公众号 TechGIS 交流