实验介绍
基于Python空间数据处理课程的Arcpy大作业，分享一些经历和心得。
题目要求：
（1）在“地质调查点基础数据表.xls”中图幅范围内增加100个随机位置的高程点。构建一个shape文件，采用自定义工具的模式，参数有两个：一个是让用户选择excel文件，一个让用户指定新生成的文件名。
（2）采用spline方法插值shape点文件中的高程信息。
（3）将高程图和点图层叠加显示。制作简单的制图模板，将图件按照图幅编号分幅输出为jpg图。
图幅范围计算
- 要生成随机点，首先应确定随机点的空间范围。从“地质调查点基础数据表.xls”文件中的高程点均在4个图幅（J50E002007、J50E002008、J50E003007、J50E003008）的范围内，依据我国基本地形图的分幅与编号规则计算各图幅的起始经纬度：
J50表示1:100万地形图的第10行第50列，每行间隔4°，每列间隔6°，其左上角坐标计算：
Ymax=10×4°=40°
E002007表示 1:5万地形图的第2行第7列，每行间隔10’，每列间隔15’， 左上角坐标计算：
Ymax=40°-1×10’=39°50’
图幅经纬度：
在excel中添加随机点，并依据“北京.img”提取对应的随机点高程，保存excel文件
- 在四个图幅范围内分别插入若干个随机高程点，高程值来源于“北京.img”，且保证随机点全部落在有高程值的区域，保存随机点的计数编号、位置信息和图幅编号。
Tip1：原题要求在图幅范围内插入100个随机点，为提高最终插值栅格的准确性，可以增加生成的高程点的数量。这里依据各图幅与“北京.img”的重合程度（下图右），大致确定各图幅范围内随机点的数量：J50E002007，J50E002008插入200个点、J50E003007插入130个点、J50E003008插入170个点。
Tip2：考虑到“北京.img”在部分图幅范围内无值（下图右），因此在提取随机点对应高程值时，删去无值的随机高程点并重新取点，保证有高程值的随机点数量，即随机点全部取在有高程值的区域。
流程图、图幅范围与北京高程图的叠加
随机高程点生成与高值提取代码：
import xlrd
import xlwt
from xlutils.copy import copy
import numpy as np
import arcpy
def open_excel(file="D:\\Desktop\\chinamap\\dizhi.xls", colindex=0):
try:
data = xlrd.open_workbook(file)
datac=copy(data)
table = data.sheet_by_index(colindex)
mysheet=datac.get_sheet(0)
#表格的行数
nrows = table.nrows
j=1
m=0
#分图幅生成不同数量的随机点
graphicode=["J50E002007","J50E002008","J50E003007","J50E003008"]
for gc in graphicode:
minx = 115.5
miny = 39.67 #J50E002007的左下角点
if gc=="J50E002007":
pass
randnum=200 #随机点的数量
elif gc=="J50E002008":
minx = minx + 0.25
randnum=200
elif gc=="J50E003007":
miny = miny -0.17
randnum=130
elif gc=="J50E003008":
minx=minx + 0.25
miny = miny - 0.17
randnum=170
while m < randnum:
#生成随机点的经纬度坐标
randx=np.random.uniform(minx, minx+0.25)
randy=np.random.uniform(miny, miny+0.17)
location = str(randx) + " " + str(randy) #建立位置参数，从北京.img中提取这个位置上的高程值
elevationvalue = arcpy.GetCellValue_management("D:\\Desktop\\chinamap\\北京.img", location, "")
valueOK = str(elevationvalue.getOutput(0)) #判断该随机点上是否有高程值
if valueOK == "NoData":
continue
#如果存在值，则写入excel中
else:
row=nrows+m #行号
code="Rand"+str(j) #随机点编号
mysheet.write(row, 3, u"随机点")
mysheet.write(row,6,randx)
mysheet.write(row,7,randy)
mysheet.write(row,8,int(elevationvalue.getOutput(0)))
mysheet.write(row,0,code)
mysheet.write(row,10,gc)
j=j+1
m=m+1
nrows=nrows+m
m=0
datac.save("D:/Desktop/rand3.xls")
print "finish rand"
except Exception, e:
print str(e)
if __name__=="__main__":
open_excel()
自定义工具，依据Excel文件构建点要素类，添加脚本文件
1. 创建如下工具箱并添加脚本，在脚本中添加脚本文件；设置脚本参数，前两个分别为输入的Excel和输出的shp文件，为必选参数；第三个是投影坐标系的shp文件，为可选参数。
Tip1：PointShpFile参数的方向需设置为Output；设置参数时可在过滤器中设置文件类型，如下图示：（如下图）：
2. 脚本文件编写
创建空的点要素类，添加点编号、名称、高程值、图幅编号等字段：
import arcpy
import time
import xlrd
from arcpy import env
import os
def creatshp(shapefilename):
env.workspace = "D:/Desktop/course/python/arcpy/chinamap"
env.extent = "115.5 39.5 116 39.84"
#分离路径与文件名，以便于输入CreateFeatureclass_management
out_path,out_name=os.path.split(shapefilename)
geometry_type = "POINT"
has_m = "DISABLED"
has_z = "DISABLED"
spatial_reference = arcpy.SpatialReference(4326)
if arcpy.Exists(shapefilename):
arcpy.Delete_management(shapefilename)
print shapefilename + " is deleted!"
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, "", has_m, has_z, spatial_reference)
print "A Shape file is created!"
fieldName1 = "pointCode"
fieldName2 = "name"
fieldName3 = "elevation"
fieldName4 = "graphcode"
fieldLength = 15
arcpy.AddField_management(shapefilename, fieldName1, "TEXT", "", "", fieldLength)
arcpy.AddField_management(shapefilename, fieldName2, "TEXT", "", "", fieldLength)
arcpy.AddField_management(shapefilename, fieldName3, "FLOAT", "", "", fieldLength)
arcpy.AddField_management(shapefilename, fieldName4, "TEXT", "", "", fieldLength)
arcpy.AddMessage("Point shp has created")
从Excel中读取高程点插入到新建的要素类中，并给字段赋值，从自定义工具中获取文件路径参数：
def xlstoarcgis(excelfilename,shapefilename):
print time.strftime('begin time: %I:%M-%S', time.localtime(time.time()))
cursor = arcpy.InsertCursor(shapefilename)
table = open_excel(excelfilename, 0)#用中文路径可能会有问题
nrows = table.nrows
ncols = table.ncols
print "num of rows:{0}, num of cols:{1}".format(nrows, ncols)
for rownum in range(1, nrows):
row = table.row_values(rownum)
if row:
feat = cursor.newRow()
pnt = arcpy.Point()
pnt.X = table.cell(rownum, 6).value
pnt.Y = table.cell(rownum, 7).value
feat.Shape = arcpy.PointGeometry(pnt)
feat.pointCode = table.cell(rownum, 1).value
feat.name = table.cell(rownum, 3).value
feat.elevation = table.cell(rownum, 8).value
feat.graphcode=table.cell(rownum,10).value
try:
print "x:{0} y:{1}".format(pnt.X, pnt.Y)
cursor.insertRow(feat)#插入点要素
except:
print arcpy.GetMessages()
print time.strftime('end time: %I:%M-%S', time.localtime(time.time()))
del cursor
def main():
excelfilename=arcpy.GetParameterAsText(0)
shapefilename=arcpy.GetParameterAsText(1)
shapefilename2 = arcpy.GetParameterAsText(2)
creatshp(shapefilename)
arcpy.env.extent = "115.5 39.5 116 39.84"
xlstoarcgis(excelfilename,shapefilename)
arcpy.Project_management(shapefilename,shapefilename2, 32650)#创建投影坐标系的点要素类
print "elevation shpfile was created successfully"
if __name__=="__main__":
main()
Tip2：此时创建的点要素类为地理坐标系（4326），若有需要，可将其转换为投影坐标系WGS_1984_UTM_Zone_50N（32650）（后续仍使用地理坐标系数据）：
arcpy.Project_management(shapefilename,shapefilename2, 32650)
Tip3：脚本文件中不能出现注释！！！，否则在Arcgis自定义工具中运行时会报错。
Tip4：由于CreateFeatureclass_management方法同时需要路径文件夹和文件名两个参数，而通过自定义工具的一个参数只能获取完整文件路径，因此需要在代码中对该完整文件路径进行分离：
out_path,out_name=os.path.split(shapefilename)
自定义工具运行结果
采用spline插值shape点文件中的高程信息
1. Spline样条函数插值
创建了一个表示插值范围的面要素类qjkExtent，代表了四个图幅的全部范围，同时也是为了方便后续SplineWithBarriers插值方法的使用：
arcpy.CreateFeatureclass_management("D:\\Desktop\\chinamap\\out", "qjkExtent.shp", "POLYGON", "", "DISABLED", "DISABLED", 4326)
coordinates = [115.5,39.5,116,39.84]
LowerLeft = arcpy.Point(coordinates[0], coordinates[1])
LowerRight = arcpy.Point(coordinates[2], coordinates[1])
UpperLeft = arcpy.Point(coordinates[0], coordinates[3])
UpperRight = arcpy.Point(coordinates[2], coordinates[3])
array = arcpy.Array()
array.add(LowerLeft)
array.add(LowerRight)
array.add(UpperRight)
array.add(UpperLeft)
array.add(LowerLeft)
polygon = arcpy.Polygon(array)
arcpy.CopyFeatures_management(polygon, "D:\\Desktop\\chinamap\\out\\qjkExtent.shp")
outSR = arcpy.SpatialReference(4326)
arcpy.DefineProjection_management("D:\\Desktop\\chinamap\\out\\qjkExtent.shp", outSR)
env.extent="D:\\Desktop\\chinamap\\out\\qjkExtent.shp"
Spline插值与SplineWithBarriers插值的对比：由于Spline插值的效果不是很好，因此尝试使用刚才创建的面要素类作为障碍进行SplineWithBarriers插值，并对比二者结果：
env.extent="D:\\Desktop\\chinamap\\out\\qjkExtent.shp"
#arcpy.env.mask = "e:/chinamap/北京.img"
# 之前创建的高程点数据集
inPointFeatures = "D:\\Desktop\\chinamap\\out\\ele1.shp"
zField = "elevation"
splineType = "REGULARIZED"
weight = 0.1
arcpy.CheckOutExtension("Spatial")
# Execute Spline
#使用Spline方法进行插值
#outSpline = Spline(inPointFeatures, zField,0.00083333333)
#采用SplineWithBarriers插值效果更好
outSpline = SplineWithBarriers(inPointFeatures, zField,"D:\\Desktop\\chinamap\\out\\qjkExtent.shp",0.00083333333)
fileName="D:\\Desktop\\chinamap\\out\\sp2"
if arcpy.Exists(fileName):
arcpy.Delete_management(fileName)
outSpline.save("D:\\Desktop\\chinamap\\out\\sp2")
#outSpline.save("D:\\Desktop\\chinamap\\out\\sp1")
#可进行投影转换
#arcpy.DefineProjection_management("D:\\Desktop\\chinamap\\out\\sp4", 32650)
print "processing is over"
结果对比：如下图所示，下方为“北京.img”高程图，左上和右上分别为Spline插值与SplineWithBarriers插值的结果（像元大小一致），可以明显看出SplineWithBarriers插值的结果与原始的高程图更为接近，因此使用SplineWithBarriers插值的结果。
2.裁剪插值栅格
使用ExtractByMask方法提取“北京.img”范围内的插值栅格：
env.workspace = "D:\\Desktop\\chinamap"
arcpy.CheckOutExtension("Spatial")
inRaster = "D:\\Desktop\\chinamap\\out\\sp2"
inMaskData ="D:\\Desktop\\chinamap\\北京.img"
arcpy.CheckOutExtension("Spatial")
outExtractByMask = ExtractByMask(inRaster, inMaskData)
outExtractByMask.save("D:\\Desktop\\chinamap\\out\\elemask")
裁剪结果：
将高程图和点图层叠加显示。制作简单的制图模板，将图件按照图幅编号分幅输出为jpg图。
1. 创建索引格网
因插值时已指定图幅范围，因此在建立格网索引时可直接设置行数和列数，如下图所示，面宽度和面高度自动变为1:6万地形图格网间隔，即15’和10’。
编辑索引格网的属性，将pageName字段修改为对应的图幅编号，并将索引格网添加到数据驱动页面中，启动数据驱动页面:
2. 创建制图模板
将索引格网设置为不可见，在布局视图中添加标题、指北针、比例尺、图例等元素，并调整各元素的位置和样式，使制图模板更加美观。
Tip1：为保证在各驱动页面范围内的图像能够单独完整显示，需要设置数据框的裁剪范围，将其设置为“裁剪至当前数据驱动页面范围”：
3.按图幅编号输出
按创建的制图模板，依据数据驱动页面将高值插值图按图幅编号输出，同时设置标题变化，将标题设置为"北京市分幅高程插值图（分幅编号）"，图片名称也设置为分幅编号：
mxd = arcpy.mapping.MapDocument("D:\\Desktop\\drivepage.mxd")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
mxd.dataDrivenPages.currentPageID = pageNum
graphiccode=str(mxd.dataDrivenPages.pageRow.getValue("pageName"))
print "printing {}".format(mxd.dataDrivenPages.pageRow.getValue("pageName"))
#设置图像标题为图幅编号
for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
if elm.name == "title":
elm.text = "北京市分幅高程插值图"+"("+graphiccode+")"
name="D:/Desktop/chinamap/out/" + graphiccode + ".jpg"
print name
arcpy.mapping.ExportToJPEG(mxd, name, resolution = 300, jpeg_quality =80)
del mxd
结果如下图所示：
代码还是有很多可以改进的地方，望批评指正。
3869
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
