Python进行.tif格式数据波段组合;多波段tif数据波段分离
Python进行.tif格式数据波段组合;多波段tif数据波段分离.tif格式数据波段组合多波段.tif数据波段分离阔别已久,我又来了;啊啊啊……读了研究生,忙到天昏地暗(哈哈哈,这都是借口,主要是花了一定时间适应研究生生活)这方面的内容我也主要是应用在遥感领域,我们知道在遥感应用中,大多数的图像数据格式为.tif格式的。进行波段组合使用ArcGIS里面的Composite Bands工具是可以实
·
Python进行.tif格式数据波段组合;多波段tif数据波段分离
阔别已久,我又来了;啊啊啊……
读了研究生,忙到天昏地暗(哈哈哈,这都是借口,主要是花了一定时间适应研究生生活)
这方面的内容我也主要是应用在遥感领域,我们知道在遥感应用中,大多数的图像数据格式为.tif格式的。进行波段组合使用ArcGIS里面的Composite Bands工具是可以实现的,波段分离也可以使用一些工具或巧妙的方法分开。但当你使用Python程序的话,就可以简化这个软件点击和等待的过程,只要修改2个路径就可以了,岂不乐载?
.tif格式数据波段组合
import os
import numpy as np
from osgeo import gdal
from skimage._shared.utils import safe_as_int
class IMAGE:
# 读图像文件
def read_img(self,filename):
dataset = gdal.Open(filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
# im_bands = dataset.RasterCount # 波段数
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵,左上角像素的大地坐标和像素分辨率
im_proj = dataset.GetProjection() # 地图投影信息,字符串表示
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
del dataset #关闭对象dataset,释放内存
# return im_width, im_height, im_proj, im_geotrans, im_data,im_bands
return im_proj, im_geotrans, im_data, im_width,im_height
# 遥感影像的存储
# 写GeoTiff文件
def write_img(self,filename, im_proj, im_geotrans, im_data):
# 判断栅格数据的数据类型
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_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape #没看懂
# 创建文件时 driver = gdal.GetDriverByName("GTiff"),数据类型必须要指定,因为要计算需要多大内存空间。
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(im_proj) # 写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
if __name__ == "__main__":
# os.chdir(r'') #切换路径到待处理图像所在文件夹
run=IMAGE()
# 第一步
proj, geotrans, data1, row1, column1 = run.read_img(r'***.tif') # 读数据
proj, geotrans, data2, row2, column2 = run.read_img(r'***.tif') # 读数据
proj, geotrans, data3, row3, column3 = run.read_img(r'***.tif') # 读数据
# 第二步:将上述读取的3个波段放到一个数组中
data = np.array((data1,data2), dtype=data1.dtype) # 按序将3个波段像元值放入
# data = np.array((data1[0,::],data1[1,::],data2), dtype=data1.dtype) #这个是应对要进行波段组合的图像原本就是多波段的特点,因此我们可以用data[0,::],data[1,::],data[2,::]等表示,前面的0,1,2是图像波段号的索引
# 第三步
run.write_img(r'***.tif', proj, geotrans, data) # 写为3波段数据,假彩色,nir,red,green
多波段.tif数据波段分离
import numpy as np
import os
import random
from osgeo import gdal
# import cv2
import tifffile as tif
from skimage import data,exposure
from sklearn import preprocessing
# 读取图像像素矩阵
# fileName 图像文件名
def readTif(fileName):
dataset = gdal.Open(fileName)
# width = dataset.RasterXSize
# height = dataset.RasterYSize
# GdalImg_data = dataset.ReadAsArray(0, 0, width, height)
# return GdalImg_data
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 bandsclip(path1,path2):
dataset_img=readTif(path1)
width=dataset_img.RasterXSize
height=dataset_img.RasterYSize
proj=dataset_img.GetProjection()
geotrans=dataset_img.GetGeoTransform()
img=dataset_img.ReadAsArray(0,0,width,height)
img_out=[]
#依次将各波段输出
for i in range(img.shape[0]):
img_out=np.array(img[i,::])
#保存tiff格式文件数据
writeTiff(img_out,geotrans,proj,path2+str(i)+'.tif') #输出波段的名称命名格式可以修改,结合传递的path2参数
os.chdir(r'G:/20211216/202109021')
path1=r'202109021.tif' #要分离波段的原始图像数据名称
path2=r'202109021' #分离的各波段结果图像部分名称
bandsclip(path1,path2) #调用上面定义的波段分离函数
print('Bandsclip END!')
更多推荐
已为社区贡献1条内容
所有评论(0)