Arcpy新增随机高程点、空间插值及批量制图

最新推荐文章于 2024-04-21 10:10:48 发布
最新推荐文章于 2024-04-21 10:10:48 发布
阅读量2.2k 收藏 21
点赞数 6
分类专栏： ArcGIS 文章标签： python
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/wlh2067/article/details/127996535
版权

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

自定义工具

1.计算图幅边界：

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

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

image-20221121142122413

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

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

image-20221121142429470

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

image-20221121143350191

2.创建边界

代码如下：

import os
import arcpy

def creat_box():
	try:
		#workspace
		ws = r"D:\OneDrive\UCAS\courses\python2\final2"
		outfile = os.path.join(ws,"output","region.shp")
		 
		#WGS84
		outSR = 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 point
		array = arcpy.Array()
		array.add(LowerLeft)
		array.add(LowerRight)
		array.add(UpperRight)
		array.add(UpperLeft)
		array.add(LowerLeft)
		#creat polygon
		polygon = arcpy.Polygon(array)	 
		#save
		arcpy.CopyFeatures_management(polygon, outfile)
		#define projection
		arcpy.DefineProjection_management(outfile, outSR)
		print('OK')
	except Exception, e:
		print str(e)

if __name__=="__main__":
     creat_box()

image-20221121152121089

3.边界相交

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

用Raster Domain工具比较快：

image-20221121150456756

最后再裁剪一下：

image-20221121150638577

代码如下：

import os
import arcpy

def 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()

结果如图：

image-20221121152027330

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

4.生成随机点

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

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

image-20221122192251906

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

分别对应

调查点编号调查点名称X坐标Y坐标地面高程地理位置地图编号
codename1xyDEMname2map

建议全用英文

最后用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 = 200
	    outName = '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 = 20
	    arcpy.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 rec
	pStart = 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()

结果如图：

image-20221122193041145

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个高程点结果如图：

image-20221122193357969

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 = 30
    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, "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.nrows
    ncols = sheet.ncols
    print "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 rows
    spatial_reference = arcpy.SpatialReference(4326)
    arcpy.DefineProjection_management(shapefilename, spatial_reference)

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

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

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

总共520个点，结果如图：

image-20221122194656715

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") 
        #WGS84
        outSR = 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 point
        array = arcpy.Array()
        array.add(LowerLeft)
        array.add(LowerRight)
        array.add(UpperRight)
        array.add(UpperLeft)
        array.add(LowerLeft)
        #creat polygon
        polygon = arcpy.Polygon(array)   
        #save
        arcpy.CopyFeatures_management(polygon, outfile)
        #define projection
        arcpy.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 = ws
        numPoints = 200
        outName = '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 = 30
        arcpy.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 rec
    pStart = 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 = 30
    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, "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.nrows
    ncols = sheet.ncols
    print "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 rows
    spatial_reference = arcpy.SpatialReference(4326)
    arcpy.DefineProjection_management(shapefilename, spatial_reference)

def merge_table(shapefilename, shapefilename1):
    excel_table = shapefilename
    random_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

image-20221122195202255

再在ToolBox里Add一个Script

image-20221122195228992

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

image-20221122195609297

Source里面设置creat_table完整代码

image-20221122195657248

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

image-20221122195745515

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

等待几秒后就成功了：

image-20221122200001299

spline方法插值

插值方法见官方文档：

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

ell_size(Optional)The cell size of the output raster that will be created.This parameter can be defined by a numeric value or obtained from an existing raster dataset. If the cell size hasn’t been explicitly specified as the parameter value, the environment cell size value will be used if specified; otherwise, additional rules will be used to calculate it from the other inputs. See the usage for more detail.Analysis Cell Size
spline_type(Optional)The type of spline to be used.REGULARIZED —Yields a smooth surface and smooth first derivatives.TENSION —Tunes the stiffness of the interpolant according to the character of the modeled phenomenon.String

这里用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 <https://opensource.org/licenses/MIT>.

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.00083333333
		splineType = "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是插值后的结果：

