Python地理数据处理 十八：arcpy批量处理数据之栅格图像的统计分析

1. 初步运行代码

目的：获取文件夹中所有tif格式的图像的统计数据，包括平均值、最大值、最小值、标准差等。
如果想在arcgis中获取，会非常的繁琐且耗时：

废话不多说直接上代码

# -*- coding: cp936 -*-
import arcpy
import os
import glob
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("ImageAnalyst")
arcpy.CheckOutExtension("spatial")

input =  r"D:/arcgisDatasetFrost/"
Output = open('D:/arcgisDatasetFrost/data.csv', 'w')

rasters = glob.glob(os.path.join(input, "*.tif"))
whereClause = "VALUE = -32556"

for ras in rasters:
    outSetNull = SetNull(ras, ras, whereClause)
    meanValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'MEAN')
    stdValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'STD')
    meanValue = meanValueInfo.getOutput(0)
    stdValue = stdValueInfo.getOutput(0)
    print os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n'
    Output.write(os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n')

for ras in rasters:
    outSetNull = SetNull(ras, ras, whereClause)
    stdValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'STD')
    stdValue = stdValueInfo.getOutput(0)
    print os.path.basename(ras).split('_')[0] + ',' + str(stdValue) + '\n'
    Output.write(os.path.basename(ras).split('_')[0] + ',' + str(stdValue) + '\n')

Output.close()
print("All project is OK！")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

批量裁剪

# -*- coding: cp936 -*-
import arcpy
arcpy.CheckOutExtension("ImageAnalyst")
arcpy.CheckOutExtension("spatial")

arcpy.env.workspace = "D:\\arcgisDatasetFrost\\clip"#栅格文件路径
rasters = arcpy.ListRasters("*","tif")#也可以是其他的文件格式，如grd
mask = "D:\\arcgisDatasetFrost\\a.shp"# 填写你的shp文件位置
for raster in rasters:
    print str(raster)
    out = "D:\\arcgisDatasetFrost\\center\\" + raster
    arcpy.gp.ExtractByMask_sa(raster, mask, out)
    print("ma_" + raster + " has done!")
print("All project is OK!!!")
1
2
3
4
5
6
7
8
9
10
11
12
13
14

2. 2023/3/16代码改进

2.1 一次性对多个文件夹中的图层进行统计，并分别生成对应的.csv文件

# -*- coding: cp936 -*-
import arcpy
import os
import glob
import arcpy
import time
from arcpy.sa import *

start = time.clock()
arcpy.CheckOutExtension("ImageAnalyst")
arcpy.CheckOutExtension("spatial")

# 指定文件夹路径
folder_path = r"D:\data" # 填写你的文件路径
# 获取文件夹下所有文件名
file_names = os.listdir(folder_path)

# 将文件名存入向量
for file_name in file_names:
    inws =  r"D:\data\\" + file_name # 输入文件夹路径,与folder_path一致
    OutputFile = open(os.path.join(folder_path, file_name + '.csv'), 'w')#修改文件名后缀为.csv
    #OutputFile = open(r'D:\data\1.csv', 'w')

    rasters = glob.glob(os.path.join(inws, "*.tif"))
    whereClause = "VALUE = -32556"
    for ras in rasters:
        outSetNull = SetNull(ras, ras, whereClause)
        meanValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'MEAN')
        #stdValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'STD')
        meanValue = meanValueInfo.getOutput(0)
        #stdValue = stdValueInfo.getOutput(0)
        print os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n'
        OutputFile.write(os.path.basename(ras).split('_')[0] + ',' + str(meanValue) + '\n')

	for ras in rasters:
    	outSetNull = SetNull(ras, ras, whereClause)
    	stdValueInfo = arcpy.GetRasterProperties_management(outSetNull, 'STD')
    	stdValue = stdValueInfo.getOutput(0)
    	print os.path.basename(ras).split('_')[0] + ',' + str(stdValue) + '\n'
    	OutputFile.write(os.path.basename(ras).split('_')[0] + ',' + str(stdValue) + '\n')
    
    OutputFile.close()
    print("Subprojects is OK! ")
    print("--------------------------------------------------------------------------------")
    print("--------------------------------------------------------------------------------")

end = time.clock()  
print("All project is OK!!!")
print (str(end-start))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49