arcpy基础篇（4）-处理几何对象

1.处理几何对象

要素类的每一个要素都包含一组用来定义、线和多边形的折点。这些可以通过几何对象(例如Point、Polyline、PointGeometry和MultiPoint来访问折点)来访问，并以点对象数组的形式返回。
如果仅仅需要该几何的某几个属性，可以定义一个”几何短语“来获取几何对象的属性。例如，SHAPE@XY会返回一组表示要素几何中心的x、y坐标值；SHAPE@LENGTH会以双精度类型返回要素的长度；SHAPE@会返回整个几何对象。
下面的代码用于计算一个线要素类中所有线要素的总长度。

import arcpy
fc = "C:/Data/roads.shp"
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@LENGTH"])
length = 0
for row in cursor:
	length += row[0]
print length

2.读取几何

要素类中的每一个要素都包含一组用来定义点、线和多边形的折点。在点要素类中，一个点要素代表一个折点。
折点可以通过搜索游标进行访问，访问点要素时，将会返回一个点对象；访问其他类型的要素时，将返回一个点对象数组。

2.1读取点对象

下面的例子使用搜索游标和for循环遍历了点要素类hospital.shp中的每一个点要素，并通过几何短语获取所有点要素的x、y坐标。脚本程序如下：

import arcpy
fc = "C:/Data/hosptials.shp"
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@XY"])
for row in cursor:
	x,y =row[0]
	print("{0}, {1}",format(x,y))

2.2读取线、多边形

线、多边形对象，访问其中的每一个要素会返回一个点对象。相应处理这些数组，需要使用两次迭代。
下面的例子，第一个for循环遍历shapefile文件中的行对象。在每一行中，需要输出OID(对象标识符)字段，用来指明每一个点对象数组的起始和结束位置。访问每一行对象都会获得一个几何对象，该对象是一个点对象数组。getPart方法可以用来获取几何图形第一部分的点对象数组。第二个循环遍历点对象，输出其x、y坐标值。

import arcpy
from arcpy import env
env.workspace = "C:/Data"
fc = "roads.shp"
cursor = arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"])
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))

运行这个脚本，将以列表的形式输出三个要素中每个折点的x、y坐标。
上述脚本中getPart方法中使用的索引值是0，这时getPart方法只返回几何对象的第一部分。如果没有指定索引值，将会返回一个包含所有点对象数组的数组。

3.处理多部分要素

要素类中的要素可以分为多个部分，我们将其称为多部分要素。
可以使用Describe函数中的shapeType属性来确定一个要素类是否为多部分要素。shapeType的有效返回值时Point、Polygon、Polyline、Multipoint、Multipatch，其中Multipatch表示三位数据。一个要素类时多部分的，并不代表要素类中的每个要素都是多部分的。几何对象的isMultipart属性可以用来确定某个要素是否为多部分的。partCount属性将返回组成某个要素所有部分的数目。
下面的代码说明了如何访问多部分的线要素和多边形要素：

import arcpy
from arcpy import env
fc = "roads.shp"
cursor = arcpy.SearchCursor(fc, ["OID@", "SHAPE@"])
for row in cursor:
	print （”Feature {0}：“.format(row[0])）
	partnum = 0
	for part in row[1]:
		print ("{Part {0}:".format(partnum))
		for point in part:
			print ("{0}, {1}".format(point.X, point.Y))
		partnum += 1

这段脚本对但部分要素和多部分要素都有效。对于单部分要素，所有部分的数目是1，并且for循环中的语法快只运行一次。

4.处理有空洞的多边形

如果一个多边形内含孔洞，那说明这个多边形包含很多环。读取含有孔洞的多边形和之前读取多部分要素的脚本很相似。只需要将读取多部分要素脚本中的第三个for循环替代即可：

import arcpy
from arcpy import env
fc = "roads.shp"
cursor = arcpy.SearchCursor(fc, ["OID@", "SHAPE@"])
for row in cursor:
	print （”Feature {0}：“.format(row[0])）
	partnum = 0
	for part in row[1]:
		print ("{Part {0}:".format(partnum))
		for point in part:
			if point:
				print ("{0}, {1}".format(point.X, point.Y))
			else:
				print "Interior Ring"
		partnum += 1

5.写入几何

通过使用插入和更新游标，可以在要素类中创建新要素或更新现有要素。
可以使用CreateFeature函数来创建一个空的要素类，并用来存储新的点对象，创建空要素类函数的语法如下：

CreateFeatureclass_management(out_path, out_name, {geometry_type}, {template}, {has_m}, {has_z}, {spatial_reference}, {config_keyword}, {spatial_grid_1}, {spatial_grid_2}, {spatial_grid_3})

其中，两个必选参数分别是新要素类的路径和名称。空间参考没有默认值，因此如果没有设置，坐标系统会输出”unknow“。

import arcpy, fileinput, string
from arcpy import env
env.overwriteOutput = True
infile = "C:/Data/newpoly.txt"
fc = "c:/Data/newpoly.shp"
arcpy.CreateFeatureclass_management("C:/Data", fc, "Polygon")

