arcpyIDW插值生成等值线等值面

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.SelectLayerByAttribute_management("smoothL", "ADD_TO_SELECTION", "Shape_Length<5")
        # 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]]]"