多角度偏振遥感云相态分类(python实现)

多角度遥感云相态分类做完了，毕业设计的编程实现是利用python编程语言，实现了遥感中云相态分类，代码比较粗糙，不足之处，还望大家多多指教。下面直接上代码：

from osgeo import gdal
import os
import numpy as np
import math
from time import time

im_width = 6480    #栅格矩阵的列数
im_height = 3240  #栅格矩阵的行数

class GRID:
    #读图像文件
    def read_img(self,filename):
        dataset=gdal.Open(filename)       #打开文件
        im_geotrans = dataset.GetGeoTransform()  #仿射矩阵
        ds = gdal.Open('data.h5')
        im_data = ds.GetMetadata()
        im_proj = im_data.get('Geographic_Projection')  # 地图投影信息
        im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
        del dataset
        return im_proj,im_geotrans,im_data
    #将三个角度文件读为数组，并返回值
    def angle(self,filename):
        gdal.PushErrorHandler('CPLQuietErrorHandler')
        dataset = gdal.Open(filename)
        im_angle = dataset.GetRasterBand(1).ReadAsArray(0,0,im_width,im_height)
        del dataset
        return im_angle
    #计算COS角度值
    def cos_angle(self,a):
        '''
        传入二维数组，然后求取COS之后再
        return计算的值
        '''
        X, Y = 3240, 6480  # 宽高
        # ===================准备参数开始==================
        start = time()
        r, pi2 = 0.0015 * math.pi / 180, 2 * math.pi  # 常量
        temp, tx = [], 0  # 临时变量
        while tx < pi2:
            temp.append(math.cos(tx))
            tx += r
        # 常量
        tempLen = len(temp)
        end = time()
        print("准备参数耗时:", (end - start) * 1000, "毫秒")

        # ===================准备参数结束==================
        tempA = a.tolist()
        y = len(tempA)
        if y == 0: return [[]]
        x = len(tempA[0])
        b = [[0 for i in range(x)] for j in range(y)]
        for m in range(y):
            for n in range(x):
                if tempA[m][n] == 65535:
                    b[m][n] = 65535
                else:
                    b[m][n] = temp[tempA[m][n] % tempLen]
        return np.array(b)
    #计算SIN角度值
    def sin_angle(self,a):
        '''
        return计算的值
        '''
        X, Y = 3240, 6480  # 宽高
        # ===================准备参数开始==================
        start = time()
        r, pi2 = 0.0015 * math.pi / 180, 2 * math.pi  # 常量
        temp, tx = [], 0  # 临时变量
        while tx < pi2:
            temp.append(math.sin(tx))
            tx += r
        # 常量
        tempLen = len(temp)
        end = time()
        print("准备参数耗时:", (end - start) * 1000, "毫秒")

        # ===================准备参数结束==================
        tempA = a.tolist()
        y = len(tempA)
        if y == 0: return [[]]
        x = len(tempA[0])
        b = [[0 for i in range(x)] for j in range(y)]
        for m in range(y):
            for n in range(x):
                if tempA[m][n] == 65535:
                    b[m][n] = 65535
                else:
                    b[m][n] = temp[tempA[m][n] % tempLen]
        return np.array(b)
    def cos_r_a_angle(self,a):
        X, Y = 3240, 6480  # 宽高
        # ===================准备参数开始==================
        start = time()
        r, pi2 = 0.006 * math.pi / 180, 2 * math.pi  # 常量
        temp, tx = [], 0  # 临时变量
        while tx < pi2:
            temp.append(math.cos(tx))
            tx += r
        # 常量
        tempLen = len(temp)
        end = time()
        print("准备参数耗时:", (end - start) * 1000, "毫秒")

        # ===================准备参数结束==================
        tempA = a.tolist()#矩阵转列表
        y = len(tempA)#计算多少行
        if y == 0: return [[]]
        x = len(tempA[0])
        b = [[0 for i in range(x)] for j in range(y)]#声明b数组
        for m in range(y):
            for n in range(x):
                if tempA[m][n] == 65535:
                    b[m][n] = 65535
                else:
                    b[m][n] = temp[tempA[m][n] % tempLen]
        return np.array(b)
    #计算散射角
    def san_angle_1(self,cos_theta_zero, cos_theta):
        cos_multipy = np.zeros(shape=(3240, 6480))
        for m in range(3240):#这个可以改为3240，后期再说叭，先搞好大头的东西
            for n in range(6480):#这个可以改为6480
                if cos_theta_zero[m][n] == 65535 or cos_theta[m][n] == 65535:
                    cos_multipy[m][n] = 65535
                else:
                    cos_multipy[m][n] = -cos_theta_zero[m][n] * cos_theta[m][n]
        return cos_multipy
    def san_angle_2(self,sin_theta_zero,sin_theta):
        result=np.zeros(shape=(3240,6480))
        for m in range(3240):
            for n in range(6480):
                if sin_theta_zero[m][n] == 65535 or sin_theta[m][n] == 65535:
                    result[m][n] = 65535
                else:
                    result[m][n] = sin_theta_zero[m][n] * sin_theta[m][n]
        return result
    def san_angle_3(self,sin_rusult,cos_fai):
        result=np.zeros(shape=(3240,6480))
        for m in range(3240):
            for n in range(6480):
                if sin_rusult[m][n] == 65535 or cos_fai[m][n] == 65535:
                    result[m][n] = 65535
                else:
                    result[m][n] = sin_rusult[m][n] * cos_fai[m][n]
        return result
    def acos_san_angle(self,a):
        b=np.zeros(shape=(3240,6480))
        for m in range(3240):
            for n in range(6480):
                if 0< a[m][n] <math.pi:
                    b[m][n] = math.acos(a[m][n])#* math.pi/180)
        return np.array(b)

    def cunzai(self, san_she_angle, cloud_result,PR):
        # 散射角范围a:135-143.23 b:125-135 点得个数为3935
        # ********a:135-145    b:122.508-135 点的个数为：4915 # 实验采用第二个散射角范围(点的个数多，较为直观的显示)
        temp1 = []
        temp2 = []
        line_m = []
        line_p = []
        sample_m = []
        sample_q = []
        for m in range(3240):
            for n in range(6480):
                if (135 < san_she_angle[m][n] < 145 and
                        cloud_result[m][n] != 255):
                    temp1.append(PR[m][n])
                    line_m.append(m)
                    sample_m.append(n)

        for p in range(3240):
            for q in range(6480):
                if (122.508 < san_she_angle[p][q] < 135 and
                        cloud_result[p][q] != 255):
                    temp2.append(PR[p][q])
                    line_p.append(p)
                    sample_q.append(q)

        line_m = np.array(line_m)
        line_p = np.array(line_p)
        sample_m = np.array(sample_m)
        sample_q = np.array(sample_q)
        dim = line_m.shape[0]

        mn_matrix = np.zeros((2, dim))
        pq_matrix = np.zeros((2, dim))
        mn_matrix[0, :] = line_m
        mn_matrix[1, :] = sample_m
        pq_matrix[0, :] = line_p
        pq_matrix[1, :] = sample_q

        numpy1 = np.array(temp1)
        numpy2 = np.array(temp2)

        if numpy1.shape[0] == numpy2.shape[0]:
            ratio_1 = numpy1 / numpy2
        else:
            print("尺寸不对应！")
        print('dim:',dim)
        print('mn_matirx:',mn_matrix.shape)
        print('pq_matrix:',pq_matrix.shape)
        print("test over!")
        return ratio_1,mn_matrix,pq_matrix
    #写文件，以写成tif为例
    def write_img(self,filename,im_proj,im_geotrans,im_data,san_she_angle,PR):
        #gdal数据类型包括
        #gdal.GDT_Byte,
        #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        #gdal.GDT_Float32, gdal.GDT_Float64

        #判断栅格数据的数据类型Q
        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
        else:
            im_bands, (im_height, im_width) = 1,im_data.shape

        #创建文件
        driver = gdal.GetDriverByName("GTiff")            #数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)              #写入仿射变换参数
        dataset.SetProjection(im_proj)                    #写入投影
        #改变数据
        # 散射角范围a:135-143.23 b:125-135 点得个数为3935
        # ********a:135-145    b:122.508-135 点的个数为：4915 # 实验采用第二个散射角范围(点的个数多，较为直观的显示)
        temp1 = []
        temp2 = []
        line_m = []
        line_p = []
        sample_n = []
        sample_q = []
        for m in range(3240):
            for n in range(6480):
                if (135 < san_she_angle[m][n] < 145 and
                        im_data[m][n] != 255):
                    temp1.append(PR[m][n])
                    line_m.append(m)
                    sample_n.append(n)

        for p in range(3240):
            for q in range(6480):
                if (122.508 < san_she_angle[p][q] < 135 and
                        im_data[p][q] != 255):
                    temp2.append(PR[p][q])
                    line_p.append(p)
                    sample_q.append(q)

        line_m = np.array(line_m)
        line_p = np.array(line_p)
        sample_n = np.array(sample_n)
        sample_q = np.array(sample_q)
        dim = line_m.shape[0]
        # 矩阵的列没有对齐，一般情况用这个方法就行了：arrayname = np.concatenate(arrayname, axis=0) 。  arrayname是自己定义的矩阵名字
        mn_matrix = np.zeros((2, dim))
        pq_matrix = np.zeros((2, dim))
        mn_matrix[0, :] = line_m
        mn_matrix[1, :] = sample_n
        # pq_matrix[0, :] = np.concatenate(line_p,axis=0)
        pq_matrix[1, :] = sample_q
        numpy1 = np.array(temp1)
        numpy2 = np.array(temp2)
        if numpy1.shape[0] == numpy2.shape[0]:
            ratio_1 = numpy1 / numpy2
        else:
            print("尺寸不对应！")
        print("test over!")





        for i in range(dim):
           # mn = mn_matrix[:, i]#该矩阵中的数据为[a,b]
           # pq = pq_matrix[:, i]#该矩阵中的数据为[c,d]
           m=int(mn_matrix[0,i])
           n=int(mn_matrix[1,i])
           p=int(pq_matrix[0,i])
           q=int(pq_matrix[1,i])
           if ratio_1[i]>1 :#and  im_data[m,n]==100 and  im_data[p,q]==100
                #print('mn_matrix[0,i]:',mn_matrix[0,i])
                #print('mn_matrix[1,i]:',mn_matrix[1,i])
               im_data[m,n]=40#标记为冰云
               im_data[p,q]=40
               #im_data[pq_matrix[0,i],pq_matrix[1,i]]=200
           if ratio_1[i] < 1 :#and im_data[m, n] == 100 and im_data[p, q] == 100
               im_data[m,n]= 200#标记为水云
               im_data[p,q]=200
                #im_data[pq_matrix[0,i], pq_matrix[1,i]]=20
           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__":
    path=os.path.abspath('')
    run = GRID()
    proj,geotrans,data = run.read_img('data\\cloud_indicator.tif')#读数据

    relative_azimuth_angle_1=run.angle('data\\phi.tif')
    solar_zenith_angle_1=run.angle('data\\thetas.tif')
    view_zenith_angle_1=run.angle('data\\thetav.tif')

    cos_solar_zenith_angle_1 = run.cos_angle(solar_zenith_angle_1)
    cos_view_zenith_angle_1 = run.cos_angle(view_zenith_angle_1)
    sin_solar_zenith_angle_1 = run.sin_angle(solar_zenith_angle_1)
    sin_view_zenith_angle_1 = run.sin_angle(view_zenith_angle_1)
    cos_relative_azimuth_angle_1 = run.cos_r_a_angle(relative_azimuth_angle_1)

    san_she_1 = run.san_angle_1(cos_solar_zenith_angle_1, cos_view_zenith_angle_1)
    san_she_2 = run.san_angle_2(sin_solar_zenith_angle_1, sin_view_zenith_angle_1)
    san_she_3 = run.san_angle_3(san_she_2, cos_relative_azimuth_angle_1)
    san_she_4 = san_she_1 + san_she_3
    result = run.acos_san_angle(san_she_4)
    result_1=np.array(result)
    sanshe_angle=result_1*180

    I865P = gdal.Open('data\\865.tif')  # fill_value=32767,-32767
    PR=I865P.ReadAsArray(0, 0, im_width, im_height)
    PR=np.array(PR)
    #提取865nm波段第一个角度的数据
    PR0=PR[0]

    #ratio,mn_matrix,pq_matrix=run.cunzai(sanshe_angle,data,PR0)
    run.write_img('data\cloud_indicator_Rewrite.tif',proj,geotrans,data,sanshe_angle,PR0)  # 写数据
    print('数据处理结束!')

实现结果：

黄色：水云
紫色：冰云
蓝色：混合云
绿色：未定