gdal读取了遥感数据后就变成了numpy数组了,之后就是numpy的操作了,numpy操作了之后再加坐标保存就行了,实际上很简单读取和写入看下面的连接就够了
(4条消息) 5. gdal实现对遥感影像的读写,信息统计和波段选择_xiaotiig的博客-CSDN博客
https://blog.csdn.net/xiaotiig/article/details/118721427
1 提取4波段的前3波段
# 将4波段的遥感影像提取出前3波段
import os
import gdal
import numpy as np
# 读取tif数据集
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")
# 这个写入函数比较特殊,跟数据的shape不一样,不一样才对,数据对应行列,这里必须是列行,先写列,再写行
# https://vimsky.com/examples/detail/python-method-gdal.GetDriverByName.html
dataset = driver.Create(path, int(im_width), int(im_height), int(3), datatype)
if (dataset != None):
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(im_proj) # 写入投影
# 这里写入前3个波段,全部输出im_bands
for i in range(3):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
def bands4to3(TifPath, SavePath):
i = 0
filenames = os.listdir(TifPath)
print(len(filenames))
for filename in filenames:
i +=1
if i%1000==0:
print("number:%d"%i)
# 输出文件名
file_path = os.path.join(TifPath, filename)
# 读取栅格文件
dataset_img = readTif(file_path)
proj = dataset_img.GetProjection()
geotrans = dataset_img.GetGeoTransform()
img = dataset_img.ReadAsArray() # 获取数据
writeTiff(img, geotrans, proj, os.path.join(SavePath, filename))
if __name__ == '__main__':
bands4to3(r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgislLuotian\3imageDataset\1lt497908\images",
r"H:\01HTutorWork\3GF2\2DataAndLabel\4Gaofen2\2arcgislLuotian\3imageDataset\1lt497908\images3"
)
2 进行波段的组合计算
在4波段的基础上加入ndvi波段,如果之后要保存为tif文件,看
(4条消息) 9.将数组存储为带坐标的tif文件_xiaotiig的博客-CSDN博客
https://blog.csdn.net/xiaotiig/article/details/122308138
def imgread(fileName, addNDVI=True):
dataset = gdal.Open(fileName)
width = dataset.RasterXSize
height = dataset.RasterYSize
data = dataset.ReadAsArray(0, 0, width, height) # data, (4, 36786, 37239) ,波段,行,列
# 如果是image的话,才进行处理,因为label是单通道,
if(len(data.shape) == 3):
# 添加归一化植被指数NDVI特征
if(addNDVI):
# 这景影像的波段顺序为,R,G,B,NIR
nir, r = data[3], data[0]
ndvi = (nir - r) / (nir + r + 0.00001) * 1.0 # 这里已经计算好了ndvi了
# 和其他波段保持统一,归到0-255,在深度学习中后面的totensor会/255统一归一化
# -1到1,先变成0到2,再乘
ndvi = (ndvi+1) * 127.5
# 去除小于0和大于255的值,这步不需要
# ndvi = np.clip(ndvi, 0, 255)
data_add_ndvi = np.zeros((5, 256, 256), np.uint8)
data_add_ndvi[0:4] = data
# 先进行四舍五入,再转换为8位深度
data_add_ndvi[4] = np.uint8(np.around(ndvi))
data = data_add_ndvi
# (C,H,W)->(H,W,C) 这里进行了波段转换,为了深度学习处理,如果单纯的要保存,不用转换用(C,H,W)的数组
data = data.swapaxes(1, 0).swapaxes(1, 2)
return data