image-20221122203845278

image-20221122203853191

制图模板

1.创建格网索引

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

代码如下：

import arcpy
from arcpy import env
import os

ws = 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修改为图幅号：

image-20221123091944184

2.制作制图模板

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

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

image-20221123092113062

将Enable Data Driven Pages勾选

image-20221123092144247

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

image-20221123092454463

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

image-20221123092354311

代码如下:

# 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 <https://opensource.org/licenses/MIT>.
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 = 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 = "Beijing DEM(m)"+"("+graphiccode+")"
    name="D:/OneDrive/UCAS/courses/python2/final2/map/" + graphiccode + ".jpg"
    print name
    arcpy.mapping.ExportToJPEG(mxd,  name, resolution = 300, jpeg_quality =80)
del mxd

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

image-20221123091350698

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

image-20221123091412677

image-20221123091622017

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

地学万事屋
关注
专栏目录
Arcpy实现批量制图
qq_34315872的博客
01-07 3215
最近有学习到Arcgis python编程，对其中的mapping模块比较感兴趣，由于每次制作遥感影像专题图都差不多的做法，大都通过arcgis添加标题、比例尺、指北针和图例等，基本上都是对这几个元素的调整。而借助arcpy.mapping就可以基于一个规整的专题图模板自动递归实现大批量数据的地图制图，对于项目中数据量比较大的情况非常适用。 为了实现整个流程，...
python使用arcpy.mapping模块批量出图
09-21
出图是项目里常见的任务，有的项目甚至会要上百张图片，所以批量出土工具很有必要，这篇文章主要介绍了python使用arcpy.mapping模块批量出图，有兴趣的可以了解一下。
Arcpy自动化批量制图
TJLCY的博客
10-25 1460
import arcpy mxd = arcpy.mapping.MapDocument('./demo_bak.mxd') lyrs = arcpy.mapping.ListLayers(mxd) time = arcpy.mapping.ListLayoutElements(mxd, 'TEXT_ELEMENT', 'time')[0] for i in range(121): shp = './shps/{}.shp'.format(i) cursor = arcpy.SearchC
（四）Arcpy 批量出图
Lucky_云佳的博客
04-21 1427
本文主要介绍了如何利用Arcpy进行批量出图。从已有的江苏省2021年人口分布图，绘制江苏各地级市的人口分布图，一个一个出图重复性操作且耗费时间，本文主要介绍数据驱动等两种方法进行批量出图。
Arcpy新增随机高程点、Spline插值及地图分幅输出
qq_42335815的博客
12-17 2066
实验介绍 基于Python空间数据处理课程的Arcpy大作业，分享一些经历和心得。 题目要求： （1）在“地质调查点基础数据表.xls”中图幅范围内增加100个随机位置的高程点。构建一个shape文件，采用自定义工具的模式，参数有两个：一个是让用户选择excel文件，一个让用户指定新生成的文件名。 （2）采用spline方法插值shape点文件中的高程信息。 （3）将高程图和点图层叠加显示。制作简单的制图模板，将图件按照图幅编号分幅输出为jpg图。 图幅范围计算 要生成随机点，首先应确定随机点的空间范围。
interpolation_Arcpy_批量插值_
10-03
`interpolation_Arcpy_批量插值_`这个标题表明我们将探讨如何使用Python中的`arcpy`模块进行批量的空间插值操作。`arcpy`是Esri公司开发的Geoprocessing工具包，它提供了与ArcGIS Desktop和ArcGIS Pro交互的接口，...
ArcPy结合数据驱动模块的批量制图（GIS思维）
06-19
本文将深入探讨如何利用ArcPy，Esri的Python库，结合数据驱动页面功能来实现高效的批量制图过程。ArcGIS是业界领先的GIS软件，其强大的功能和灵活性使得它成为大数据和云计算环境中的首选工具。 首先，我们要理解...
Arcpy 提取质心、添加属性值、插值、批量出图脚本文件汇总
最新发布
04-21
【内容概要】：文件包含所写Arcpy文章内的python代码脚本以及自定义工具箱Geo_Tools.tbx文件，在WorkSpace文件夹下也包含原始数据JS矢量数据和人口Excel数据，以及批量出图的绘图模板.mxd文件和图层.lyr文件。...
IDW插值批处理arcpy
08-10
本程序基于arcpy编写的Idw批处理程序，程序中具体的参数设置需根据需要自行设置
ArcGIS（ArcPy）脚本excel批量添加字段
07-27
本文将深入探讨如何使用ArcPy脚本来批量地在Excel表格中添加字段，这对于处理大量地理属性数据非常有用。 首先，了解ArcPy的基本结构是必要的。ArcPy包含了多个模块，如arcpy.env用于设置环境变量，arcpy....
高程点内插
04-18
解压密码：nlsh，根据等高线内插高程点，有图文解说。
内插高程点
09-29
利用autolisp语言在两高程点内插入高程点
ArcGIS批量投影栅格工具和批量克里金插值工具
12-04
在ArcGIS10.3版本中，基于模型构建器建立了批量投影栅格工具和批量克里金插值工具，批量投影栅格工具可实现批量投影栅格，批量克里金插值工具可实现对点图层批量插值。ArcGIS10.3及以上版本可以直接打开和使用此工具。
基于自然邻点插值的数据处理方法
09-16
基于自然邻点的插值方法，也请大家尽情下载，希望能帮助大家
arcgispython空间插值_python2 使用arcpy自定义GIS工具箱，实现excel读取、空间插值
weixin_39851977的博客
12-10 703
一、实验目的1.掌握Arcpy基本语法、功能函数2.学会使用Arcpy自定义GIS工具箱3.使用arcpy进行空间插值4.自定义创建的Arcmap制图模板，使用arcpy调节图例、注记、比例尺等位置二、实验内容1.arcpy包自定义工具箱，选择文件，生成新的文件名，按照第4列经、第5列纬度，生成点矢量图形.2.按照第6列地面高程数据进行spline空间插值，生成空间插值后的图层，在第十一章Arcp...
python2 使用arcpy自定义GIS工具箱，实现excel读取、空间插值
Never Settle
01-20 4450
一、实验目的 1.掌握Arcpy基本语法、功能函数 2.学会使用Arcpy自定义GIS工具箱 3.使用arcpy进行空间插值 4.自定义创建的Arcmap制图模板，使用arcpy调节图例、注记、比例尺等位置 二、实验内容 1.arcpy包自定义工具箱，选择文件，生成新的文件名，按照第4列经、第5列纬度，生成点矢量图形. 2.按照第6列地面高程数据进行spline空间插值，生成空间插值后的图层，在第...
基于ArcGIS的批量克里金插值
qq_50291870的博客
07-06 2960
获取到地下水位点位数据，需要进行空间插值以获得空间连续的栅格数据。记录一下操作流程。
利用Arcpy批量图斑生成图片
漓涂
05-10 1416
批量生成图片任务需求环境与数据准备代码实现运行结果 任务需求 实际作业中，往往需要导出一些经典的案例图斑用作参考。而能称之为经典的图斑大部分都会来源于多个不同的数据源中。接下来就实现跨数据源导出指定图斑图片： 环境与数据准备 1、arcgis10.2.2 2、准备一份数据源，此处使用GDB作为数据源。 3、准备提供地图要素参考信息的地图文档.mxd文件,可以在里面设置图斑信息。 代码实现 # encoding:utf-8 import os.path import arcpy import sys imp
arcpy在图幅范围内增加200个随机位置的高程点。
06-08
以下是使用ArcPy在图幅范围内增加200个随机位置的高程点的示例代码： ``` import arcpy import random # 设置工作空间和输出路径 arcpy.env.workspace = r"C:\data" out_path = r"C:\output" # 定义输入DEM和输出...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

地学万事屋

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值