空间资料处理综合实习

深空探测实习

平面椭圆

1.1 原理

思路：只考虑平面椭圆轮廓的反射，各向同性单次反射，假设光照方向和观测视角固定。对于固定的光照方向，与光照方向平行的椭圆的两条切线所框起来的受到光照的那部分就是光所能照到椭圆的全部；同理，对固定的观测视角，所能观测到的椭圆的的全部也是如此。

在某些情况，两个部分会产生交集，这个部分既能被光照到，又可能被观测到，即有效部分。有效部分相对于光照方向的垂直高度为d，光照强度为n（j/m），光反射之后能被观测到的概率为p，则光照强度L=n * d * p。p的变化很小，n是一个常数，令n*p=1，则L=d。

1.2 代码实现

1.2.1 思路

利用matplotlib.animation的FuncAnimation函数实现动态效果。函数的原理：

from matplotlib.animation import FuncAnimation

ani = FuncAnimation(fig, update, frames, interval)

Fig是画布对象，frames是一个数组，interval是控制时间间隔的参数，单位是ms，updata是用于更新画布的函数，它只有一个参数，每隔一定时间接收frames的数作为输入。

先创建椭圆生成函数，这里使用椭圆的参数方程，为了可以控制旋转，还要再加上旋转矩阵：

x'(t) = a * cos(t) * cos(θ) - b * sin(t) * sin(θ)

y'(t) = a * cos(t) * sin(θ) + b * sin(t) * cos(θ)

这里，θ 是旋转角度，a 和 b 分别是椭圆的长轴半径和短轴半径，t 是参数，范围为 [0, 2π]。

对于垂直高度d的计算，可通过计算切点的坐标来获得；这里我设置光照方向为水平向右，观测方向垂直向上，对应的切点就是椭圆轮廓的极点，在这种特殊情况，d可以通过左极点的y坐标减去下极点的y坐标获得。

1.2.2 代码

'''
目标：只考虑平面椭圆轮廓，只考虑反射，各向同性且反射光线均分，假设光照方向和观测视角固定，模拟椭圆的光变曲线
原理：假设光照方向水平向右，观测视角垂直向上，对于观测者来说，L=n*d*alpha/180，其中n为每米的光照强度，alpha为观测者能接收到的反射范围（忽略微小偏差）
，则光度变化只与左极点和下极点的垂直高度差有关，假设令n*alpha/180=1，则光度L=d。

将椭圆的参数方程代入，得到旋转后的椭圆参数方程：
x'(t) = a * cos(t) * cos(θ) - b * sin(t) * sin(θ)
y'(t) = a * cos(t) * sin(θ) + b * sin(t) * cos(θ)
这里，θ 是旋转角度，a 和 b 分别是椭圆的长轴半径和短轴半径，t 是参数，范围为 [0, 2π]。

运行代码遇到警告：dline.set_xdata(x_axis)
MatplotlibDeprecationWarning: Setting data with a non sequence type is deprecated since 3.7 and will be remove two minor releases later
不推荐使用非序列数据类型，不过并不影响代码的正确性,可以把代码改为dline.set_xdata([x_axis])即可
'''

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


#生产椭圆的函数,返回椭圆line2D对象，和极点位置信息，以及有效值d和位置信息
def ellipse(xc, yc, a, b, theta, angle):  #a,b分别是短、长半轴
    angle = np.deg2rad(angle) #转换弧度
    x_unrotated=xc + a * np.cos(theta)
    y_unrotated=yc + b * np.sin(theta)
    # 应用旋转矩阵
    x = x_unrotated * np.cos(angle) - y_unrotated * np.sin(angle)
    y = x_unrotated * np.sin(angle) + y_unrotated * np.cos(angle)
    line1.set_data(x,y)
    #获取极点信息
    x_min,x_max=min(x),max(x)
    y_min,y_max=min(y),max(y)
    #算d值
    indices = np.where(x == x_min)
    d=y_max-min(y[indices]) #这里min（）是防止对应y有两个
    #计算d的位置信息，用于显示
    index=np.where(y==y_min)
    x_axis = min(x[index])
    temp=(7-min(y[indices]))/14
    return line1,[x_min,x_max,y_min,y_max],d,temp,(7-y_max)/14,x_axis

#用于更新图像
def update(angle):
    #设置基础参数
    theta = np.linspace(0, 2 * np.pi, 360)
    xc, yc = 0, 0
    a, b = 3, 5
    # 更新数据
    line1,points,d,ymin,ymax,x_axis=ellipse(xc, yc, a, b, theta, angle)
    #ymin,ymax,x_axis是为了显示d值线设置的
    dline.set_xdata([x_axis])
    dline.set_ydata([ymin,ymax])
    light.append(d)
    t.append(len(light))
    line2.set_data((t,light))
    for i in range(2):
        x_lines[i].set_xdata([points[i]])
        y_lines[i].set_ydata([points[i+2]])