运行这段脚本，就可以创建一个名为newpoly.shp的空要素类。
用于表示多边形折点的点对象可以通过ArcPy的Point类创建。这些点对象需要存储在一个数组中。对象数组可以通过ArcPy的Array类创建。
下面的例子，数组将存储点对象，此外，需要使用插入游标来创建一个新的行对象，即一个新的要素。

cursor = arcpy.da.InsertCursor(fc, ["SHAPE@"])
array = arcpy.Array()
point = arcpy.Point()

接下来使用fileinput模块读取文档

for line in fileinput.input(infile):
	point.ID, point.X, point.Y = line.split()

最后脚本需要遍历文本文档每行的数据，并为每行创建一个点对象。结果会得到一个包含一个21点对象的数组。

import arcpy, fileinput, os
from arcpy import env

env.workspace = "C:/Data"
infile = "C:/Data/points.txt"
fc  = "newpoly.shp"
arcpy.CreateFeatureclass_management("C:/Data", fc, "Polygon")
cursor = arcpy.da.InsertCursor(fc, ["SHAPE@"])
array = arcpy.Array()
point = arcpy.Point()
for line in fileinput.input(infile):
	point.Id, point.X, point.Y = line.split()
	line_array.add(point)
polygon = arcpy.Polygon(array)
cursor.insertRow([polygon])
fileinput.close()
del cursor

6.使用游标设置空间参考

要素类的空间参考描述了要素类的坐标系、空间域和精度。如果要素类中没有设置空间参考，可以使用Define Projection工具来定义该要素类的空间参考。

如果搜索游标的空间参考与要素类的空间参考不一致，那么搜索游标会将要素类的空间参考转换为游标的空间参考。
下面的例子以十进制的形式输出了一个采用国家平面坐标系的点要素类的x、y坐标SearchCursor函数是用来在国家坐标系上创建一个只读游标，但是这个只读游标的空间参考要设置为以十进制表示的地理坐标系。代码如下：

import arcpy
fc = "C:/Data/hospitals.shp"
prjfile = "C:/projections/GCS_NAD_1983.prj"
spatialref = arcpy.SpatialReference(prjfile)
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@"], "", spatialref)

然后，使用open函数创建一个输出文件。通过写模式打开文件，然后写入新的代码就行了

output =open("result.txt", "w")

下一步遍历行，为每行创建一个几何对象，使用write方法把x、y坐标写入到输出文件中。
for row in cursor:
point = row[0]
output.write(str(point.X) + " " +str(point.Y) + “\n”)
坐标以十进制字符串的形式写入。完整代码如下：

import arcpy
from arcpy import env
env.workspace = "C:/Data"
fc = "C:/Data/hospitals.shp"
prjfile = "C:/projections/GCS_NAD_1983.prj"
spatialref = arcpy.SpatialReference(prjfile)
cursor = arcpy.da.SearchCursor(fc, ["SHAPE@"], "", spatialref)
output = open("result.txt", "w")
for row in cursor:
	point = row[0]
	output.write(str(point.X) + " " + str(point.Y) + "\n")
output.close()

在这个例子中，空间参考是使用以及存在的投影文件，它也可以从已有的要素类中获取。

7.使用地理处理工具处理几何对象

7.1几何对象作为输入

地理处理工具的参数通常需要输入要素类，然后优势要素类并不存在，需要通过相关坐标数据创建。这种情况下线创建一个空的要素类，然后使用游标设置其属性值，最后再将其作为地理处理工具的输入要素。这种方法很繁琐。另一种方法是用几何对象代替输入或输出要素，这样就可以使地理处理更简便。
例如，下面的代码根据一系列坐标值，创建了一个几何对象，然后使用这些几何对象作为Buffer工具的输入要素。代码如下：

import arcpy
from arcpy import env
env.workspace = "C:/Data"
coordlist = [[17.0, 20.0], [125.0, 32.0], [4.0, 87.0]]
pointlist = []
for x,y in coordlist:
	point = arcpy.Point(x,y)
	pointgeometry = arcpy.PointGeometry(point)
	pointlist.append(pointgeometry)
arcpy.buffer_analysis(pointlist, "buffer.shp", "10 METERS")

在这个例子中，几何对象通过点对象列表被创建。使用PointGeometry类将点对象创建成几何对象，几何对象存放在列表中，这个列表就成了Buffer工具的输入要素。

7.2几何对象作为输出

几何对象也可以作为地理处理工具的输出。例如，下面的代码使用空的几何对象作为Copy Features的工具输出，其输出结果为几何对象列表》

import arcpy
fc = "C:/Data/roads.shp"
geolist = arcpy.CopyFeatures_management(fc, arcpy.Geometry())
length = 0
for geometry in geolist:
	length += geometry.length
print "Total length:" + length

使用几何对象可以提高代码运行效率，因为使用几何对象可以避免创建临时要素类，同时也可以利用游标遍历所有要素。