标签：raster zone Python value 波段 im output data 影像
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）
import ogr
from osgeo import gdal
import numpy as np
from rasterstats import zonal_stats
def 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_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
elif 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 dataset
def 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.RasterCount
im_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 = 0
for 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 +=1
if __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', 'range
zonal(input_shp, input_raster, Stats_type, path)
标签：raster,zone,Python,value,波段,im,output,data,影像 来源： https://blog.csdn.net/x2434417239/article/details/111353655
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。