多波段影像数据区域统计（Python代码）

Tips ：本代码实现ArcGIS中的区域统计功能（Zonal），有两点特色：
1）能够将统计得到的属性值直接赋予矢量文件中；
2）可实现多波段影像区域统计（方法比较笨，先拆分波段，然后进行统计）。

区域统计的类型参考ArcGIS：
MEAN— Calculates the average of all cells in the value raster that belong to the same zone as the output cell.
MAJORITY — Determines the value that occurs most often of all cells in the value raster that belong to the same zone as the output cell.
MAX — Determines the largest value of all cells in the value raster that belong to the same zone as the output cell.
MEDIAN — Determines the median value of all cells in the value raster that belong to the same zone as the output cell.
MIN — Determines the smallest value of all cells in the value raster that belong to the same zone as the output cell.
MINORITY — Determines the value that occurs least often of all cells in the value raster that belong to the same zone as the output cell.
RANGE — Calculates the difference between the largest and smallest value of all cells in the value raster that belong to the same zone as the output cell.
STD — Calculates the standard deviation of all cells in the value raster that belong to the same zone as the output cell.
SUM — Calculates the total value of all cells in the value raster that belong to the same zone as the output cell.
VARIETY — Calculates the number of unique values for all cells in the value raster that belong to the same zone as the output cell. （unique）

# -*- coding: utf-8 -*-
"""
Created on Thu Dec 17 11:39:53 2020@author: xiao_gf"""
import ogr
from osgeo import gdal
import numpy as np
from rasterstats import zonal_statsdef readTif(fileName):dataset = gdal.Open(fileName)if dataset == None:print(fileName+"文件无法打开")return dataset# 保存tif文件
def writeTiff(im_data, im_geotrans, im_proj, path):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif len(im_data.shape) == 2:im_data = np.array([im_data])im_bands, im_height, im_width = im_data.shape#创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype)if(dataset!= None):dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数dataset.SetProjection(im_proj) #写入投影for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del datasetdef zonal(input_vector,in_raster,Stats_type,path):    in_ds =readTif(in_raster)im_width = in_ds.RasterXSize  # 栅格矩阵的列数im_height = in_ds.RasterYSize  # 栅格矩阵的行数im_bands = in_ds.RasterCountim_data = in_ds.ReadAsArray(0, 0, im_width, im_height)im_geotrans = in_ds.GetGeoTransform()im_proj = in_ds.GetProjection()for band in range(0,im_bands):data = im_data[band]rasters = path+'\B{}'.format(band+1)+'.tif'  #输出单波段名称为B1、B2、...,可根据自己需求修改设定writeTiff(data, im_geotrans, im_proj, rasters)shp = ogr.Open(input_vector,1)lyr = shp.GetLayer()   # 添加字段zonal_Field = ogr.FieldDefn("B{}".format(band+1), ogr.OFTReal) #添加的字段名称和单波段文件名称一致，为B1、B2、...,可根据自己需求修改设定zonal_Field.SetPrecision(9)lyr.CreateField(zonal_Field)zonal_method = zonal_stats(input_vector, rasters, stats=[Stats_type])FID = 0for feat in lyr:Index = zonal_method[FID][Stats_type]# print(Index)feat.SetField("B{}".format(band+1), Index)feat.SetFID(FID)lyr.SetFeature(feat)FID +=1if __name__ =="__main__":input_shp = r"D:\DATA\test.shp"  #输入统计的矢量文件input_raster = r"D:\DATA\test_r.tif"  #需要统计的栅格数据（遥感影像）path = r"D:\DATA\band" #获取的单波段数据存储地址Stats_type = 'mean' #区域统计的类型# Stats_type表示统计的类型，包括：'min', 'max', 'mean', 'count', 'sum', 'std', 'median', 'majority', 'minority', 'unique', 'rangezonal(input_shp, input_raster, Stats_type, path)

多波段影像数据区域统计（Python代码）相关推荐

python 代码行数统计工具_使用Python设计一个代码统计工具
问题设计一个程序,用于统计一个项目中的代码行数,包括文件个数,代码行数,注释行数,空行行数.尽量设计灵活一点可以通过输入不同参数来统计不同语言的项目,例如: # type用于指定文件类型 pytho ...
区域统计批处理功能（Python代码）
Tips :本代码实现ArcGIS中的区域统计功能(Zonal),将统计的属性值直接赋予矢量文件属性表,可实现批量统计. # -*- coding: utf-8 -*- ""&qu ...
ArcGIS+Python读取flt文件并进行区域统计
功能描述:通过ARCGIS的Python脚本,读取flt文件格式,然后根据一个shape文件的区域进行统计其平均值,最后读取出来保存到txt文件中. 1.由于统计后生成的为dbf文件,所以读取需要一个 ...
python代码统计字符串中大写字符、小写字符、特殊字符以及数值字符出现的次数
python代码统计字符串中大写字符.小写字符.特殊字符以及数值字符出现的次数 #python代码统计字符串中大写字符.小写字符.特殊字符以及数值字符出现的次数 import restring = & ...
Python代码统计工具
目录 Python代码统计工具声明一. 问题提出二. 代码实现三. 效果验证 Python代码统计工具标签: Python 代码统计声明本文将对<Python实现C代码统计工具(一 ...
使用PYTHON统计项目代码行数
目录一使用PYTHON统计项目代码行数二应用实例注:原创不易,转载请务必注明原作者和出处,感谢支持! 一使用PYTHON统计项目代码行数遇到一个非常小的需求:统计一个项目里头的各类源代码 ...
python乘法表代码注释_Python统计python文件中代码,注释及空白对应的行数示例【测试可用】...
本文实例讲述了Python实现统计python文件中代码,注释及空白对应的行数.分享给大家供大家参考,具体如下: 其实代码和空白行很好统计,难点是注释行 python中的注释分为以#开头的单行注释或 ...
python统计英文单词个数_统计英文单词的个数的python代码及字符串分割
字符串分割 str="a|and|hello|||ab" alist = str.split('|') print alist结果 str="a hello{这里换成5个 ...
c语言字符统计2sdut,山东理工大学SDUT - ACM OJ 题： Python代码及分析
Python基础语法学习完成,先刷基础题100道巩固 ,附题目.代码.知识分析题目:http://acm.sdut.edu.cn/onlinejudge2/index.php/Home/Index ...

多波段影像数据区域统计（Python代码）

多波段影像数据区域统计（Python代码）相关推荐

最新文章

热门文章