（1）在“地质调查点基础数据表.xls”中图幅范围内增加200个随机位置的高程点。构建一个shape文件，采用自定义工具的模式，参数有两个：一个是让用户选择excel文件，一个让用户指定新生成的文件名。
（2）采用spline方法插值shape点文件中的高程信息。
（3）将高程图和点图层叠加显示。制作简单的制图模板，将图件按照1:1万图幅编号分幅输出为jpg图。

自定义工具

1.计算图幅边界：

打开“地质调查点基础数据表.xls”，发现一共有四种图幅号：

一共有四种：J50E002007、J50E002008、J50E003007、J50E003008

接下来就是根据图幅号计算经纬度了，这又让我回忆起大一下学期的地图学考试，抓耳挠腮计算的我。

But，四年后我带着科技回来了，没想到吧：

算完大概是这样，感兴趣的同学可以自己算一下，我觉得没什么技术含量就不算了。或者也可以自己写个python小程序搞一下：

2.创建边界

代码如下：

import os
import arcpydef creat_box():try:#workspacews = r"D:\OneDrive\UCAS\courses\python2\final2"outfile = os.path.join(ws,"output","region.shp")#WGS84outSR = arcpy.SpatialReference(4326)coordinates = [115.5, 39.5, 116, 39.83]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])#creat pointarray = arcpy.Array()array.add(LowerLeft)array.add(LowerRight)array.add(UpperRight)array.add(UpperLeft)array.add(LowerLeft)#creat polygonpolygon = arcpy.Polygon(array)	 #savearcpy.CopyFeatures_management(polygon, outfile)#define projectionarcpy.DefineProjection_management(outfile, outSR)print('OK')except Exception, e:print str(e)if __name__=="__main__":creat_box()

3.边界相交

为了避免采样点搞到北京外面去，还需要北京与其相交。在这之前还要DEM转为矢量。

用Raster Domain工具比较快：

最后再裁剪一下：

代码如下：

import os
import arcpydef clip_region():try:region = r'D:\OneDrive\UCAS\courses\python2\final2\output\region.shp'inRaster = r'D:\OneDrive\UCAS\courses\python2\final2\input\beijingDEM.tif'outPolygons = r'D:\OneDrive\UCAS\courses\python2\final2\output\beijing.shp'outPolygons2 = r'D:\OneDrive\UCAS\courses\python2\final2\output\region_final.shp'beijings = arcpy.RasterDomain_3d(inRaster, outPolygons, "POLYGON")arcpy.Clip_analysis(beijings, region, outPolygons2)print('OK')except Exception, e:print str(e)if __name__=="__main__":clip_region()

结果如图：

接下来采样的矢量设置这里就行

4.生成随机点

接下来“增加200个随机位置的高程点”

用CreateRandomPoints_management工具即可，主要是创建点后，需要添加一系列字段来匹配已有的表

由于“调查点编号”和“调查点原编号”是一个，就都设置为code

分别对应

建议全用英文

最后用CalculateField_management来为各字段赋值，code是连续的编号，name1和name2都是统一的直接赋值，x，y是通过计算几何获得，map要写个if判断

代码如下：

