Python实现栅格批量统计

Python实现栅格批量统计

---

利用ZonalStatisticsAsTable工具实现按照矢量数据特定字段进行栅格批量统计。

---

#This code is used for the  raster batch statistics.
#Author: Changqing Guo
# coding: utf-8
import numpy as np
import arcpy
from arcpy import sa
from arcpy.sa import *
from arcpy import env
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
# 工作空间路径
env.workspace = "F:/Data/extend"

# 修改Raster数据的路径
#in_ULC = r"F:/Data/extend/ISA.tif"

# 修改输出结果的路径
out_Path = "F:/Data/Table/"
# 输入所用的行政分区
Region = "F:/Data/Global.shp"

tiffiles = arcpy.ListFiles("*.tif")
for filename in tiffiles:
    ULC = env.workspace + "/" + filename
    field = filename[:-4]
    shp_qu = out_Path + filename[:-4]
    # 关联字段
    joinID = "value"
    # 统计类型很有多种，自行按需要选择
    arcpy.gp.ZonalStatisticsAsTable_sa(Region, joinID, ULC, shp_qu, "DATA", "MEAN")
    # 在shp文件中增加列的名称，根据需要选择字段类型
    arcpy.AddField_management(shp_qu, field, "DOUBLE")
    arcpy.CalculateField_management(shp_qu, field, "[MEAN]", "VB", "")
    arcpy.JoinField_management(Region, joinID, shp_qu, joinID, [field])

---

先裁剪，后实现按照矢量数据特定字段进行栅格批量统计。

---

#This code is used for the  raster batch statistics.
#Author:Shu Zhang,Changqing Guo
# coding: utf-8
import numpy as np
import arcpy
from arcpy import sa
from arcpy.sa import *
from arcpy import env
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

env.workspace = "F:/Data/Table/"

# 修改栅格数据路径
ULC = r"F:/Data/extend/GUGS_2010.tif"
#扩展区边界
in_GUB = "F:/Data/result/GUB_0010.tif"

# 修改输出结果的路径
out_Path = "F:/Data/Table/"
# 输入行政分区
Region = "F:/Global.shp"

# 关联字段
joinID = "value"
# 在shp文件中增加列的名称
field = "GUGS_0010"

#裁剪
ULC_file = arcpy.Raster(in_ULC)
ULC_clip = ExtractByMask(ULC_file,in_GUB)
ULC_clip.save("F:/Data/GUGS_0010.tif")

shp_qu = out_Path + "GUGS_0010"
# 关联字段
joinID = "value"
# 统计类型很有多种，自行按需要选择
arcpy.gp.ZonalStatisticsAsTable_sa(Region, joinID, ULC, shp_qu, "DATA", "MEAN")
# 在shp文件中增加列的名称，根据需要选择字段类型
arcpy.AddField_management(shp_qu, field, "DOUBLE")
arcpy.CalculateField_management(shp_qu, field, "[MEAN]", "VB", "")
arcpy.JoinField_management(Region, joinID, shp_qu, joinID, [field])