记录最近的学习
参考博客：参考博客
目录
遥感图像预处理
数据介绍
本次实验利用合肥市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 = 25
self.files_arr = []
# 只要7个波段文件 因为在ENVI里也只展示了可见光的七个波段
self.bans = 7
first_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 image
del 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):
...
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).DataType
x_ize = dataset.GetRasterBand(1).XSize
y_ize = dataset.GetRasterBand(1).YSize
del 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] < 10
if len(clip.shape) <= 2: im_bands, (im_height, im_width)="1," clip.shape else: if is_gdal_array: im_height, im_width="clip.shape" im_width, im_bands="clip.shape" print(clip.shape) gtif_driver="gdal.GetDriverByName(format_name)" 'int8' in clip.dtype.name: datatype="gdal.GDT_Byte" elif 'int16' # （第四个参数代表波段数目，最后一个参数为数据类型，跟原文件一致） output_tif="gtif_driver.Create(output," datatype) output_tif.setgeotransform(img_trans) output_tif.setprojection(img_prj) print("开始写入") 1: output_tif.getrasterband(1).writearray(clip) for i range(im_bands): output_tif.getrasterband(i + 1).writearray(clip[i]) 1).writearray(clip[:, :, i]) print(i) 写入磁盘 output_tif.flushcache() 计算统计信息 range(1, im_bands): output_tif.getrasterband(i).computestatistics(false) 创建概视图或金字塔图层 output_tif.buildoverviews('average', [2, 4, 8, 16, 32]) 关闭output文件 del print("成功写入" output)< code></=>
对于灰阶图像而言，也就是我们所说的单波段影像，它的存储形式一个二维数组，也就是一个平面 格网上，每一个格网储存着一个像元值。举个例子：
x = np.arange(9.).reshape(3, 3)
print(x)
[[0. 1. 2.]
[3. 4. 5.]
[6. 7. 8.]]
获取第一行第二个数据可以作如下索引
x[0,1] 数组的行列从索引0开始
1.0
比如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_ysize
clip = 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做的，右边使用上述代码做的，基本差别不大
大气校正
回忆我们在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 = 7
Aerosol_Model = 3
# 设置 50nm气溶胶光学厚度 从这个网站查找 https://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_new
aot550 = 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 = 0
s.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
经过对比，发现与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 ndvi
except 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_out
gray = 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 output
del src_data
输入np数组
def np_conversion_gdal(array):
if len(array.shape) == 2:
gdal_arr = array
else:
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.shape
np_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 = 25
self.files_arr = []
# 只要7个波段文件 因为在ENVI里也只展示了可见光的七个波段
self.bans = 7
first_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 image
del 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 = 7
Aerosol_Model = 3
# 设置 50nm气溶胶光学厚度 从这个网站查找 https://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_new
aot550 = 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 = 0
s.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 1), (ndvi, 0, 0))
return ndvi
except 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_out
gray = 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 output
del src_data
输入np数组
def np_conversion_gdal(array):
if len(array.shape) == 2:
gdal_arr = array
else:
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.shape
np_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])
Original: https://blog.csdn.net/qq_47730141/article/details/125560605
Author: 舒南风
Title: 【gdal学习笔记】利用python 的gdal，以及相关库进行遥感图像处理（影像裁剪，辐射定标，大气校正，异常值去除）——以基于landsat8数据提取NDVI为例
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/757199/
转载文章受原作者版权保护。转载请注明原作者出处！