import os
import arcpy
#arcpy.CreateFileGDB_management('D:/OneDrive/UCAS/courses/python2/final2/output', 'table.gdb')
def creat_shp():try:outGDB = "D:/OneDrive/UCAS/courses/python2/final2/output"numPoints = 200outName = 'random_point'fcExtent1 = r"D:\OneDrive\UCAS\courses\python2\final2\output\region_final.shp"arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4326)random = arcpy.CreateRandomPoints_management(outGDB, outName, fcExtent1, '', numPoints)arcpy.env.outputCoordinateSystem = ""print "A Shape file is created!"fieldName1 = "code"fieldName2 = "name1"fieldName3 = "x"fieldName4 = "y"fieldName5 = "DEM"fieldName6 = "name2"fieldName7 = 'map'fieldLength = 20arcpy.AddField_management(random, fieldName1, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName2, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName3, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName4, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName5, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName6, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName7, "TEXT", "", "", fieldLength)arcpy.AddMessage("Point shp has created")arcpy.CalculateField_management(random,"code",expression,"PYTHON_9.3",codeblock)arcpy.CalculateField_management(random,"name1",expression2,"PYTHON_9.3")arcpy.CalculateField_management(random,"x",'!shape.extent.XMax!',"PYTHON_9.3")arcpy.CalculateField_management(random,"y",'!shape.extent.YMax!',"PYTHON_9.3")arcpy.CalculateField_management(random,"name2",expression3,"PYTHON_9.3")arcpy.CalculateField_management(random,"map",expression4,"PYTHON_9.3",codeblock2)except Exception, e:print str(e)codeblock = '''
rec=0
def autoIncrement():global recpStart = 1 pInterval = 1 if (rec == 0): rec = pStart else: rec = rec + pInterval return rec
'''
codeblock2 = '''
def getmap(x,y):if (x < 115.75) & (y > 39.67):return 'J50E002007'elif (x < 115.75) & (y < 39.67):return 'J50E003007'elif (x > 115.75) & (y > 39.67):return 'J50E002008'else:return 'J50E003008'
'''
expression = "'Drandom' + str(autoIncrement())"
expression2 = "'随机点北'"
expression3 = "'fangshan_north'"
expression4 = 'getmap(!x!, !y!)'if __name__=="__main__":creat_shp()

结果如图：

5.提取高程值

直接ExtractValuesToPoints提取高程值，注意这会生成新的一列。通过字段计算器把这列赋给DEM，再删除这列，代码如下：

import os
import arcpy
from arcpy.sa import *
def extract_DEM():try:arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4326)in_features = ExtractValuesToPoints("D:/OneDrive/UCAS/courses/python2/final2/output/random_point.shp", "D:/OneDrive/UCAS/courses/python2/final2/input/beijingDEM.tif","D:/OneDrive/UCAS/courses/python2/final2/output/DEM_points.shp","INTERPOLATE", "ALL")arcpy.env.outputCoordinateSystem = ""arcpy.CalculateField_management(in_features,"DEM",'!RASTERVALU!',"PYTHON_9.3")
arcpy.DeleteField_management(in_features, 'beijingDEM')print('OK')except Exception, e:print str(e)
if __name__=="__main__":extract_DEM()

200个高程点结果如图：

6.将excel表转shp

这一步开始有了自定义工具：【自定义工具的模式，参数有两个：一个是让用户选择excel文件，一个让用户指定新生成的文件名。】

所以有了额外的输入变量表和输出变量shp，之后配置到工具中：

首先创建一张空表，代码如下，定义好字段，使之与生成的随机点一致：

import arcpy
import xlrd
from arcpy import env
import os
from arcpy.sa import *def creatshp(shapefilename):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 = "code"fieldName2 = "name1"fieldName3 = "x"fieldName4 = "y"fieldName5 = "DEM"fieldName6 = "name2"fieldName7 = 'map'fieldLength = 30arcpy.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, "FLOAT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName5, "FLOAT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName6, "TEXT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName7, "TEXT", "", "", fieldLength)arcpy.AddMessage("Point shp has created")

然后excel转点：

我这里用了几何令牌的方法’SHAPE@XY’，并加入多字段一并写入。

几何令牌可用于访问特定几何信息，通常访问完整几何更加耗时，所以指定返回项目更加关心的内容。下面 是常用的几何令牌。