if __name__=='__main__':
    # 解决坐标轴中不能出现负号的问题
    plt.rcParams['axes.unicode_minus'] = False
    #解决中文显示问题
    plt.rcParams['font.sans-serif'] = ['SimHei']
    # 创建fig，axes对象
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6), facecolor='gray')
    # 创建line2D和垂直线对象
    line1, = ax1.plot([], [], 'r', lw=2)  # 用于接收椭圆坐标
    line2, = ax2.plot([], [], 'b', lw=2,label='brightness')  # 用于接收光强
    # 表示视野 垂直从下往上看
    xline_min = ax1.axvline(0, color='black', linestyle='-', linewidth=1, label='观测视野')
    xline_max = ax1.axvline(0, color='black', linestyle='-', linewidth=1)
    # 表示辐射范围
    yline_min = ax1.axhline(0, color='yellow', linestyle='-', linewidth=1, label='光照视野')
    yline_max = ax1.axhline(0, color='yellow', linestyle='-', linewidth=1)
    #用于显示d值
    dline = ax1.axvline(0, color='g', linestyle='-', linewidth=3,label='d值')
    x_lines = [xline_min, xline_max]
    y_lines = [yline_min, yline_max]
    # 用于接收数据
    light = []
    t = []
    #显示动画，100帧，360度旋转
    ani = FuncAnimation(fig, update, frames=np.linspace(0, 720), interval=110)

    ax1.set(xlim=[-7,7],ylim=[-7,7],title='模拟图',yticks=[],xticks=[])
    ax2.set(xlim=[0, 70], ylim=[0, 10], title='result', yticks=[], xticks=[],xlabel='time',ylabel='brightness')
    # # 设置坐标轴范围，便于显示
    # ax1.set_xlim(-7, 7)
    # ax1.set_ylim(-7, 7)
    # ax2.set_xlim(0, 70)
    # ax2.set_ylim(0, 10)
    # #设置label，便与显示
    # ax1.set_title('模拟图')
    # ax1.set_xticks([])
    # ax1.set_yticks([])
    # ax2.set_title('result')
    # ax2.set_xlabel('time')
    # ax2.set_ylabel('brightness')
    ax1.legend(loc=(2))
    ax2.legend(loc=('upper left'))

    # ani.save('result.gif','pillow',fps=100) #保存文件
    plt.show()

1.3 结果展示 

 

 

 空间资料处理综合实习

前言：本次实验的目的是使用ENVI和代码程序对影像进行分类，对结果进行比较（软件：ENVI5.31；程序语言：python3.11，使用sklearn库进行操作）。主要做了以下操作：选好了成都市成华区作为研究区，在地理空间数据云上下载好数据，完成预处理，包括辐射定标、大气校正、图像融合和感兴趣区裁剪；在研究区选取了2组建筑/道路、白色建筑、裸地、水体、植被5种类别的roi分别作为训练样本和检验样本；用ENVI和sklearn（python包）分别进行支持向量机分类；之后用5折交叉验证法来优化分类器参数，再进行分类；最后进行分类后处理，用多种分类方法的分类结果中分的最好的类别进行集成。分别与ENVI分类结果进行对比和在arcgis上做差异性分析，最后制图。 

 

第一章 数据预处理

1.1数据源

数据来自地理空间数据云，2014年8月13日LandSat8影像，云量少于10%。

1.2数据预处理

1.2.1 几何处理

由于此类产品分辨率较低且没有明显的扭曲变形，几何处理对于后续操作的作用不大，所以略去此操作。

1.2.2像素处理

1.2.2.1 辐射定标

ENVI5.31：Radiation correction/Radiation calibration

将仪器所测量的辐射值转换为表征地物表面和大气的物理参数，如地表反射率、辐射亮温度、大气透过率等，以便于后续的遥感数据分析和应用。

1.2.2.2 大气校正

ENVI5.31：Radiation correction/Atmospheric correction

去除大气的干扰，让数据反应的是真实的地表特征。

1.2.2.3 图像融合 

ENVI5.31：:Image Sharpening/Gram-Schmidt Pan Sharp

让多光谱和全色融合，输出结果会既有多光谱的波段属性又有全色的高分辨率

1.2.2.4 裁剪  

ENVI5.31：Regions of Interest/Subset Data from ROIs

根据行政区边界裁剪感兴趣区域。

1.3预处理结果显示

 

 最终结果