# arcpyIDW插值生成等值线等值面

## 功能说明

1. 读取grid浓度、温度、风场等数据
2. IDW插值
3. 插值数据根据插值范围生成等值线
4. 等值线平滑处理
5. 等值线根据插值范围clip生成等值面
6. 等值面与等值线spatial_join赋等值面范围
7. 输出json结果

## 操作环境

1. arcpy环境下(安装了ArcMap，并且已授权所有功能包括3D Analysis等)
2. python环境变量配置为ArcMap附带安装的目录下，我的是C:\Python27\ArcGIS10.4

## 代码实现

``````# encoding: utf-8
import sys
import arcpy
import os
def main():
try:
#python脚本位置
filePath = sys.argv[1]
#等值线范围"0;10;100;1000;10000;100000;1000000;10000000;100000000"
contourValues = sys.argv[2]
#文件名前缀,防止中间生成的文件重复冲突
pre = sys.argv[3]
#数据四至 [[13245.9678, 2879.7319], [13345.9678, 2879.7319], [13345.9678, 2979.7319], [13245.9678, 2979.7319], [13245.9678, 2879.7319]]"
extentCoords = sys.argv[4]
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\Users\Jessie\Documents\ArcGIS\Default.gdb"
# print(sys.argv[1])
pre = "a"+pre
points = pre+"points"
idwraster = pre+"idwraster"
contourL = pre + "contourL"
smoothL = pre + "smoothL"
polygon = pre + "polygon"
polygonPath = "C://Users//Jessie//Documents//ArcGIS//"+polygon+".shp"
test = pre + "test"
extent = pre+"extent"

# create new extent
arcpy.Copy_management("extent", extent)
# extentCursor = arcpy.UpdateCursor("extent", "1=1", "", "", "")
# for _ex in extentCursor:
#     extentCursor.deleteRow(_ex)
# del extentCursor
rows = arcpy.InsertCursor(extent)
geojson_polygon = {
"type": "Polygon",
"coordinates": eval(extentCoords)
}
print('extent created')
print(geojson_polygon)
_extent = arcpy.AsShape(geojson_polygon)
row = rows.newRow()
row.setValue("SHAPE", _extent)
rows.insertRow(row)
del row
del rows

# print(points+idwraster+contourL+smoothL+polygon)

# arcpy.Delete_management("points", "")
# arcpy.Delete_management("idwraster", "")
# arcpy.Delete_management("contourL", "")
# arcpy.Delete_management("smoothL", "")
# arcpy.Delete_management("polygon", "")
# arcpy.Delete_management("test", "")

# print(filePath)
arcpy.RasterToPoint_conversion(filePath, points, "")
# print("toPoint")
arcpy.Idw_3d(points, "grid_code", idwraster, 0.4, 2, 12, "")
# print("IDW")
arcpy.ContourList_3d(idwraster, contourL, contourValues)
# print("ContourL")
arcpy.SmoothLine_cartography(
contourL, smoothL, "PAEK", 10, "NO_FIXED", "NO_CHECK")
# print("Smooth")
# arcpy.DeleteFeatures_management("smoothL")
delCursor = arcpy.UpdateCursor(smoothL, "Shape_Length<5")
for r in delCursor:
delCursor.deleteRow(r)
del delCursor
# print("ContourPolygon")
arcpy.FeatureToPolygon_management(
smoothL+";"+extent, polygon+".shp", "", "ATTRIBUTES", "")

# print("Attribute")
arcpy.SpatialJoin_analysis(polygonPath, smoothL, test,
"JOIN_ONE_TO_MANY", "KEEP_ALL")
# 获取值
cursor = arcpy.SearchCursor(polygonPath, "1=1", "", "", "")
r = []
for _r in cursor:
fid = str(_r.getValue("FID"))
_cursor = arcpy.SearchCursor(test, "TARGET_FID="+fid, "", "", "")
if _cursor != None:
contour = ""
for interFs in _cursor:
contour = contour+str(interFs.getValue("Contour"))+";"
r.append({"contour": contour, "geom": _r.getValue("Shape").JSON})
# r.append({"contour": contour, "geom": _r.getValue("FID")})
del _cursor
print(r)
arcpy.Delete_management(points, "")
arcpy.Delete_management(idwraster, "")
arcpy.Delete_management(contourL, "")
arcpy.Delete_management(smoothL, "")
arcpy.Delete_management(polygon, "")
arcpy.Delete_management(test, "")
arcpy.Delete_management(extent, "")
sys.exit(0)
except Exception as e:
print(e.message)
sys.exit(1)

if __name__ == '__main__':
main()
``````

## 使用方法

``````#在控制台下执行下面代码
Python2.7 route.py "D:\\work\\taiyuan\\data\\000338AIRI.GRD" "0;10;100;1000;10000;100000;1000000;10000000;100000000" "000338AIRI" "[[[13245.9678, 2879.7319], [13345.9678, 2879.7319], [13345.9678, 2979.7319], [13245.9678, 2979.7319], [13245.9678, 2879.7319]]]"``````