def xlstoarcgis(excelfilename,shapefilename):pathdir = os.path.split(shapefilename)[0]filename = os.path.split(shapefilename)[1]table = xlrd.open_workbook(excelfilename)sheet = table.sheet_by_name('data')nrows = sheet.nrowsncols = sheet.ncolsprint "num of rows:{0}, num of cols:{1}".format(nrows, ncols)x = sheet.col_values(6, 1, nrows)y = sheet.col_values(7, 1, nrows)code = sheet.col_values(1, 1, nrows)name1 = sheet.col_values(3, 1, nrows)DEM = sheet.col_values(8, 1, nrows)name2 = sheet.col_values(9, 1, nrows)Map = sheet.col_values(10, 1, nrows)rows = arcpy.da.InsertCursor(shapefilename, ['name1', 'x', 'y', 'code','name2', 'DEM', 'Map', 'SHAPE@XY'])for i in range(0,nrows-1):pointrow = (name1[i], x[i], y[i], code[i], name2[i], DEM[i], Map[i], (x[i], y[i]))rows.insertRow(pointrow)del rowsspatial_reference = arcpy.SpatialReference(4326)arcpy.DefineProjection_management(shapefilename, spatial_reference)

最后将随机点（200个）与Excel点（320个）合并为一个

由于这两个shp点的所有字段名称、格式、大小都相同，合并会自动对应的：

def merge_table(shapefilename, shapefilename1):excel_table = shapefilenamerandom_table = os.path.join(ws,"random_point.shp")arcpy.Merge_management([excel_table, random_table], shapefilename1)

总共520个点，结果如图：

7.主函数：

主函数只需要将上面提到的函数依次调用：

if __name__=="__main__":excelfilename=arcpy.GetParameterAsText(0)shapefilename1=arcpy.GetParameterAsText(1)shapefilename=os.path.join(ws,"excel_table.shp")creat_box()clip_region()creat_random()extract_DEM()creatshp(shapefilename)xlstoarcgis(excelfilename,shapefilename)merge_table(shapefilename, shapefilename1)

8.完整代码：

