arcpy练习（污染排放点污染浓度预测）

一、题目：

二、代码（内有注释，易懂）

# coding:utf-8
import arcpy
from arcpy.sa import *
import math

# 烟囱高度（m）
H = 50
# 排放源强（mg/s）
Q = 150
# 风向，风速（m/s），正西风不用进行坐标轴变换
u = 2
# 点源坐标x0，y0
x0, y0 = 0.0, 0.0
# 存放预测点坐标的list
point_geo_list = []

'''
注意：
    地面到任意污染点的高度（m），z = 0即地面点，z = 50即为烟囱高度;
    大气稳定度等级为A级，My，Mz为y和z方向的弥漫系数,  x为预测点在X方向上到污染源的距离（m）;
    c为预测点点的污染物浓度;
    
    我们以烟囱排放源为中心，向四周进行污染浓度计算，当污染浓度 <50mg/s 时停止计算，在x、y和z方向上计算步长为 1m;
    最后我们对生成的点关于x轴对称添加，对生成shp进行克里金插值显示
'''

# 计算弥漫系数
def get_My_Mz(x):
    My = 0.22 * x / math.sqrt(1 + 0.0001 * x)
    Mz = 0.2 * x
    return My, Mz

# 计算预测点c的污染浓度
def get_c(My, Mz, y, z):
    c = (Q/(u*My*Mz*math.pi)) * (math.exp(-y*y/(2*My*My))) * (math.exp(-(z-H)*(z-H)/(2*Mz*Mz)) + math.exp(-(z+H)*(z+H)/(2*Mz*Mz)))
    return round(c, 4)

# 将坐标放入point_geo_list
def create_shp_point(point_x, point_y, point_z):
    point = arcpy.Point()
    point.X = point_x
    point.Y = point_y
    point.Z = point_z
    point_geo = arcpy.PointGeometry(point)
    point_geo_list.append(point_geo)

# 根据坐标点创建shp，添加预测点浓度字段和值
def create_shp(c):
    shppath = r"C:\Users\Administrator\Desktop\point.shp"
    # 使用point_geo_list创建shp
    if arcpy.Exists(shppath):
        arcpy.Delete_management(shppath)
        arcpy.CopyFeatures_management(point_geo_list, shppath)
    else:
        arcpy.CopyFeatures_management(point_geo_list, shppath)

    # 向point.shp中添加浓度值字段
    arcpy.AddField_management(shppath, "c", "FLOAT", 100, 50)
    rows = arcpy.UpdateCursor(shppath)

    # 添加字段值
    cursor = arcpy.UpdateCursor(shppath)
    i = 1
    for my_row in cursor:
        my_row.setValue('c', c[i])
        cursor.updateRow(my_row)
        i += 1

# 对point.shp进行克里金插值输出
def Kriging_shp():
    outKrig = Kriging(r"C:\Users\Administrator\Desktop\point.shp", "c", KrigingModelOrdinary("CIRCULAR"))
    outKrig.save(r"C:\Users\Administrator\Desktop\Kriging_point")

# 主函数
if __name__ == "__main__":
    point_list = []
    c_list = []
    c = Q
    for z in (sorted(range(1, 51, 1), reverse=True)):       # 对1到50list倒序，从高度 z 为50开始计算
        while c > 10:                                       # 浓度<10时 y 就不循环了
            
            while c > 10:                                   # 浓度<10时 x 就不循环了
                point_list.append([x0, y0, z])
                c_list.append(c)
                x0 = round(x0 + 1, 2)
                My, Mz = get_My_Mz(x0)
                c = get_c(My, Mz, y0, z)
            
            x0 = 1                                          # x不能在0出，如果在那么My，Mz为0，c 的计算就有问题
            y0 = round(y0 + 0.2, 2)
            My, Mz = get_My_Mz(x0)
            c = get_c(My, Mz, y0, z)
            
    for point in point_list[1:]:
        create_shp_point(point[0], point[1], point[2])
        
    for point in point_list[1:]:                            # 锦上添花，y在负半轴上的浓度计算，关于x轴对称
        create_shp_point(point[0], -point[1], point[2])
        
    create_shp(c_list + c_list)
    
    Kriging_shp()

三、结果图：