Landsat8的悬浮泥沙反演实例（By python）

最新推荐文章于 2024-04-07 18:58:53 发布

YukAgame

最新推荐文章于 2024-04-07 18:58:53 发布

阅读量2.3k

点赞数 6

分类专栏：杂记文章标签：计算机视觉可视化 opencv

本文链接：https://blog.csdn.net/weixin_38686737/article/details/116243464

版权

杂记专栏收录该内容

4 篇文章 0 订阅

订阅专栏

根据文献查阅的一次线性反演方程来进行反演
在这里插入图片描述

大概思路就是通过Gdal读取TIF像素并进行计算操作之后再可视化。
网上下载12月22日杭州湾地区的Landsat8图像（冬季泥沙含量高，云量低）

from osgeo import gdal
import os
import numpy as np
import matplotlib.pyplot as plt
import math


class Inversion():
    def __init__(self, type, pos, name):
        self.type = type
        os.chdir(pos)  # 改变默认文件夹位置到对应文件夹
        self.name = name  # 文件名
        # self.run()

    def load_img(self):  # 读入图片
        img = gdal.Open(self.name)
        return img

    def ImgtoArray(self, img, band_num=1):  # 读入图片并将像素值转换成对应矩阵
        data_array = img.GetRasterBand(band_num).ReadAsArray()
        return data_array

    def SSC_INV_LANDSAT8(self, data1, data2, data3):
        '''
        1.此函数是反演悬浮泥沙浓度的
        2.公式为
        2.data1,data2,data3是反演算法当中使用的波段
        '''
        ssc = []
        shape = np.shape(data1)
        for i,v,w in zip(np.reshape(data1,-1),np.reshape(data2,-1),np.reshape(data3,-1)):
            z = v+w
            if z!=0:
                # print(z)
                calc = i/(v+w)
                # print(calc)
                calc = (0.9339*calc-0.3514)*100000
                # calc = Decimal(calc).quantize(0.000000)
#                 calc ='%.6f'%calc
                ssc.append(calc)
            elif z==0:
                ssc.append(14000)
        ssc = np.asarray(ssc).reshape(shape)
        print(ssc.shape)
        return ssc


    def Visualization(self, ssc):  # 可视化经过反演处理过后的图像
        # plt.figure(figsize=(16,9))
        plt.imshow(ssc,cmap='BrBg',vmin = 14000,vmax = 17500)
        plt.legend()
        plt.title('Inversion of %s'%type)
        plt.show()
        # cv2.waitKey(0)


# if __name__ == "main":
os.chdir(r'G:/')
# img = gdal.Open('LC08_L1TP_118039_20201222_20201222_01_RT_B2.TIF')
# data_array = img.GetRasterBand(1).ReadAsArray()


ssc_landsat1 = Inversion('SSC', r'G:/', 'LC08_L1TP_118039_20201222_20201222_01_RT_B2.TIF')
ssc_landsat2 = Inversion('SSC', r'G:/', 'LC08_L1TP_118039_20201222_20201222_01_RT_B3.TIF')
ssc_landsat3 = Inversion('SSC', r'G:/', 'LC08_L1TP_118039_20201222_20201222_01_RT_B4.TIF')
# if ssc_landsat.type =='SSC':
img1 = ssc_landsat1.load_img()
dataset1 = ssc_landsat1.ImgtoArray(img1)
img2 = ssc_landsat2.load_img()
dataset2 = ssc_landsat2.ImgtoArray(img2)
img3 = ssc_landsat3.load_img()
dataset3 = ssc_landsat3.ImgtoArray(img3)
calc = ssc_landsat1.SSC_INV_LANDSAT8(dataset1, dataset2, dataset3)
ssc_landsat1.Visualization(calc)