本文主要是介绍【gdal学习笔记】利用python 的gdal，以及相关库进行遥感图像处理（影像裁剪，辐射定标，大气校正，异常值去除）——以基于landsat8数据提取NDVI为例，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
记录最近的学习
参考博客：参考博客
- AtmosphericCorrection大气校正_landsat8_见贤思齐547的博客-CSDN博客
目录
遥感图像预处理
数据介绍
遥感数字图像存储格式
图像裁剪：
辐射定标：
大气校正
计算NDVI
其他处理函数：
测试
完整代码以及用例
遥感图像预处理
数据介绍
本次实验利用合肥市landsat8影像作归一化植被指数NDVI的提取。
元数据文档：
landsat8数据说明（参考官方使用手册）：
关于landsat8数据说明可以参考此文
https://www.cnblogs.com/icydengyw/p/12056211.html
通过上图可以查看你的landsat数据说明。参阅官方文档，依据自己的数据等级，确定要处理哪些操作，哪些操作不进行处理。
开始着手进行操作：
虽然landsat8有许多波段，但一般我们只用可见光波段，也就是1-7波段。在编写相关函数之前，我们应该先编写存取影像的相关类。
如下：
import numpy as np
import shapefile
import glob
from osgeo import gdal_array, gdal
import re# reader 类
class landsat_reader(object):def __init__(self, path):# 研究区海拔高度self.Altitude = 25self.files_arr = []# 只要7个波段文件 因为在ENVI里也只展示了可见光的七个波段self.bans = 7first_file_name = glob.glob("{}/*.tif".format(path))[0]for i in range(1, self.bans + 1):self.files_arr.append(first_file_name.replace('B1', "B" + str(i)))def index(self, i):return self.files_arr[i]# format_name can choose "GTiff"or "ENVI"def mul_com_fun(self, indexs, out_name, format_name="GTiff", mask=1):if len(indexs) > 2:print("波段合成中")data = read_img(self.band(1))image = np.zeros((data[3], data[4], len(indexs)), dtype='f4')for i in tqdm(range(len(indexs))):data2 = read_img(self.band(indexs[i]))image[:, :, i] = data2[9]del data2# 替换所有为零处的数据image = np.choose(image == 0, (image, np.nan))# 标记代表要不要进行辐射定标与大气校正if mask == 1:image = self.rad_cai(image, indexs)image = self.Atmospheric_correction(image, indexs)write_img(out_name, image, data[1], data[2])del imagedel datadef band(self, i):return self.files_arr[i - 1]def indexs(self):return self.files_arrdef print_filename(self):for f in self.files_arr:print(f)def read_img(self, fileindex):read_img(self.files_arr[fileindex])def rad_cai(self, image, indexs):...def py6s(self, index):...# 进入的是一个np数组def Atmospheric_correction(self, image, indexs):...
读取图像数据的函数：
# 输入tif文件名
# return im_data, im_proj, im_geotrans, im_height, im_width, im_bands
def read_img(dataset):dataset = gdal.Open(dataset) im_width = dataset.RasterXSize im_height = dataset.RasterYSize im_bands = dataset.RasterCount im_geotrans = dataset.GetGeoTransform()im_proj = dataset.GetProjection() im_data = dataset.ReadAsArray() # 下面这个读取的是图像第一波段的的矩阵，是个二维矩阵im_data_2 = dataset.GetRasterBand(1).ReadAsArray()im_dy = dataset.GetRasterBand(1).DataTypex_ize = dataset.GetRasterBand(1).XSizey_ize = dataset.GetRasterBand(1).YSizedel dataset return im_data, im_proj, im_geotrans, im_height, im_width, im_bands, im_dy, x_ize, y_ize, im_data_2由数组生成图像的保存函数如下：这里需要注意的是，gdal库的数组与np的数组形状不一样。gdal数组的形状是深度在前，行列数在后，而np类型的数组形状波段数在行列之后
# 要输入的是数组化的图像
# 生成指定路径的tif图像
def write_img(output, clip, img_prj, img_trans, format_name):# 从该区域读出波段数目，区域大小# 根据不同数组的维度进行读取Is_GDAL_array = clip.shape[0] < 10if len(clip.shape) <= 2:im_bands, (im_height, im_width) = 1, clip.shapeelse:if Is_GDAL_array:im_bands, im_height, im_width = clip.shapeelse:im_height, im_width, im_bands = clip.shapeprint(clip.shape)gtif_driver = gdal.GetDriverByName(format_name)if 'int8' in clip.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in clip.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# （第四个参数代表波段数目，最后一个参数为数据类型，跟原文件一致）output_tif = gtif_driver.Create(output, im_width, im_height, im_bands, datatype)output_tif.SetGeoTransform(img_trans)output_tif.SetProjection(img_prj)print("开始写入")if im_bands == 1:output_tif.GetRasterBand(1).WriteArray(clip)else:for i in range(im_bands):if Is_GDAL_array:output_tif.GetRasterBand(i + 1).WriteArray(clip[i])else:output_tif.GetRasterBand(i + 1).WriteArray(clip[:, :, i])print(i)# 写入磁盘output_tif.FlushCache()# 计算统计信息for i in range(1, im_bands):output_tif.GetRasterBand(i).ComputeStatistics(False)# 创建概视图或金字塔图层output_tif.BuildOverviews('average', [2, 4, 8, 16, 32])# 关闭output文件del output_tifprint("成功写入" + output)对于灰阶图像而言，也就是我们所说的单波段影像，它的存储形式一个二维数组，也就是一个平面 格网上，每一个格网储存着一个像元值。举个例子：
x = np.arange(9.).reshape(3, 3)
print(x)
[[0. 1. 2.][3. 4. 5.][6. 7. 8.]]
# 获取第一行第二个数据可以作如下索引
x[0,1] 数组的行列从索引0开始
1.0对于多波段图像，只要波段大于二，就采用三维数组存储，波段的数目可以叫做数组的深度。可以理解为多张平面矩阵叠加在一起形成了一个三维空间，我们把深度当作z坐标，xy坐标和原来一样不变。下图用一张三通道图像示意图为例，比如说我们要获取多波段图像上的某一个值，应该这么索引
比如x是一个多波段图像的数组，获取第i行第j列红光波段的dn值，可以这么表示：
（对于一般图像转换成的数组，都是第三位坐标为深度）
x[i, j, 0]
#假设红光波段在第一波段，从这个例子可以看出来，深度的索引也是从0开始的在遥感图像处理中我们利用GDAL库来将.tif 格式的遥感数组转换为数组，
dataset = gdal.Open("shiyan.tif") # 打开文件# ReadAsArray()有参数，如果缺省的话，就返回整个范围的数组，不指定波段的话有可能返回三维数组
im_data = dataset.ReadAsArray() # 读取栅格图像的像元数组# 下面这个读取的是图像第一波段的的矩阵，是个二维矩阵
# GetRasterBand(1) ，波段的索引从一开始
im_data_2 = dataset.GetRasterBand(1).ReadAsArray()而对于写入的库，其方法和读取的库的类似，便不多介绍。
gdal库默认保存的图像是按照像素排列的BIP， 通过这行代码改变形式，改变为BSQ存储格式的
output_tif = gtif_driver.Create(output,
im_width, im_height, im_bands,
datatype, options=["INTERLEAVE=BAND"])遥感数字图像存储格式
BSQ（波段顺序格式）
BSQ（band sequential format）是按波段保存，也就是一个波段保存后接着保存第二个波段。该格式最适于对单个波谱波段中任何部分的空间（X，Y）存取。
BIL（波段按行交叉格式）
BIL（band interleaved by line format）是按行保存，就是保存第一个波段的第一行后接着保存第二个波段的第一行，依次类推。该格式提供了空间和波谱处理之间一种折衷方式。
BIP（波段按像元交叉格式）
BIP（band interleaved by pixel format）是按像元保存，即先保存第一个波段的第一个像元，之后保存第二波段的第一个像元，依次保存。该格式为图像数据波谱（Z） 的存取提供最佳性能。
我们利用GDAL 库读取的数组形状与使用numpy创建的数组或者其他库创建的数组并不相同：
GDAL数组形状为（波段数，行数，列数），
而numpy数组为（行数，列数，深度）
下面这行代码展示一个（栅格行数，波段数，栅格列数的）
a = np.array(range(27)).reshape(3, 3, 3)
print(a)
"""
[[[ 0 1 2] 第一波段 的第一行 三个数分别为从左往右 第一个像元值,,,,第二个像元值[ 3 4 5] 第二波段 的第一行 第一个像元值,,,,[ 6 7 8]] 第三波段 的第一行 第一个像元值,,,,[[ 9 10 11] 第一波段 的第二行 第一个像元值,,,,[12 13 14] 第二波段 的第二行 第一个像元值,,,,[15 16 17]] 第三波段 的第二行 第一个像元值,,,,[[18 19 20] 第一波段 的第三行 第一个像元值,,,,[21 22 23] 第二波段 的第三行 第一个像元值,,,,[24 25 26]]] 第三波段 的第三行 第一个像元值,,,,
"""
# 提取第二波段的像元矩阵
print(a[:, 1, :]) # 栅格行数，波段数，栅格列数
"""
[[ 3 4 5][12 13 14][21 22 23]]
"""图像裁剪：
原理就是读取一个范围内的子图层的数组，再重新保存就好了，这是对于利用.shp文件边界裁剪的情况，如果依据形状裁剪，还需要将不在边界内的像元值变为-999.0或者nan（not a number）。代码展示依据边界的裁剪：
# 三个参数依次为 输入的tif路径 ， 输出路径 ， 用于裁剪的shp文件路径
def clip_function(input, output, shp):input_tif = gdal.Open(input)# 读取裁剪的shape文件的外接矩形大小shp = shapefile.Reader(shp)minX, minY, maxX, maxY = shp.bbox# 定义切图的起始点和终点坐标(相比原点的横坐标和纵坐标偏移量)offset_x, offset_y = geotoimagexy(input_tif, minX, minY)endset_x, endset_y = geotoimagexy(input_tif, maxX, maxY)block_xsize = int(endset_x - offset_x + 1)block_ysize = int(abs(endset_y - offset_y) + 1)# 需要注意的是图像原点在左上角，因此我们要将原点的y向上移动，因为y轴向下，所以符号为减offset_y = offset_y - block_ysizeclip = input_tif.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)img_prj = input_tif.GetProjection()img_trans = transform(input_tif.GetGeoTransform(), offset_x, offset_y)# 调用保存函数print("开始裁剪")write_img(output, clip, img_prj, img_trans, "GTiff")辐射定标：
辐射定标是根据一定的参数将遥感的图像的DN值，就是遥感图像的数据量化值，转换为具有物理意义的值
- 转换为绝对辐射亮度值（辐射率）（对应ENVI辐射定标的定标类型（Calibration Type）：辐射率数据Radiance）
- 或者转换与地表（表现）反射率（对应ENVI辐射定标的定标类型（Calibration Type）：第二个参数Reflectance）、表面温度等物理量有关的相对值（Temperature
我们将转化为绝对辐射亮度值，因为公式简单好记，公式为L＝gain*DN+bias ，而且后续的6s模型也需要输入辐亮度。
我们查看landsat8使用手册：
图中给出了每个波段的参数如何获取，可以通过正则表达式在它的元数据文档中获取对应波长参数
第二种，大气反射率
以及辐射亮温：
以定标为辐亮度的方法为例：
代码如下：（这里需要注意的是元数据是16位整型，但是我们要转为话32位浮点型处理，也就是”f4“，在np里代表Float32）
def rad_cai(self, image, indexs):print("辐射定标")# 读取'MTL.txt'内容with open(self.band(1).replace('B1.TIF', 'MTL.txt')) as f:content = f.read()# 利用正则表达式匹配参数gain = re.findall(r'RADIANCE_MULT_BAND_\d\s=\s(\S*)', content)bias = re.findall(r'RADIANCE_ADD_BAND_\d\s=\s(\S*)', content)new_image = np.zeros_like(image, dtype=np.float32)for i in range(len(indexs)):new_image[:, :, i] = float(gain[indexs[i]-1]) * image[:, :, i] + float(bias[indexs[i]-1])return new_image
大气校正
回忆我们在ENVI里进行大气校正时，我们使用的是BIL格式的遥感图像，因为很多教程都说要转为为BIL格式。BIL和BIP格式转换也很简单，新建一个同型数组，在把对应波段储存进去。
有许多的公共大气校正模型，其中6s的效果较好。再python中，可以利用Py6s这个库来进行大气校正，不过这个库安装有些麻烦，然而用
Anaconda 安装十分简单。
conda install gdal
conda install -c conda-forge py6s
代码如下：
def py6s(self, index):# 一些要手动输入的参数# avg海拔高度 单位为千米self.Altitude = 0.030# 气溶胶模型# NoAerosols = 0Continental = 1Maritime = 2Urban = 3 Desert = 5BiomassBurning = 6Stratospheric = 7Aerosol_Model = 3# 设置 50nm气溶胶光学厚度 从这个网站查找 https://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_newaot550 = 0.271# 添加 py6s 预定义的wavelength = Add_wavelength()# 打开landsat8元数据文档with open(self.band(1).replace('B1.TIF', 'MTL.txt')) as f:content = f.read()# 初始化模型，寻找可执行的exe文件s = SixS()s.geometry = Geometry.User()# 设置太阳天顶角和方位角solar_z = re.findall(r'SUN_ELEVATION = (\S*)', content)solar_a = re.findall(r'SUN_AZIMUTH = (\S*)', content)s.geometry.solar_z = 90 - float(solar_z[0])s.geometry.solar_a = float(solar_a[0])# 卫星天顶角和方位角s.geometry.view_z = 0s.geometry.view_a = 0# 获取 影像 范围b_lat = re.findall(r'CORNER_\w*LAT\w*\s=\s(\S*)', content)b_lon = re.findall(r'CORNER_\w*LON\w*\s=\s(\S*)', content)# 求取影像中心经纬度center_lon = np.mean([float(i) for i in b_lon])center_lat = np.mean([float(i) for i in b_lat])# print(center_lon)# print(center_lat)# 正则匹配时间,返回月日time = re.findall(r'DATE_ACQUIRED = (\d{4})-(\d\d)-(\d\d)', content)# print("成像时间是{}年{}月{}日".format(time[0][0], time[0][1], time[0][2]))s.geometry.month = int(time[0][1])s.geometry.day = int(time[0][2])# 大气模式类型if -15 < center_lat <= 15:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)if 15 < center_lat <= 45:if 4 < s.geometry.month <= 9:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)else:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)if 45 < center_lat <= 60:if 4 < s.geometry.month <= 9:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)else:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)s.aero_profile = AtmosProfile.PredefinedType(Aerosol_Model)# 这个参数不是很明白s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)# ENVI中这个默认就是40，单位千米 # s.visibility = 40.0# 550nm气溶胶光学厚度s.aot550 = aot550# 研究区海拔、卫星传感器轨道高度s.altitudes = Altitudes()s.altitudes.set_target_custom_altitude(self.Altitude)# 将传感器高度设置为卫星高度s.altitudes.set_sensor_satellite_level()"""PredefinedWavelengths.LANDSAT_OLI_B1预定义的波长，根据点后面的关键字查找 ，py6s库里面列出了B1 到 B9的波长"""# 设置b波普响应函数s.wavelength = Wavelength(wavelength[index])# 下垫面非均一、朗伯体s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)# 运行6s大气模型s.run()xa = s.outputs.coef_xaxb = s.outputs.coef_xbxc = s.outputs.coef_xcx = s.outputs.valuesreturn xa, xb, xc# 进入的是一个np数组def Atmospheric_correction(self, image, indexs):print("大气校正开始")new_image = np.zeros_like(image, dtype=np.float32)for i in tqdm(range(len(indexs))):a, b, c = self.py6s(indexs[i])x = image[:, :, i]y = a * x - bnew_image[:, :, i] = y / (1 + y * c) * 10000return new_image经过对比，发现与ENVI做的同一幅大气校正遥感影像相比，曲线走向类似，不过不同波段对应的最大值有些不同，可能是因为参数不同导致。顺带一提，辐射定标系数用的是1.0
进行下一步：
计算NDVI
开始着手进行操作：
由上图我们可以得知红光波段是band4 ，近红外波段是band5 ，因此得出
计算公式应该为
计算NDVI之后，我们还需要剔除值小于0和大于1的情况
def getndvi(nir_data, red_data):try:denominator = np.array(nir_data + red_data, dtype='f4')numerator = np.array(nir_data - red_data, dtype='f4')ndvi = np.divide(numerator, denominator, where=denominator != 0.0)# 去除异常值 使得ndvi的值在0与1之间ndvi = np.choose((ndvi < 0) + 2 * (ndvi > 1), (ndvi, 0, 0))return ndviexcept BaseException as e:print(str(e))其他处理函数：
线性拉伸，单波段阈值法分类，建立颜色卡，查看灰度直方图之类的
# 此函数进入的是gdal类型数组
# 包括多波段与单波段
def linear_stretch(gray, truncated_value, max_out=255, min_out=0):truncated_down = np.percentile(gray, truncated_value)truncated_up = np.percentile(gray, 100 - truncated_value)gray = (gray - truncated_down) / (truncated_up - truncated_down) * (max_out - min_out) + min_outgray = np.choose(gray < min_out + 2 * (gray > max_out), (gray, min_out, max_out))gray = np.uint8(gray)return gray# 传入tif文件
def show_hist(src):# 直方图src_array = read_img(src)[0]# 灰度拉伸src_array = linear_stretch(src_array, 2)plt.hist(src_array)plt.title("Single_band_histogram of %s" % src)plt.xlabel("x axis caption")plt.ylabel("y axis caption")print("Histogram display")plt.show()# Image classification by single band threshold method
# Input parameter TIF file, target classification picture, threshold
def classification(src, tgt, threshold):print("Image classification start")threshold = list(threshold)src_data = read_img(src)src_array = np.array(src_data[0])# 灰度拉伸src_array = linear_stretch(src_array, 2)Max = np.max(src_array)Min = np.min(src_array)print("The maximum value of data is%s" % Max)print("The minimum value of data is%s" % Min)print(src_array)# 保存提取结果rgb = Create_color_slice(src_array)# color = [0, 0, 175]# # 掩膜# mask = gdal_array.numpy.logical_and(src_array > threshold[0], src_array < threshold[1])# rgb = gdal_array.numpy.zeros((3, src_array.shape[0], src_array.shape[1],), gdal_array.numpy.uint8)# for i in range(3):# rgb[i] = np.choose(mask, (255, color[i]))output = gdal_array.SaveArray(rgb.astype(gdal_array.numpy.uint8), tgt, format="JPEG")del outputdel src_data# 输入np数组
def np_conversion_gdal(array):if len(array.shape) == 2:gdal_arr = arrayelse:gdal_arr = gdal_array.numpy.zeros((array.shape[2], array.shape[0], array.shape[1],), gdal_array.numpy.uint8)for i in range(array.shape[2]):gdal_arr[i] = array[:, :, i]return gdal_arr# 输入gdal类型数组
def gdal_conversion_np(array):array2 = array.shapenp_arr = np.zeros(shape=(array2[1], array2[2], array2[0]))for i in range(array2[0]):np_arr[:, :, i] = array[i]return np_arr# 输入数组 返回一个gdal类型的数组 ,输入的是灰度直方图
def Create_color_slice(array, bins=10):# 获取区间# random.randint(0, 255)classes = gdal_array.numpy.histogram(array, bins=bins)[1]rgb = gdal_array.numpy.zeros((3, array.shape[0], array.shape[1],), gdal_array.numpy.uint8)for i in range(bins):mask = gdal_array.numpy.logical_and(array > classes[i], array < classes[i + 1])for j in tqdm(range(3)):rgb[j] = np.choose(mask, (rgb[j], random.randint(0, 255)))rgb = rgb.astype(gdal_array.numpy.uint8)return rgb
测试
最后我们实际测试一下效果：
# 主函数
if __name__ == "__main__":print("主函数开始运行")# 读取路径下的landsat文件夹file = landsat_reader("./LC81210382021028LGN00")# 先对第4波段和第5波段进行辐射校正和大气校正file.mul_com_fun([4], "red.tif")file.mul_com_fun([5], "nir.tif")# 根据裁剪文件进行裁剪clip_function("red.tif", "landsat_red.tif", "裁剪用数据/utm50n_mask.shp")clip_function("nir.tif", "landsat_nir.tif", "裁剪用数据/utm50n_mask.shp")# 读取裁剪后文件red_signal = read_img("landsat_red.tif")nir_signal = read_img("landsat_nir.tif")# 计算植被指数ndvi = getndvi(nir_signal[0], red_signal[0])write_img("ndvi.tif", ndvi, red_signal[1], red_signal[2])运行完美：查看效果：
完整代码以及用例
import glob
import random
import re
import numpy as np
import shapefile
from Py6S import *
from matplotlib import pyplot as plt
from osgeo import gdal_array, gdal
from tqdm import tqdm# 这个文件的作用是进行栅格读写 与其他功能
# 遥感影像下载接口 ： 地理空间数据云 http://www.gscloud.cn/
# 国外接口： EarthExplorer https://earthexplorer.usgs.gov/# reader 类
class landsat_reader(object):# 提供landsat8数据所在目录def __init__(self, path):# 研究区海拔高度self.Altitude = 25self.files_arr = []# 只要7个波段文件 因为在ENVI里也只展示了可见光的七个波段self.bans = 7first_file_name = glob.glob("{}/*.tif".format(path))[0]for i in range(1, self.bans + 1):self.files_arr.append(first_file_name.replace('B1', "B" + str(i)))# 返回文件索引def index(self, i):return self.files_arr[i]# format_name can choose "GTiff"or "ENVI"def mul_com_fun(self, indexs, out_name, format_name="GTiff", mask=1):if len(indexs) > 2:print("波段合成中")data = read_img(self.band(1))image = np.zeros((data[3], data[4], len(indexs)), dtype='f4')for i in tqdm(range(len(indexs))):data2 = read_img(self.band(indexs[i]))image[:, :, i] = data2[0]# 替换所有为零处的数据image = np.choose(image == 0, (image, np.nan))# 标记代表要不要进行辐射定标if mask == 1:image = self.rad_cai(image, indexs)image = self.Atmospheric_correction(image, indexs)write_img(out_name, image, data[1], data[2])del imagedel data# 用波段号作为索引def band(self, i):return self.files_arr[i - 1]# 返回整个索引def indexs(self):return self.files_arr# 打印文件名def print_filename(self):for f in self.files_arr:print(f)# 读取数组def read_img(self, fileindex):read_img(self.files_arr[fileindex])# 辐射定标def rad_cai(self, image, indexs):print("辐射定标")# 读取'MTL.txt'内容with open(self.band(1).replace('B1.TIF', 'MTL.txt')) as f:content = f.read()# 利用正则表达式匹配参数gain = re.findall(r'RADIANCE_MULT_BAND_\d\s=\s(\S*)', content)bias = re.findall(r'RADIANCE_ADD_BAND_\d\s=\s(\S*)', content)new_image = np.zeros_like(image, dtype=np.float32)for i in tqdm(range(len(indexs))):new_image[:, :, i] = (float(gain[indexs[i] - 1]) * image[:, :, i] + float(bias[indexs[i] - 1]))return new_image# 6s模型参数获取def py6s(self, index):# 一些要手动输入的参数# avg海拔高度 单位为千米self.Altitude = 0.030# 气溶胶模型# NoAerosols = 0Continental = 1Maritime = 2Urban = 3 Desert = 5BiomassBurning = 6Stratospheric = 7Aerosol_Model = 3# 设置 50nm气溶胶光学厚度 从这个网站查找 https://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_newaot550 = 0.271# 添加 py6s 预定义的wavelength = Add_wavelength()# 打开landsat8元数据文档with open(self.band(1).replace('B1.TIF', 'MTL.txt')) as f:content = f.read()# 初始化模型，寻找可执行的exe文件s = SixS()s.geometry = Geometry.User()# 设置太阳天顶角和方位角solar_z = re.findall(r'SUN_ELEVATION = (\S*)', content)solar_a = re.findall(r'SUN_AZIMUTH = (\S*)', content)s.geometry.solar_z = 90 - float(solar_z[0])s.geometry.solar_a = float(solar_a[0])# 卫星天顶角和方位角s.geometry.view_z = 0s.geometry.view_a = 0# 获取 影像 范围b_lat = re.findall(r'CORNER_\w*LAT\w*\s=\s(\S*)', content)b_lon = re.findall(r'CORNER_\w*LON\w*\s=\s(\S*)', content)# 求取影像中心经纬度center_lon = np.mean([float(i) for i in b_lon])center_lat = np.mean([float(i) for i in b_lat])# print(center_lon)# print(center_lat)# 正则匹配时间,返回月日time = re.findall(r'DATE_ACQUIRED = (\d{4})-(\d\d)-(\d\d)', content)# print("成像时间是{}年{}月{}日".format(time[0][0], time[0][1], time[0][2]))s.geometry.month = int(time[0][1])s.geometry.day = int(time[0][2])# 大气模式类型if -15 < center_lat <= 15:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)if 15 < center_lat <= 45:if 4 < s.geometry.month <= 9:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)else:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)if 45 < center_lat <= 60:if 4 < s.geometry.month <= 9:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)else:s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)s.aero_profile = AtmosProfile.PredefinedType(Aerosol_Model)# 这个参数不是很明白s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)# ENVI中这个默认就是40，单位千米 ，听说Py6s中不能用这个# s.visibility = 40.0# 550nm气溶胶光学厚度s.aot550 = aot550# 研究区海拔、卫星传感器轨道高度s.altitudes = Altitudes()s.altitudes.set_target_custom_altitude(self.Altitude)# 将传感器高度设置为卫星高度 ，非常疑惑，它怎么知道我卫星高度多少s.altitudes.set_sensor_satellite_level()"""PredefinedWavelengths.LANDSAT_OLI_B1预定义的波长，根据点后面的关键字查找 ，py6s库里面列出了B1 到 B9的波长"""# 设置b波普响应函数s.wavelength = Wavelength(wavelength[index])# 下垫面非均一、朗伯体s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)# 运行6s大气模型s.run()xa = s.outputs.coef_xaxb = s.outputs.coef_xbxc = s.outputs.coef_xcx = s.outputs.valuesreturn xa, xb, xc# 大气校正# 进入的是一个np数组def Atmospheric_correction(self, image, indexs):print("大气校正开始")new_image = np.zeros_like(image, dtype=np.float32)for i in tqdm(range(len(indexs))):a, b, c = self.py6s(indexs[i])x = image[:, :, i]y = a * x - bnew_image[:, :, i] = y / (1 + y * c) * 10000return new_imagedef Add_wavelength():wavelength = [0, PredefinedWavelengths.LANDSAT_OLI_B1,PredefinedWavelengths.LANDSAT_OLI_B2,PredefinedWavelengths.LANDSAT_OLI_B3,PredefinedWavelengths.LANDSAT_OLI_B4,PredefinedWavelengths.LANDSAT_OLI_B5,PredefinedWavelengths.LANDSAT_OLI_B6,PredefinedWavelengths.LANDSAT_OLI_B7,PredefinedWavelengths.LANDSAT_OLI_B8,PredefinedWavelengths.LANDSAT_OLI_B9,]return wavelength# 此函数将地图坐标转换为影像坐标，通过gdal的六参数模型
def geotoimagexy(dataset, x, y):trans = dataset.GetGeoTransform()a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])b = np.array([x - trans[0], y - trans[3]])return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解# 读取原图仿射变换参数值
def transform(ori_transform, offset_x, offset_y):# 以列表的形式返回变换数组ori_transform = list(ori_transform)top_left_x = ori_transform[0] + offset_x * ori_transform[1]top_left_y = ori_transform[3] + offset_y * ori_transform[5]# 根据仿射变换参数计算新图的原点坐标ori_transform[0] = top_left_xori_transform[3] = top_left_yreturn ori_transform# 输入tif文件名
# return im_data, im_proj, im_geotrans, im_height, im_width, im_bands
def read_img(dataset):dataset = gdal.Open(dataset) # 打开文件im_width = dataset.RasterXSize # 栅格矩阵的列数im_height = dataset.RasterYSize # 栅格矩阵的行数im_bands = dataset.RasterCount # 波段数'''GDAL仿射矩阵，包含六个参数 im_geotrans[0]为图像左上角的x坐标 ，im_geotrans[3]为左上角的y坐标im_geotrans[1] x方向像素分辨率 ， im_geotrans[5] 为y方向的分辨率im_geotrans[2]和 im_geotrans[4] 都为旋转角度，图像朝上的话一般值为0'''im_geotrans = dataset.GetGeoTransform()im_proj = dataset.GetProjection() # 地图投影信息，字符串表示im_data = dataset.ReadAsArray() # 读取栅格图像的像元数组# 下面这个读取的是图像第一波段的的矩阵，是个二维矩阵im_data_2 = dataset.GetRasterBand(1).ReadAsArray()im_dy = dataset.GetRasterBand(1).DataTypex_ize = dataset.GetRasterBand(1).XSizey_ize = dataset.GetRasterBand(1).YSizedel dataset # 关闭对象dataset，释放内存return im_data, im_proj, im_geotrans, im_height, im_width, im_bands, im_dy, x_ize, y_ize, im_data_2# 要输入的是数组化的图像
# 生成指定路径的tif图像
def write_img(output, clip, img_prj, img_trans, format_name="GTiff"):# 从该区域读出波段数目，区域大小# 根据不同数组的维度进行读取Is_GDAL_array = clip.shape[0] < 10if len(clip.shape) <= 2:im_bands, (im_height, im_width) = 1, clip.shapeelse:if Is_GDAL_array:im_bands, im_height, im_width = clip.shapeelse:im_height, im_width, im_bands = clip.shape# 获取Tif的驱动，为创建切出来的图文件做准备gtif_driver = gdal.GetDriverByName(format_name)if 'int8' in clip.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in clip.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# （第四个参数代表波段数目，最后一个参数为数据类型，跟原文件一致）output_tif = gtif_driver.Create(output, im_width, im_height, im_bands, datatype, options=["INTERLEAVE=BAND"])output_tif.SetGeoTransform(img_trans)output_tif.SetProjection(img_prj)print("开始写入")if im_bands == 1:if len(clip.shape) == 3:clip = clip[:, :, 0]output_tif.GetRasterBand(1).WriteArray(clip)else:for i in tqdm(range(im_bands)):if Is_GDAL_array:output_tif.GetRasterBand(i + 1).WriteArray(clip[i])else:output_tif.GetRasterBand(i + 1).WriteArray(clip[:, :, i])# 写入磁盘output_tif.FlushCache()# 计算统计信息# for i in range(1, im_bands):# output_tif.GetRasterBand(i).ComputeStatistics(False)# # 创建概视图或金字塔图层# output_tif.BuildOverviews('average', [2, 4, 8, 16, 32])# 关闭output文件del output_tifprint("成功写入" + output)# 三个参数依次为 输入的tif路径 ， 输出路径 ， 用于裁剪的shp文件路径
def clip_function(input, output, shp):# 打开需要裁剪的遥感影像input_tif = gdal.Open(input)# 读取裁剪的shape文件的外接矩形大小shp = shapefile.Reader(shp)minX, minY, maxX, maxY = shp.bbox# 定义切图的起始点和终点坐标(相比原点的横坐标和纵坐标偏移量)offset_x, offset_y = geotoimagexy(input_tif, minX, minY)endset_x, endset_y = geotoimagexy(input_tif, maxX, maxY)# 定义切图的大小（矩形框）block_xsize = int(endset_x - offset_x + 1)block_ysize = int(abs(endset_y - offset_y) + 1)# 将裁剪区域的影像转换为数组# 需要注意的是图像原点在左上角，因此我们要将原点的y向上移动，因为y轴向下，所以符号为减offset_y = offset_y - block_ysizeclip = input_tif.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)img_prj = input_tif.GetProjection()img_trans = transform(input_tif.GetGeoTransform(), offset_x, offset_y)# 调用保存函数print("开始裁剪")write_img(output, clip, img_prj, img_trans, "GTiff")# 获取ndvi ，返回数组 ，去除异常值
def getndvi(nir_data, red_data):try:# 将红外数组与近红外数组相加，得到新数组，每一个行列确定的像元上，都会带有一个像素信息，这些就是波段存储的值'''np里面的数组运算需要行列相等，也就是shape相等，但是对于行或者列向量可以豁免，计算时将扩充为同型矩阵 ，矩阵的数学乘法用np.dot()函数矩阵乘法，前行X后列，需要满足前列等于后行'''denominator = np.array(nir_data + red_data, dtype='f4')numerator = np.array(nir_data - red_data, dtype='f4')ndvi = np.divide(numerator, denominator, where=denominator != 0.0)# 去除异常值 使得ndvi的值在0与1之间ndvi = np.choose((ndvi < 0) + 2 * (ndvi > 1), (ndvi, 0, 0))return ndviexcept BaseException as e:print(str(e))# 此函数进入的是gdal类型数组
# 包括多波段与单波段
def linear_stretch(gray, truncated_value, max_out=255, min_out=0):truncated_down = np.percentile(gray, truncated_value)truncated_up = np.percentile(gray, 100 - truncated_value)gray = (gray - truncated_down) / (truncated_up - truncated_down) * (max_out - min_out) + min_outgray = np.choose(gray < min_out + 2 * (gray > max_out), (gray, min_out, max_out))gray = np.uint8(gray)return gray# 传入tif文件
def show_hist(src):# 直方图src_array = read_img(src)[0]# 灰度拉伸if src_array.shape[0] < 10:src_array = linear_stretch(src_array[0, :, :], 2)print(src_array)plt.hist(src_array)plt.title("Single_band_histogram of %s" % src)plt.xlabel("x axis caption")plt.ylabel("y axis caption")print("Histogram display")plt.show()# Image classification by single band threshold method
# Input parameter TIF file, target classification picture, threshold
def classification(src, tgt, threshold):print("Image classification start")threshold = list(threshold)src_data = read_img(src)src_array = np.array(src_data[0])# 灰度拉伸if src_array.shape[0] < 10:src_array = linear_stretch(src_array[0, :, :], 2)Max = np.max(src_array)Min = np.min(src_array)print("The maximum value of data is%s" % Max)print("The minimum value of data is%s" % Min)print(src_array)# 保存提取结果# rgb = Create_color_slice(src_array)color = [0, 0, 175]# 掩膜mask = gdal_array.numpy.logical_and(src_array > threshold[0], src_array < threshold[1])rgb = gdal_array.numpy.zeros((3, src_array.shape[0], src_array.shape[1],), gdal_array.numpy.uint8)for i in range(3):rgb[i] = np.choose(mask, (255, color[i]))output = gdal_array.SaveArray(rgb.astype(gdal_array.numpy.uint8), tgt + ".jpg", format="JPEG")del outputdel src_data# 输入np数组
def np_conversion_gdal(array):if len(array.shape) == 2:gdal_arr = arrayelse:gdal_arr = gdal_array.numpy.zeros((array.shape[2], array.shape[0], array.shape[1],), gdal_array.numpy.uint8)for i in range(array.shape[2]):gdal_arr[i] = array[:, :, i]return gdal_arr# 输入gdal类型数组
def gdal_conversion_np(array):array2 = array.shapenp_arr = np.zeros(shape=(array2[1], array2[2], array2[0]))for i in range(array2[0]):np_arr[:, :, i] = array[i]return np_arr# 输入数组 返回一个gdal类型的数组 ,输入的是灰度直方图
def Create_color_slice(array, bins=10):# 获取区间# random.randint(0, 255)classes = gdal_array.numpy.histogram(array, bins=bins)[1]rgb = gdal_array.numpy.zeros((3, array.shape[0], array.shape[1],), gdal_array.numpy.uint8)for i in range(bins):mask = gdal_array.numpy.logical_and(array > classes[i], array < classes[i + 1])for j in tqdm(range(3)):rgb[j] = np.choose(mask, (rgb[j], random.randint(0, 255)))rgb = rgb.astype(gdal_array.numpy.uint8)return rgb
不过有个问题需要解决，也就是大气校正精度的问题。
# 引入写好的文件
from landsat8_fun import *
from Py6S import SixS# 主函数
if __name__ == "__main__":# 这行代码测试你的6s有没有正常运行# 如果你的py6s 库找不到路径请注释掉 辐射定标和大气校正的代码SixS.test()print("主函数开始运行")# bind_math("(a-b)/(c+d)", 1, 5, 25, 29)# 读取路径下的landsat文件夹file = landsat_reader("./LC81210382021028LGN00")# 先对第4波段和第5波段进行辐射校正和大气校正，第一个参数是一个列表，数字代表landsat第几个波段，一旦长度超过二，就会进行波段融合# 如果你的py6s安装不上注释下述代码file.mul_com_fun([3], "./测试结果/green.tif")file.mul_com_fun([6], "./测试结果/nir.tif")# 根据裁剪文件进行裁剪clip_function("./测试结果/green.tif", "./测试结果/landsat_green.tif", "裁剪用数据/utm50n_mask.shp")clip_function("./测试结果/nir.tif", "./测试结果/landsat_nir.tif", "裁剪用数据/utm50n_mask.shp")# 跳过辐射定标和大气校正# clip_function(file.band(4), "./测试结果/landsat_green.tif", "裁剪用数据/utm50n_mask.shp")# clip_function(file.band(5), "./测试结果/landsat_nir.tif", "裁剪用数据/utm50n_mask.shp")# 读取裁剪后文件green_signal = read_img("./测试结果/landsat_green.tif")nir_signal = read_img("./测试结果/landsat_nir.tif")# 计算植被指数ndsi = getndvi(green_signal[0], nir_signal[0])write_img("./测试结果/ndsi.tif", ndsi, green_signal[1], green_signal[2])这篇关于【gdal学习笔记】利用python 的gdal，以及相关库进行遥感图像处理（影像裁剪，辐射定标，大气校正，异常值去除）——以基于landsat8数据提取NDVI为例的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
