>>> import arcpy ... from arcpy import env ... import numpy ... from arcpy.sa import * ... import os ... ... #rasterPath = arcpy.GetParameterAsText(0) #载入栅格数据 ... rasterPath_LST = r'C:\Users\Administrator\Desktop\LstDem_essay\data.gdb\JJJ_LST_1000' #载入温度栅格数据 ... rasterPath_DEM = r'C:\Users\Administrator\Desktop\LstDem_essay\data.gdb\JJJ_DEM_1000' #载入温度栅格数据 ... ... # outraster = arcpy.GetParameterAsText(1) #定义输出地址 ... # Output raster (outraster) #定义输出栅格数据 ... #outraster = os.path.join(os.path.dirname(rasterPath),"LST_Result.tif") #定义输出栅格数据 ... ... #打开栅格数据 ... openFile_LST = Raster(rasterPath_LST) ... openFile_DEM = Raster(rasterPath_DEM) ... ... workpath=r'C:\Users\Administrator\Desktop\LstDem_essay\data.gdb' #指定工作空间路径 ... env.workpath = workpath #定义工作空间 ... env.scratchWorkspace = workpath #定义临时工作空间 ... env.overwriteOutput = True #允许对原始文件覆盖 ... arcpy.env.extent = arcpy.Extent(4728039.31346127, 3896494.23063293, 5242039.31346127, 4661494.23063293) #定义分析研究的范围 ... #mx_min = rasterPath.extent.XMin ... #mx_max = rasterPath.extent.XMax ... #my_min = rasterPath.extent.YMin ... #my_max = rasterPath.extent.YMax ... ... prjfile = r'C:\Users\Administrator\Desktop\LstDem_essay\data\JJJ_Board.prj' #输入栅格图像的空间参考 ... spatialRef = arcpy.SpatialReference(prjfile) ... myref = spatialRef.name ... print myref ... arcpy.env.outputCoordinateSystem = spatialRef #设置环境中的投影坐标系 ... arcpy.env.cellSize = rasterPath_LST #设置环境中栅格的分辨率30米 ... ... #获取栅格中的NoData值 ... noDataValue = openFile_LST.noDataValue ... noDataValue = openFile_DEM.noDataValue ... ... #将栅格数据转换为Numpy ... array_LST = arcpy.RasterToNumPyArray(openFile_LST) ... array_DEM = arcpy.RasterToNumPyArray(openFile_DEM) ... my_array = arcpy.RasterToNumPyArray(openFile_LST) ... ... #波段、行数、列数 ... rowNum_LST,colNum_LST = array_LST.shape ... rowNum_DEM,colNum_DEM = array_DEM.shape ... ... #设置栅格分辨率 ... CELL_size = 1000.0 ... ... #计算地表温度平均值以及S值 ... countN_SUM = 0.0 ... LST_sum1 = 0.0 #命名地表温度累加和 ... LST_sum2 = 0.0 #命名地表温度的平方后累加和 ... LST_average = 31.578 #命名区域内的平均地表温度 ... S = 0.0 #命名S ... for i in range(1,rowNum_LST-1): ... for j in range(1,colNum_LST-1): ... if array_LST[i][j] != 0 : ... countN_SUM += 1 ... LST_sum1 = LST_sum1 + array_LST[i][j] ... LST_sum2 = LST_sum2 + array_LST[i][j] * array_LST[i][j] ... LST_average = LST_sum1 / countN_SUM #区域内的平均地表温度 ... print countN_SUM ... print LST_average ... print LST_sum1 ... print LST_sum2 ... S = numpy.sqrt((LST_sum2/countN_SUM) - (LST_average * LST_average)) #区域内S值 ... ... #计算G值 ... for i in range(1,rowNum_LST-1): ... for j in range(1,colNum_LST-1): ... if array_LST[i][j] != 0: ... ... Distance_1 = 0.0 #三维距离定义 ... Distance_2 = 0.0 #三维距离定义 ... Distance_3 = 0.0 #三维距离定义 ... Distance_4 = 0.0 #三维距离定义 ... Distance_5 = 0.0 #三维距离定义 ... Distance_6 = 0.0 #三维距离定义 ... Distance_7 = 0.0 #三维距离定义 ... Distance_8 = 0.0 #三维距离定义 ... Distance_sum = 0.0 #单位面积内的三维距离和定义 ... ... WX_sum = 0.0 #权重与X的乘积和 ... WS_sum = 0.0 #权重平方的乘积和 ... K = 0.0 #指代S旁边根号下的内容 ... K1 = 0.0 #指代K的第一子部分内容 ... G = 0.0 #指代K的第一子部分内容 ... ... count = 0 ... #计算单位区域内的反距离 ... for m in range(-1,2): ... for n in range(-1,2): ... ... if m==-1 and n==-1: ... H_1 = 0.0 # ... d_1 = 0.0 # ... D_1 = 0.0 # ... ... H_1 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_1 = numpy.square(H_1)+numpy.square(1.0)+numpy.square(1.0) ... D_1 = numpy.sqrt(d_1) ... Distance_1 = 1.0/(D_1*D_1) ... Distance_sum = Distance_sum + Distance_1 ... count += 1 ... elif m==-1 and n==0: ... H_2 = 0.0 # ... d_2 = 0.0 # ... D_2 = 0.0 # ... ... H_2 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_2 = numpy.square(H_2)+numpy.square(1.0) ... D_2 = numpy.sqrt(d_2) ... Distance_2 = 1.0/(D_2*D_2) ... Distance_sum = Distance_sum + Distance_2 ... count += 1 ... elif m==-1 and n==1: ... H_3 = 0.0 # ... d_3 = 0.0 # ... D_3 = 0.0 # ... ... H_3 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_3 = numpy.square(H_3)+numpy.square(1.0)+numpy.square(1.0) ... D_3 = numpy.sqrt(d_3) ... Distance_3 = 1.0/(D_3*D_3) ... Distance_sum = Distance_sum + Distance_3 ... count += 1 ... elif m==0 and n==-1: ... H_4 = 0.0 # ... d_4 = 0.0 # ... D_4 = 0.0 # ... ... H_4 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_4 = numpy.square(H_4)+numpy.square(1.0) ... D_4 = numpy.sqrt(d_4) ... Distance_4 = 1.0/(D_4*D_4) ... Distance_sum = Distance_sum + Distance_4 ... count += 1 ... elif m==0 and n==1: ... H_5 = 0.0 # ... d_5 = 0.0 # ... D_5 = 0.0 # ... ... H_5 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_5 = numpy.square(H_5)+numpy.square(1.0) ... D_5 = numpy.sqrt(d_5) ... Distance_5 = 1.0/(D_5*D_5) ... Distance_sum = Distance_sum + Distance_5 ... count += 1 ... elif m==1 and n==-1: ... H_6 = 0.0 # ... d_6 = 0.0 # ... D_6 = 0.0 # ... ... H_6 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_6 = numpy.square(H_6)+numpy.square(1.0)+numpy.square(1.0) ... D_6 = numpy.sqrt(d_6) ... Distance_6 = 1.0/(D_6*D_6) ... Distance_sum = Distance_sum + Distance_6 ... count += 1 ... elif m==1 and n==0: ... H_7 = 0.0 # ... d_7 = 0.0 # ... D_7 = 0.0 # ... ... H_7 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_7 = numpy.square(H_7)+numpy.square(1.0) ... D_7 = numpy.sqrt(d_7) ... Distance_7 = 1.0/(D_7*D_7) ... Distance_sum = Distance_sum + Distance_7 ... count += 1 ... elif m==1 and n==1: ... H_8 = 0.0 # ... d_8 = 0.0 # ... D_8 = 0.0 # ... ... H_8 = float(array_DEM[i+m][j+n]-array_DEM[i][j])/CELL_size ... d_8 = numpy.square(H_8)+numpy.square(1.0)+numpy.square(1.0) ... D_8 = numpy.sqrt(d_8) ... Distance_8 = 1.0/(D_8*D_8) ... Distance_sum = Distance_sum + Distance_8 ... count += 1 ... elif m==0 and n==0: ... count += 1 ... ... ... #计算反距离权重 ... if count != 0: ... for p in range(-1,2): ... for q in range(-1,2): ... if p==-1 and q==-1: ... Weight_1 = 0.0 ... Weight_1Square = 0.0 ... Weight_1 = Distance_1/Distance_sum ... Weight_1Square = Weight_1*Weight_1 ... WX_sum = WX_sum + Weight_1 * array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_1Square ... elif p==-1 and q==0: ... Weight_2 = 0.0 ... Weight_2Square = 0.0 ... Weight_2 = Distance_2/Distance_sum ... Weight_2Square = Weight_2*Weight_2 ... WX_sum = WX_sum + Weight_2*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_2Square ... elif p==-1 and q==1: ... Weight_3 = 0.0 ... Weight_3Square = 0.0 ... Weight_3 = Distance_3/Distance_sum ... Weight_3Square = Weight_3*Weight_3 ... WX_sum = WX_sum + Weight_3*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_3Square ... elif p==0 and q==-1: ... Weight_4 = 0.0 ... Weight_4Square = 0.0 ... Weight_4 = Distance_4/Distance_sum ... Weight_4Square = Weight_4*Weight_4 ... WX_sum = WX_sum + Weight_4*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_4Square ... elif p==0 and q==1: ... Weight_5 = 0.0 ... Weight_5Square = 0.0 ... Weight_5 = Distance_5/Distance_sum ... Weight_5Square = Weight_1*Weight_5 ... WX_sum = WX_sum + Weight_5*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_5Square ... elif p==1 and q==-1: ... Weight_6 = 0.0 ... Weight_6Square = 0.0 ... Weight_6 = Distance_6/Distance_sum ... Weight_6Square = Weight_6*Weight_6 ... WX_sum = WX_sum + Weight_6*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_6Square ... elif p==1 and q==0: ... Weight_7 = 0.0 ... Weight_7Square = 0.0 ... Weight_7 = Distance_7/Distance_sum ... Weight_7Square = Weight_7*Weight_7 ... WX_sum = WX_sum + Weight_7*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_7Square ... elif p==1 and q==1: ... Weight_8 = 0.0 ... Weight_8Square = 0.0 ... Weight_8 = Distance_8/Distance_sum ... Weight_8Square = Weight_8*Weight_8 ... WX_sum = WX_sum + Weight_8*array_LST[i+p][j+q] ... WS_sum = WS_sum + Weight_8Square ... elif p==0 and q==0: ... Weight_9 = 1.0 ... Weight_8Square = 1.0 ... WX_sum = WX_sum + 1*array_LST[i+p][j+q] ... WS_sum = WS_sum + 1 ... #计算G值 ... if count >=1: ... K1 = (float(countN_SUM) * WS_sum -4.0)/float(countN_SUM-1) ... K = numpy.sqrt(K1) ... G = (WX_sum - LST_average*2)/(S*K) ... my_array[i][j] = float(G) ... outraster = arcpy.NumPyArrayToRaster(my_array,arcpy.Point(4728039.31346127, 3896494.23063293),1000.0,1000.0) ... outraster.save(os.path.dirname(rasterPath_LST) + "_LST_Result.tif") ... openraster = Raster(os.path.dirname(rasterPath_LST) + "_LST_Result.tif") ... env.mask = r'C:\Users\Administrator\Desktop\LstDem_essay\data.gdb\JJJ_Board' #设置环境中掩膜 ... openrasterClip = ApplyEnvironment(openraster) #环境设置中的各种应用 ... openrasterClip.save(os.path.dirname(rasterPath_LST) + "_LST_Result_30.tif")