import arcpy
import xlrd
from arcpy import env
import os
from arcpy.sa import *
rs = r"D:\OneDrive\UCAS\courses\python2\final2\input"
ws = r"D:\OneDrive\UCAS\courses\python2\final2\output"
def creat_box():try:#workspace#ws = r"D:\OneDrive\UCAS\courses\python2\final2"#outfile = os.path.join(ws,"output","region.shp")outfile = os.path.join(ws,"region.shp") #WGS84outSR = arcpy.SpatialReference(4326)coordinates = [115.5, 39.5, 116, 39.83]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])#creat pointarray = arcpy.Array()array.add(LowerLeft)array.add(LowerRight)array.add(UpperRight)array.add(UpperLeft)array.add(LowerLeft)#creat polygonpolygon = arcpy.Polygon(array)   #savearcpy.CopyFeatures_management(polygon, outfile)#define projectionarcpy.DefineProjection_management(outfile, outSR)print('OK')except Exception, e:print str(e)def clip_region():try:region = os.path.join(ws,"region.shp")inRaster = os.path.join(rs,'beijingDEM.tif') outPolygons = os.path.join(ws,'beijing.shp')outPolygons2 = os.path.join(ws,'region_final.shp')beijings = arcpy.RasterDomain_3d(inRaster, outPolygons, "POLYGON")arcpy.Clip_analysis(beijings, region, outPolygons2)print('OK')except Exception, e:print str(e)def creat_random():try:outGDB = wsnumPoints = 200outName = 'random_point'fcExtent1 = os.path.join(ws,"region_final.shp")arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4326)random = arcpy.CreateRandomPoints_management(outGDB, outName, fcExtent1, '', numPoints)arcpy.env.outputCoordinateSystem = ""print "A Shape file is created!"fieldName1 = "code"fieldName2 = "name1"fieldName3 = "x"fieldName4 = "y"fieldName5 = "DEM"fieldName6 = "name2"fieldName7 = 'map'fieldLength = 30arcpy.AddField_management(random, fieldName1, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName2, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName3, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName4, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName5, "FLOAT", "", "", fieldLength)arcpy.AddField_management(random, fieldName6, "TEXT", "", "", fieldLength)arcpy.AddField_management(random, fieldName7, "TEXT", "", "", fieldLength)arcpy.AddMessage("Point shp has created")arcpy.CalculateField_management(random,"code",expression,"PYTHON_9.3",codeblock)arcpy.CalculateField_management(random,"name1",expression2,"PYTHON_9.3")arcpy.CalculateField_management(random,"x",'!shape.extent.XMax!',"PYTHON_9.3")arcpy.CalculateField_management(random,"y",'!shape.extent.YMax!',"PYTHON_9.3")arcpy.CalculateField_management(random,"name2",expression3,"PYTHON_9.3")arcpy.CalculateField_management(random,"map",expression4,"PYTHON_9.3",codeblock2)except Exception, e:print str(e)codeblock = '''
rec=0
def autoIncrement():global recpStart = 1 pInterval = 1 if (rec == 0): rec = pStart else: rec = rec + pInterval return rec
'''
codeblock2 = '''
def getmap(x,y):if (x < 115.75) & (y > 39.67):return 'J50E002007'elif (x < 115.75) & (y < 39.67):return 'J50E003007'elif (x > 115.75) & (y > 39.67):return 'J50E002008'else:return 'J50E003008'
'''expression = "'Drandom' + str(autoIncrement())"
expression2 = "'random_north'"
expression3 = "'fangshan_north'"
expression4 = 'getmap(!x!, !y!)'
def extract_DEM():try:arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(4326)in_point_features = os.path.join(ws,"random_point.shp")in_rasters = os.path.join(rs,"beijingDEM.tif")in_features = ExtractMultiValuesToPoints(in_point_features, [in_rasters], "BILINEAR")arcpy.env.outputCoordinateSystem = ""arcpy.CalculateField_management(in_features,"DEM",'!beijingDEM!',"PYTHON_9.3")arcpy.DeleteField_management(in_features, 'beijingDEM')print('OK')except Exception, e:print str(e)def creatshp(shapefilename):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 = "code"fieldName2 = "name1"fieldName3 = "x"fieldName4 = "y"fieldName5 = "DEM"fieldName6 = "name2"fieldName7 = 'map'fieldLength = 30arcpy.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, "FLOAT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName5, "FLOAT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName6, "TEXT", "", "", fieldLength)arcpy.AddField_management(shapefilename, fieldName7, "TEXT", "", "", fieldLength)arcpy.AddMessage("Point shp has created")def xlstoarcgis(excelfilename,shapefilename):pathdir = os.path.split(shapefilename)[0]filename = os.path.split(shapefilename)[1]table = xlrd.open_workbook(excelfilename)sheet = table.sheet_by_name('data')nrows = sheet.nrowsncols = sheet.ncolsprint "num of rows:{0}, num of cols:{1}".format(nrows, ncols)x = sheet.col_values(6, 1, nrows)y = sheet.col_values(7, 1, nrows)code = sheet.col_values(1, 1, nrows)name1 = sheet.col_values(3, 1, nrows)DEM = sheet.col_values(8, 1, nrows)name2 = sheet.col_values(9, 1, nrows)Map = sheet.col_values(10, 1, nrows)rows = arcpy.da.InsertCursor(shapefilename, ['name1', 'x', 'y', 'code','name2', 'DEM', 'Map', 'SHAPE@XY'])for i in range(0,nrows-1):pointrow = (name1[i], x[i], y[i], code[i], name2[i], DEM[i], Map[i], (x[i], y[i]))rows.insertRow(pointrow)del rowsspatial_reference = arcpy.SpatialReference(4326)arcpy.DefineProjection_management(shapefilename, spatial_reference)def merge_table(shapefilename, shapefilename1):excel_table = shapefilenamerandom_table = os.path.join(ws,"random_point.shp")arcpy.Merge_management([excel_table, random_table], shapefilename1)if __name__=="__main__":excelfilename=arcpy.GetParameterAsText(0)shapefilename1=arcpy.GetParameterAsText(1)shapefilename=os.path.join(ws,"excel_table.shp")creat_box()clip_region()creat_random()extract_DEM()creatshp(shapefilename)xlstoarcgis(excelfilename,shapefilename)merge_table(shapefilename, shapefilename1)

代码配置：

第六行和第七行改成你的input和output，代码、环境、数据下载见文末

其中input里有北京DEM和我的Toolbox，output是一个空文件夹

9.自定义工具的配置

在catalog里new一个Toolbox

再在ToolBox里Add一个Script

右键刚刚生成的Scripts，Properties——General命名为excel2table

Source里面设置creat_table完整代码

Parameter的名字一定不要输错，excelfilename对应的应该是Table类型（我也没找到更好的类型，文件类型不可以），shapefilename对应的是Shapefile类型

配置好后双击这个工具，输入excel输入和shapefile输出（excel我把名字改成英文了，避免中文报错）

等待几秒后就成功了：

spline方法插值

插值方法见官方文档：

https://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/spline.htm

这里用REGULARIZED方法，可自行更改：

# Copyright (c) 2022 Longhao Wang. hhuwlh@163.com All rights reserved.
# This work is licensed under the terms of the MIT license.  
# or a copy, see .import arcpy
from arcpy import env
import os
from arcpy.sa import *
ws = r"D:\OneDrive\UCAS\courses\python2\final2\output"
def spline():try:inPointFeatures = os.path.join(ws,"data.shp") region = os.path.join(ws,"region_final.shp") zField = "DEM"outRaster = os.path.join(ws,"DEMspline.tif")extent = os.path.join(ws,"region_final.shp")cellSize = 0.00083333333splineType = "REGULARIZED"spline = arcpy.Spline_3d(inPointFeatures, zField, outRaster, cellSize, splineType)outExtractByMask = ExtractByMask(spline, region)extract = os.path.join(ws,"DEMspline_extract.tif")outExtractByMask.save(extract)except Exception, e:print str(e)if __name__=="__main__":spline()print('OK')

下图1是原始DEM，图2是插值后的结果：

制图模板

1.创建格网索引

格网索引工具能生成规则的矩形格网，一般用于制图：

代码如下：

import arcpy
from arcpy import env
import osws = r"D:\OneDrive\UCAS\courses\python2\final2\output"
env.outputCoordinateSystem = arcpy.SpatialReference(4326)# Set local variables
outFeatureClass = os.path.join(ws,"grid.shp") 
polygonWidth = "0.25"
polygonHeight= "0.166666667"
originCoord = "115.5 39.5"
numberRows = "2"
numberColumns = "2"# Execute GridIndexFeatures
arcpy.GridIndexFeatures_cartography(outFeatureClass, "", "", "", "",polygonWidth, polygonHeight, originCoord,numberRows, numberColumns)
expression = 'getmap(!PageNumber!)'
codeblock = '''
def getmap(x):if (x == 1):return 'J50E002007'elif (x == 3):return 'J50E003007'elif (x == 2):return 'J50E002008'else:return 'J50E003008'
'''
arcpy.CalculateField_management(outFeatureClass,"PageName",expression,"PYTHON_9.3",codeblock)

该工具的原理是创建2×2的格网，并且将PageName修改为图幅号：

2.制作制图模板

接下来制作自己的制图模板：

接下来设置驱动页面和显示范围

将Enable Data Driven Pages勾选

制作一张好看的地图模板（自行发挥），将数据裁剪至数据驱动范围：

设置文字标题的Element Name，方便在代码中批量修改

代码如下:

# Copyright (c) 2022 Longhao Wang. hhuwlh@163.com All rights reserved.
# This work is licensed under the terms of the MIT license.  
# or a copy, see .
import os
import arcpy
from arcpy.sa import *
mxd = arcpy.mapping.MapDocument("D:/OneDrive/UCAS/courses/python2/final2/workspace.mxd")
for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):mxd.dataDrivenPages.currentPageID = pageNumgraphiccode=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 = "Beijing DEM(m)"+"("+graphiccode+")"name="D:/OneDrive/UCAS/courses/python2/final2/map/" + graphiccode + ".jpg"print namearcpy.mapping.ExportToJPEG(mxd,  name, resolution = 300, jpeg_quality =80)
del mxd

接下来就可以运行代码了：

Load之后回车，就打印好了，结果如图：

所有数据和代码见下面公众号