该篇是在申请南航博士时，布置的一个探索题。由于第一次接触遥感图像，我只是机械的完成了每题，没有注重题与题之间的关系，尚未达到老师预期。本来想丢掉，但是还是上传一下吧。

博士探索题

题1 通过网络等手段找到任意某颗遥感卫星的一幅图像数据和辅助数据，完成以下工作：

（辅助数据是指在主体数据处理中具有辅助作用的数据。辅助数据不构成数据处理的主要实体，但是如果在数据处理过程中加入辅助数据，将有助于提高数据处理精度。在遥感数据输入地理信息系统的过程中，需要对遥感数据进行有效地分类。在分类过程中借助地理信息系统中的辅助数据（如该地区高程数据、地面坡度数据、地表坡向数据及土地利用数据等）的支持，可改善遥感数据的分类精度，以达到地理信息系统中数据处理的精度要求）

1.对图像数据和辅助数据进行解译，解释说明图像的各个参数含义（编程完成）

（解译：通过遥感图像所提供的各种识别目标的特征信息进行分析、推理和判断，最终达到识别目标或现象的目的。遥感图像解译可以说是遥感成像过程的逆过程，即从遥感对地面实况的模拟影响中提取遥感信息、反演地面原型的过程）

1).数据来源

地理空间数据云 (gscloud.cn)

高级检索>>Landsat 8 OLI_TIRS 卫星数字产品>>设置搜索范围

2).数据解读（下述内容部分参考【Landsat 8】遥感影像文件内容及命名规则_landsat8影像命名规则_Uyoin的博客-CSDN博客

a. 文件构成

- 压缩包名称：LC81180382020357LGN00.tar

压缩包命名分段	含义
LC8	Landsat-8 卫星
118	WRS Path
038	WRS Row
2020	影像获取年份
357	影像获取日期，1月1日起累计
LGN	接站代码
00	产品级别
.tar	压缩格式后缀

（WRS（Worldwide Reference System）：是Landsat卫星采用的全球参考系统,用以区分全球各区域对应的Landsat系列卫星影像编号,其用“Path”与“Row”两个数值确定影像的编号与位置）

- 文件包含：2个txt文档与11个tif文件(命名方LXSS_LLLL_PPPRRR_YYYYMMDD_yyyymmdd_CC_TX_BB）

子文件命名分段	含义
L	Landsat 卫星
X	传感器（’C’= OLI与TIRS组合；’O’只有OLI；’T’只有TIRS;’ E’=ETM+；’T’=TM;’M’=MSS）
SS	卫星编号，’08’指Landsat-8 卫星
LLL	处理等级
PPP	WRS Path
RRR	WRS Row
YYYYMMDD	影像获取年月日
yyyymmdd	影像处理年月日
CC	编码
TX	‘RT’实时层;’T1’第一层;’T2’第二次
BB	ANG:投影信息;MTL:影像元数据；B1~B11:不同波段遥感影像；BQA：影像质量信息

b. 文件解释

-TXT文件（MTL文件）包含以下信息：

文件元数据

命名	含义
来源 (`ORIGIN`)	数据来源。
请求ID (`REQUEST_ID`)	数据请求的唯一标识符。
Landsat场景ID (`LANDSAT_SCENE_ID`)	图像的唯一标识符。
Landsat产品ID (`LANDSAT_PRODUCT_ID`)	产品的独特标识。
收集编号 (`COLLECTION_NUMBER`)	数据收集批次的编号
文件日期 (`FILE_DATE`)	数据文件的生成日期
站点ID (`STATION_ID`)	接收遥感数据的地面站标识
处理软件版本 (`PROCESSING_SOFTWARE_VERSION`)	处理数据的软件版本

产品元数据

命名	含义/值
数据类型 (`DATA_TYPE`)	数据的类型(数据处理的级别）：本文件为Level 1 Terrain Corrected Product（L1TP)，即第一级地形校正产品。这意味着图像已经经过了地形校正，以校正地表地形的影响。这种类型的数据对于大多数地表特征分析和应用来说是最有用的，因为它提供了准确的地理位置信息。
数据类别 (`COLLECTION_CATEGORY`）	数据质量等级：本文件为T1，最高质量的数据类别，适用于科学和应用目的。
高程数据源 (`ELEVATION_SOURCE`)	地形校正过程中的高程数据源：本文为GLS2000，全球地表覆盖数据集，提供地表的高程信息，提高数据的几何精度和可靠性。
输出格式 (`OUTPUT_FORMAT`)	文件输出格式：”GEOTIFF”。
航天器ID(`SPACECRAFT_ID`)	航天器编号：”LANDSAT_8”
传感器ID (`SENSOR_ID`)	传感器型号：”OLI_TIRS”。(OLI是一种高级成像仪器，设计用于捕获地表图像；TIRS是用于捕捉地表温度或热特性的热红外传感器)
WRS路径和行 (`WRS_PATH`, `WRS_ROW`)	用于确定影像编号与位置，遥感卫星轨道上的路径为path（东西）,遥感卫星在路径上的行为row（南北）：（118，38）
观测角度（`NADIR_OFFNADIR`）	本文件为NADIR。(NADIR 观测角度为传感器以垂直地球表面方向观测地面，有助于测量地面特征。OFFNADI为有一定倾斜角度，提供更广阔覆盖范围)
目标观测区域(`TARGET_WRS_PATH，TARGET_WRS_ROW`	（118，38）
获取日期 (`DATE_ACQUIRED`)	2020-12-22
场景中心时间 (`SCENE_CENTER_TIME`)	“02:25:20.1010799Z”
角点坐标 (`CORNER__LAT_PRODUCT`, `CORNER__LON_PRODUCT`)	（32.78382，120.68188），（32.80518，123.15168），（30.66079，120.73416），（30.68047，123.14825）
在投影坐标系中的坐标(`CORNER__PROJECTION_X_PRODUCT，CORNER__PROJECTION_Y_PRODUCT`）	（282900，3629700），（514200，3629700），（282900，3394200），（514200，3394200）
全色波段图像的高度与宽度(`PANCHROMATIC_LINES，PANCHROMATIC_SAMPLES`)	（15701，15421）
反射波段图像高度与宽度(`REFLECTIVE_LINES，REFLECTIVE_SAMPLES`)	（7851，7711）
热红外波段图像高度与宽度(`THERMAL_LINES，THERMAL_SAMPLES`)	（7851，7711）
文件名称`FILE_NAME_BAND_1~11 ,FILE_NAME_BAND_QUALITY, ANGLE_COEFFICIENT_FILE_NAME, METADATA_FILE_NAME)`	命名方式见2）a
标识文件元数据集信息(`CPF_NAME，BPF_NAME_OLI,BPF_NAME_TIRS,RLUT_FILE_NAME`)	校准参数名，生成产品的偏倚参数文件(BPF)的文件名（传感器：OLI与TIRS)，生成产品的响应线性化查找表。

图像元数据

命名	含义/值
云覆盖度(`CLOUD_COVER`)	图像中被云层覆盖的比例：0.61%
陆地上的云覆盖度(`CLOUD_COVER_LAND`)	陆地地表被云层覆盖的比例:1.26%
OLI传感器图像质量评分(`IMAGE_QUALITY_OLI`)	OLI波段的图像质量评估:9(受云覆盖、大气干扰、传感器问题等影响等级通常由1到10，1最低)
TIRS传感器图像质量评分(`IMAGE_QUALITY_TIRS`)	TIRS波段的图像质量评估:9
TIRS传感器扫描镜模型(`IMAGE_QUALITY_TIRS`)	TIRS扫描镜的模型:FINAL(通过改进、校准确定的最终版本)
TIRS传感器扫描镜位置状态(`TIRS_SSM_POSITION_STATUS`)	ESTIMATED(表示估算得到)
TIRS传感器漂移光矫正源(`TIRS_STRAY_LIGHT_CORRECTION_SOURCE`)	TIRS漂移光矫正的数据源:TIRS
卫星的滚转角度(`ROLL_ANGLE`)	卫星相对于参考轨道的滚转角度:-0.001
太阳方位角(`SUN_AZIMUTH`)	太阳在地平面上的方向：156.98570422
太阳高度角(`SUN_ELEVATION`)	太阳位于地平面上的高度：31.14958435
地球到太阳的距离(`EARTH_SUN_DISTANCE`)	地球与太阳之间的平均距离：0.9836681
各个波段的饱和度(`ATURATION_BAND_1~SATURATION_BAND_9`)	对应波段是否饱和。”Y” 表示饱和，”N” 表示未饱和。
地面控制点版本(`GROUND_CONTROL_POINTS_VERSION`)	地面控制点的数据版本：4
地面控制点模型(`GROUND_CONTROL_POINTS_MODEL`)	地面控制点的模型：102
几何校正均方根误差(`GEOMETRIC_RMSE_MODEL`)	几何校正模型的均方根误差：8.773
几何校正均方根误差的垂直和水平分量(`GEOMETRIC_RMSE_MODEL_Y,GEOMETRIC_RMSE_MODEL_X`)	几何校正均方根误差的垂直和水平分量：5.847/6.541
验证地面控制点的数量(`GROUND_CONTROL_POINTS_VERIFY`)	验证地面控制点的数量：39
验证几何校正均方根误差(`GEOMETRIC_RMSE_VERIFY`)	验证几何校正均方根误差：9.995
OLI波段的截断方式(`TRUNCATION_OLI`)	UPPER（表示向上截断）

其他

命名	含义/值
最大-最小辐射值(`MIN_MAX_RADIANCE`)	1~11各个波段的辐射最大/最小值
最大-最小反射率值(`MIN_MAX_REFLECTANCE`)	1~9各个波段的反射率最大/最小值
最大-最小像素值(`MIN_MAX_PIXEL_VALUE`)	1-11各个波段的量化最大/最小值
辐射/反射系数(`RADIOMETRIC_RESCALING`)在后续的辐射定标和反射定标中需要使用	1-11各个波段的辐射系数最大/最小值（MULT:乘法因子；ADD：加法因子），1-9各个波段的反射率系数最大/最小值（MULT:乘法因子；ADD：加法因子）（乘法因子和加法因子是用于将数字值转换为物理量（如辐射亮度或反射率）的参数，以确保遥感图像与地表特性的物理量之间的关系准确和可靠）
热红外定标参数(`TIRS_THERMAL_CONSTANTS`)在后续的温度定标中需要使用	TIRS波段10/11的热红外定标参数
地图投影类型(`MAP_PROJECTION`)	图像使用的地图投影类型：UTM（UTM：横向墨卡托投影。一种常用的地图投影方法，通常用于将地球上的表面划分为小块，每块都使用平面投影来表示）
大地基准(`DATUM`)	图像使用的大地基准：WGS84（坐标系统使用经度（Longitude，东西）和纬度（Latitude，南北）来表示地球上的位置，这个坐标系统的原点位于地球的质心，并且与地球的椭球体模型相对应）
椭球体(`ELLIPSOID`)	图像使用的椭球体：WGS84
UTM投影的区域代号(`UTM_ZONE`)	图像使用的UTM投影的区域：51（UTM划分的小块区域代号51）
不同波段的栅格单元大小(`GRID_CELL_SIZE_PANCHROMATIC，GRID_CELL_SIZE_REFLECTIVE，GRID_CELL_SIZE_THERMAL`)	不同波段的栅格单元大小：（15，30，30）
图像的方向(`ORIENTATION`)	图像的方向是否朝北：NORTH_UP
重采样选项(`RESAMPLING_OPTION`)	图像重采样时所采用的方法：CUBIC_CONVOLUTION。（三次样条插值是一种数学插值方法，它使用三次多项式来逼近像素之间的值。这种插值方法在重采样过程中考虑了像素之间的连续性和平滑性，因此通常能够产生较平滑的图像）

-TXT文件（ANG文件）包含以下信息：

头文件

命名	含义/值
Landsat卫星场景的唯一标识符(`LANDSAT_SCENE_ID`)	包含传感器、日期和路径/行等信息：LC81180382020357LGN00
拍摄图像的卫星的标识符(`SPACECRAFT_ID`)	图像是由Landsat 8 卫星拍摄的LANDSAT_8
图像中的波段数量(`NUMBER_OF_BANDS`)	图像中包含的不同光谱波段的数量：11
图像中包含的波段列表(`BAND_LIST`)	列表形式展示波段

投影信息

命名	含义/值
地球椭球体的半长轴和半短轴的长度(`ELLIPSOID_AXES`)	(6378137.000000, 6356752.314200)
地图投影的类型(`MAP_PROJECTION`)	同MTL文件
投影的单位(`PROJECTION_UNITS`)	Meters(米)
地图投影所基于的大地基准(`DATUM`)	同MTL文件
UTM投影的区域代号(`UTM_ZONE`)	同MTL文件
地图投影的参数`PROJECTION_PARAMETERS`)	包括平移、旋转和缩放等。本文件全为0，表示没有变换。
四个角的地理坐标（经度和纬度）(`UL_CORNER，UR_CORNER，LL_CORNER，UR_CORNER`)	同MTL文件

星历数据（用于确定卫星的位置和运动轨迹）

命名	含义/值
星历数据的参考时刻的年份（`EPHEMERIS_EPOCH_YEAR`)	2020
星历数据的参考时刻的轮内第几天（`EPHEMERIS_EPOCH_DAY`)	357
星历数据的参考时刻的秒数（`EPHEMERIS_EPOCH_SECONDS`)	8693.716066
星历数据中包含的数据点数量（`NUMBER_OF_POINTS`)	55
一系列时间值，表示星历数据中的时间点（`EPHEMERIS_TIME`)	——（以秒为单位，从参考时刻开始逐步增加）
卫星的地心地固坐标系（ECEF坐标系）中的X、Y和Z坐标分量有关数据（`EPHEMERIS_ECEF_X、EPHEMERIS_ECEF_Y、EPHEMERIS_ECEF_Z`)	——（每个坐标对应于星历数据中的时间点，描述卫星在地球上的位置随时间的变化-）

（ECEF坐标系：以地球质心为原点，地球的自转轴为坐标系的Z轴，并与地球的自转一起固定。X轴指向本初子午线（通常是格林威治子午线）与赤道的交点，指向东方。y轴与x轴垂直，指向南方）

太阳矢量（用于确定太阳在地球上不同时间和位置的方向和强度，其中的值基本与星历数据相同）

命名	含义/值
太阳矢量数据的参考时刻的年份（`SOLAR_EPOCH_YEAR`)	2020
太阳矢量数据的参考时刻的轮内第几天（`SOLAR_EPOCH_DAY`)	357
太阳矢量数据的参考时刻的秒数（`SOLAR_EPOCH_SECONDS`)	8693.716066
地球与太阳之间的平均距离（`EARTH_SUN_DISTANCE`)	0.98366810天文单位（天文单位是地球与太阳之间的平均距离，约等于1.496 × 10^8公里）
太阳矢量数据中包含的数据点数量（`NUMBER_OF_POINTS`)	55
一系列时间值，表示太阳矢量数据的时间点（`SOLAR_TIME`)	——（以秒为单位，从参考时刻开始逐步增加）
太阳的地心地固坐标系（ECEF坐标系）中的X、Y和Z坐标分量有关数据（`SOLAR_ECEF_X、SOLAR_ECEF_Y、SOLAR_ECEF_Z`)	——（每个坐标对应于太阳矢量数据中的时间点，描述太阳在地球上的位置随时间的变化）

各个波段的信息，以波段1为例

命名	含义/值
传感器芯片组件数量（`BAND01_NUMBER_OF_SCAS`)	该频带内有14个独立的传感器芯片组件
L1T级别数据的行宽（`BAND01_NUM_L1T_LINES，BAND01_NUM_L1T_SAMPS`)	（7851，7711）
L1T级别图像的四个角点位置（`BAND01_L1T_IMAGE_CORNER_LINES，BAND01_L1T_IMAGE_CORNER_SAMPS`)	( 4.182939, 1476.124065, 7848.396736, 6352.972113) ( 1523.595155, 7705.349531, 6176.008193, 0.230545)
像素尺寸（`BAND01_PIXEL_SIZE`)	30米/像素
数据开始时间和每行的时间间隔（`BAND01_START_TIME，BAND01_LINE_TIME`)	10.495779/0.004236000
平均高度（`BAND01_MEAN_HEIGHT`)	250米
L1R和L1T级别的平均行和列（`BAND01_MEAN_L1R_LINE_SAMP，BAND01_MEAN_L1T_LINE_SAMP`)	（3763.889, 3468.214）（3930.371, 3860.961）（L1R级别：L1R代表辐射校正级别，通常是遥感数据的第一处理级别；L1T级别：L1T代表地表校正级别，通常是遥感数据的第二处理级别）
卫星平均位置矢量（`BAND01_MEAN_SAT_VECTOR`)	(-0.002261318, -0.003964814, 0.996297368)
卫星各轴位置系数（`BAND01_SAT__NUM_COEF，BAND01_SAT__DEN_COEF`)	——————
太阳平均位置矢量（`BAND01_MEAN_SUN_VECTOR`)	( 0.334537344, -0.787571723, 0.517343518)
太阳各轴位置系数（`BAND01_SUN__NUM_COEF，BAND01_SUN__DEN_COEF`)	——————
通道中的SCA列表（`BAND01_SCA_LIST`)	一共14个。[1,2….14]。
各个传感器平均高度（`BAND01_*_MEAN_HEIGHT`)	——————
各个传感器的L1R和L1T级别的平均行和列（`BAND01__MEAN_L1R_LINE_SAMP，BAND01__MEAN_L1T_LINE_SAMP`)	——————

-TIF文件包含以下信息：

TIF文件包含多个波段的信息。目的在于，不同地物对于不同波段的辐射具有不同的响应，通过捕捉多个波段的信息，可以更全面地理解地表的各种特征和属性。

波段名（Band Name）	波长范围（Bandwidth）(µm)	空间分辨率（Resolution） (m)
Band 1 Coastal （海岸波段，可见）	0.43 – 0.45	30
Band 2 Blue （蓝波段，可见）	0.45 – 0.51	30
Band 3 Green （绿波段，可见）	0.53 – 0.59	30
Band 4 Red （红波段）	0.64 – 0.67	30
Band 5 NIR （近红外波段）	0.85 – 0.88	30
Band 6 SWIR 1 （短波红外1）	1.57 – 1.65	30
Band 7 SWIR 2 （短波红外2）	2.11 – 2.29	30
Band 8 Pan （全色波段）	0.50 – 0.68	15
Band 9 Cirrus （卷云波段）	1.36 – 1.38	30
Band 10 TIRS 1 (热红外1)	10.6 – 11.19	100
Band 11 TIRS 2 (热红外2)	11.5 – 12.5	100

（波段1-9为陆地成像仪：Operational Land Imager (OLI)拍摄；波段10-11为热红外传感器Thermal Infrared (TIRS)拍摄）

3).解译过程（下述部分内容参考Python 处理遥感影像入门 - 知乎 (zhihu.com)）

a. 导入必须库

import rasterio  #主要用于空间栅格/矢量数据处理
import rasterio.plot
import pyproj    #用于坐标转换
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

b. 打印文件信息

1
2
3

image_path ='D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B4.TIF'
with rasterio.open(filepath) as src:
    print(src.profile)

1
2

{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 7711, 'height': 7851, 'count': 1, 'crs': CRS.from_epsg(32651), 'transform': Affine(30.0, 0.0, 282885.0,
       0.0, -30.0, 3629715.0), 'blockysize': 1, 'tiled': False, 'interleave': 'band'}

src.profile 包含了文件的元数据。我们得知文件格式为GeoTIFF，像元值为无符号整型，‘nodata’像元值为被定义，图片分辨率为7711乘以7851，数据仅包含一个波段，UTM坐标系，并且使用了仿射变换。数据存取以1*1的数据块为单元，存储在AWS中。

c. 显示图像

方案1：

with rasterio.open(image_path) as src:
    image = src.read(1)  
# 图片显示
#plt.figure(figsize=(5, 5),dpi=300)
plt.imshow(image)
plt.title('Overview - Band 4 {}'.format(image.shape))
plt.colorbar(label='Digital Number (DN) Values')
plt.show()

受地球自转以及卫星轨道影响，上述图片中存在着无效数据0。图像数据原格式为uint16,无效区域值为NaN。为去掉四周无效数据，需要先将unit16转为float32，最后显示图像。代码及结果展示如下。

with rasterio.open(image_path) as src:
    image = src.read(1)  
    image = image.astype('f4') #图片转数据格式，Nan转为0
    image[image==0] = np.nan #去除边界为0

# 图片显示
plt.imshow(image)
plt.title('Overview - Band 4 {}'.format(image.shape))
plt.colorbar(label='Digital Number (DN) Values')
plt.show()

方案2：

"""读取图像波段的子集"""
def read_band_subset(image_path, subset_size=(1000, 1000)):
    with rasterio.open(image_path) as src:
        width, height = src.width, src.height #获取宽与高
        start_x, start_y = width // 2 - subset_size[0] // 2, height // 2 - subset_size[1] // 2 #取中心位置
        window = rasterio.windows.Window(start_x, start_y, subset_size[0], subset_size[1])
        band_subset = src.read(1, window=window)
    return band_subset

b4_band_subset = read_band_subset(image_path)


# 展示图像和部分辅助数据的内容
plt.imshow(b4_band_subset)
plt.title('Overview - Band 4 {}'.format(b4_band_subset.shape))
plt.colorbar(label='Digital Number (DN) Values')
plt.show()

上图为选取了整张图片的中心区域并展示。

d. 波段组合

组合4、3、2波段（自然色）

#数据导入，使用上述读取中心区域的函数
b4_band_subset = read_band_subset('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B4.TIF')
b3_band_subset = read_band_subset('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B3.TIF')
b2_band_subset = read_band_subset('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B2.TIF')


# 标准化波段数据
b4_band_subset= b4_band_subset / b4_band_subset.max()
b3_band_subset = b3_band_subset / b3_band_subset.max()
b2_band_subset = b2_band_subset/ b2_band_subset.max()

# 组合成RGB图像
rgb_image = np.dstack((b4_band_subset, b3_band_subset, b2_band_subset))

# 展示彩色图像
plt.imshow(rgb_image)
plt.title('RGB Composite Image (B4, B3, B2)')
plt.show()

此外还有一些波段组合展示如下图所示：

e. 计算NDVI（归一化植被制数）

归一化植被指数（NDVI）通过测量近红外（植被强烈反射）和红光（植被吸收）之间的差异来量化植被。NDVI的范围始终为-1至+1。但是每种类型的土地覆盖并没有明确的界限。

例如，当值为负数时，很可能是水。另一方面，如果NDVI值接近+1，则很有可能是茂密的绿叶。然而，当NDVI接近于零时，没有绿叶，甚至可能是城市化区域。
$$
NDVI = (nir-red)/(nir+red)
$$
其中，nir指近红外通道，red指红色通道。与其他波长相比，健康的植被（叶绿素）反射更多的近红外（NIR）和绿光。但是它吸收更多的红色和蓝色光。

def read_band_subset(image_path, subset_size=(1000, 1000)):
    """读取图像波段的子集"""
    with rasterio.open(image_path) as src:
        # 计算要读取的子集的位置
        width, height = src.width, src.height
        start_x, start_y = width // 2 - subset_size[0] // 2, height // 2 - subset_size[1] // 2
        window = rasterio.windows.Window(start_x, start_y, subset_size[0], subset_size[1])
        band_subset = src.read(1, window=window).astype(float)  
    return band_subset

red = read_band_subset('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B4.TIF')
nir = read_band_subset('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B5.TIF')

# 计算NDVI
ndvi = (nir -red) / (nir + red)

# 展示NDVI图像
plt.imshow(ndvi, cmap='RdYlGn')
plt.title('NDVI Image')
plt.colorbar(label='NDVI Value')
plt.show()

f. 更多指数

归一化差异水体指数 (NDWI): 用于监测水体和水分含量。通常使用绿色和近红外波段来计算。
$$
NDWI = （绿光-近红外）/（绿光+近红外）
$$
归一化差异建筑区域指数 (NDBI):用于识别城市化区域或建筑物。它通常利用近红外和短波红外波段的数据。
$$
NDBI = （短板红外波-近红外）/（短板红外波+近红外）
$$
土壤调节植被指数 (SAVI): 类似于NDVI，但包含一个用于减少土壤亮度影响的调节因子。这对于低植被覆盖区域特别有用。L一般取0.5。
$$
SAVI = （1+L）*（近红外-红光）/（近红外+红光+L）
$$

g. 保存结果(NDVI)到本地

red_path =('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B4.TIF')
nir_path =('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B5.TIF')
with rasterio.open(red_path) as src:
    red = src.read(1).astype('f4')  
with rasterio.open(nir_path) as src:
    nir = src.read(1).astype('f4') 
# 计算NDVI
ndvi = (nir-red) / (nir + red)
# 展示NDVI图像
plt.imshow(ndvi, cmap='RdYlGn')
plt.title('NDVI Image')
plt.colorbar(label='NDVI Value')
plt.show()

在保存为’tif’文件时，需要修改为仿射变换（线性与平移变换的叠加）矩阵。具体操作如下：

# 计算下采样因子，本例子中下采样比例为1。
with rasterio.open(red_path) as src:
    original_width, original_height = src.width, src.height #原图尺寸
downsampled_width, downsampled_height = ndvi.shape[1], ndvi.shape[0]#NDVI尺寸
width_factor = original_width / downsampled_width
height_factor = original_height / downsampled_height
width_factor,height_factor

Out[1]:

1	(1.0, 1.0)

# 计算变换后的仿射矩阵
with rasterio.open(red_path) as src:
    transform = src.transform
transform = transform * rasterio.Affine.scale(width_factor, height_factor)
transform

Out[2]:

1 2	Affine(30.0, 0.0, 282885.0, 0.0, -30.0, 3629715.0)

# 保存
ndvi_image_path = 'D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_NDVI.TIF'#保存地址

with rasterio.open(red_path) as src: #打开源文件
    profile = src.profile
    profile.update(dtype=rasterio.float32, transform=transform, width=downsampled_width, height=downsampled_height) #更新信息

    with rasterio.open(ndvi_image_path, 'w', **profile) as dst:
        dst.write(ndvi.astype(rasterio.float32), 1) #保存

测试保存内容。

with rasterio.open(ndvi_image_path) as src:
    image = src.read(1)  
# 图片显示
#plt.figure(figsize=(5, 5),dpi=300)
plt.imshow(image)
plt.title('Overview - Ndvi {}'.format(image.shape))
plt.colorbar(label='Digital Number (DN) Values')
plt.show()

h. 空间索引及提取像素值

遥感影像有两组坐标系，数据的坐标系为UTM坐标系。一组是图片的行列坐标，另一组是空间坐标。遥感影像中的每一个像元都对应地球上的一个地点。如需获取某个地理坐标的NDVI值，需要将经纬度坐边投影为UTM。具体操作如下：

with rasterio.open(ndvi_image_path) as src:

    utm = pyproj.Proj(src.crs)
    lonlat = pyproj.Proj(init='epsg:4326') #创建一个坐标参考系统
 
    lon,lat = (119, 31) #上理经纬度
    east,north = pyproj.transform(lonlat, utm, lon, lat) #坐标转换
 
    print('Fresno NDVI\n-------')
    print(f'lon,lat=\t\t({lon:.2f},{lat:.2f})')
    print(f'easting,northing=\t({east:g},{north:g})')
     
    row, col = src.index(east, north) # spatial --> image coordinates
    print(f'row,col=\t\t({row},{col})')
 

    value = ndvi[row, col]
    print(f'ndvi=\t\t\t{value:.2f}')

output[3]

Fresno NDVI
-------
lon,lat=		(119.00,31.00)
easting,northing=	(117989,3.43648e+06)
row,col=		(6441,-5497)
ndvi=			0.08

i. 计算一段时间内NDVI值的变化

原先的数据是2022年12月22号的数据，这里又下载了2022年8月16号的数据，计算其NDVI值。

red_path1 =('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data2/LC08_L1TP_118038_20200222_20200225_01_T1_B4.TIF')
nir_path1 =('D:/YANJIUSHENG/南航博士申请探索题/Dataset/data2/LC08_L1TP_118038_20200222_20200225_01_T1_B5.TIF')
with rasterio.open(red_path1) as src:
    red1 = src.read(1).astype('f4')  
with rasterio.open(nir_path1) as src:
    nir1 = src.read(1).astype('f4') 
# 计算NDVI
ndvi1 = (nir1-red1) / (nir1 + red1)
# 展示NDVI图像
plt.imshow(ndvi, cmap='RdYlGn')
plt.title('NDVI Image')
plt.colorbar(label='NDVI Value')
plt.show()

而后两者相减，展示变化。

# 计算NDVI的变化
ndvi_change = ndvi - ndvi1

# 可视化NDVI的变化
plt.imshow(ndvi_change, cmap='RdYlGn')
plt.colorbar(label='NDVI Change')
plt.title('NDVI Change between Two Time Periods')
plt.show()

分析：差值可以看出基本在0值附近，处于0值下方。这说明绿色植被减少，但也可能受季节影响，12.22属秋冬季，绿色植被较少，08.16为夏季，绿色植被较多。

2. 利用STK仿真载荷成像时刻轨道位置与视场

(STK软件提供分析引擎用于计算数据、并可显示多种形式的二维地图，显示卫星和其它对象如运载火箭、导弹、飞机、地面车辆、目标等。 STK的核心能力是产生位置和姿态数据、获取时间、遥感器覆盖分析。)

1）联合python与STK创建应用获取接口

from comtypes.client import CreateObject
# 创建STK11桌面应用
app=CreateObject("STK11.Application")
# STK11可见
app.Visible = 1
# 获取IAgStkObjectRoot接口
root = app.Personality2
# 查看root属性
print(type(root))

2）创建场景

from comtypes.client import CreateObject
from comtypes.gen import STKObjects

# 获取IAgStkObject接口
rootIAF = root.QueryInterface(STKObjects.IAgStkObject)
# 利用Children属性进行场景创建
scenario = rootIAF.Children.New(19,'Shanghai')

3）设置场景时间

scenario2 = scenario.QueryInterface(STKObjects.IAgScenario)
root.UnitPreferences.SetCurrentUnit('DateFormat', 'UTCG')
# 设置场景时间
scenario2.StartTime = '22 Dec 2020 02:25:20.101'
scenario2.StopTime =  '22 Dec 2020 04:00:00.000'
scenario2.CurrentTime =scenario2.StartTime

4) 创建卫星，设置参数

同样可以python创建或者内部创建。内部创建共有6种方法。这里通过搜索ID创建卫星：

5）结果展示

3D图与2D图，但角点坐标以及成像信息与之前的TXT文档及TIF显示结果有略微差别（需要继续探索原因）。

3. 利用辅助数据对图像数据进行辐射定标（编程完成）

1）简介

a. 辐射定标定义

将图像的亮度灰度值转换为绝对的辐射亮度（统一值）的过程。单位通常是 W/(m²·sr·μm)。

其中，W能量单位，用在此处为辐射的功率。m²面积单位，表示辐射亮度是如何分布在一个特定的面积上的。sr三维空间中的角度单位，用于描述辐射在空间中的分布范围，1立体角是定义在球面上、以球心为顶点、截面积为半径的平方的球冠所对应的角。μm长度单位，指定波长的范围或特定波长。因此辐射亮度为描述每单位面积、每单位立体角、在特定波长范围内的辐射功率的量。即从一个方向在特定波长上到达或离开一个表面的光的强度。

b. 为什么要进行辐射定位

在第一个问题的分析中，我们使用了像元值对同景图像内部进行了比较。然而全球资源和环境变化需要提供多时域、多区域、多种传感器且相互之间具有可比性的遥感数据。因此需要将像元值转换为统一值。

c. 还可以进行哪些定标

多光谱数据（B1~B9)可以进行辐射率定标和反射率定标。公式如下：

辐射率定标：
$$
L = M_L * Q_L + B_L
$$
其中，M_L为一个波段特定的辐射定标系数，Q_L为原图像量化值，B_L为一个波段特定的辐射定标加数。

反射率定标（像素的地面反射率）ρ ：
$$
ρ = M_ρ * Q_ρ + B_ρ
$$
其中，Mρ 为反射率定标系数，Qρ原图像量化值， Bρ 为反射率定标加数。

热红外数据(B10~B11)可以进行辐射率定标和亮度温度定标。

亮度温度定标（传感器记录的辐射值转换为等效的亮度温度）：
$$
T = K_2/(ln(K_1/L+1))
$$
其中T是亮度温度，通常以开尔文（K为单位）。K1和K2是卫星特定的定标常数(K1传感器特性有关的放射常数，K2波长有关的温度转换常数)。L是辐射亮度。通过辐射定标公式计算得到。

2) 编程实现

a. 在TXT文件中找到相关参数

相关系数在MTL.txt文件中，在前文的文件解释中MTL.txt文件的其他表格内，对应于RADIOMETRIC_RESCALING和TIRS_THERMAL_CONSTANTS 。值如下：

波段	辐射乘数/加数	反射率乘数/加数	K1 常数	K1 常数
1	1.2976E-02/-64.88072	2.0000E-05/-0.100000	————	————
2	1.3288E-02/-66.43865	2.0000E-05/-0.100000	————	————
3	1.2245E-02/-61.22265	2.0000E-05/-0.100000	————	————
4	1.0325E-02/-51.62639	2.0000E-05/-0.100000	————	————
5	6.3186E-03/-31.59277	2.0000E-05/-0.100000	————	————
6	1.5714E-03/-7.85684	2.0000E-05/-0.100000	————	————
7	5.2963E-04/-2.64817	2.0000E-05/-0.100000	————	————
8	1.1685E-02/-58.42686	2.0000E-05/-0.100000	————	————
9	2.4694E-03/-12.34717	2.0000E-05/-0.100000	————	————
10	3.3420E-04/0.10000	————	774.8853	1321.0789
11	3.3420E-04/0.10000	————	480.8883	1201.1442

b. 代码

定义计算函数>>遍历图像>>调用函数>>图像显示

import rasterio
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# 辐射定标函数
def radiance_calibration(band_data, mult, add):
    return mult * band_data + add
# 反射率定标函数,反射率的加减法因子都一样
def reflectance_calibration(band_data, mult = 2.0000E-05, add = -0.100000):
    return mult * band_data + add  # Assuming no atmospheric correction
# 亮度温度定标函数
def brightness_temperature_calibration(band_data, k1, k2):
    return k2 / np.log((k1 / band_data) + 1)

#给定参数
radiance_p = {
    1: (1.2976E-02, -64.88072),
    2: (...)
    ...
}

brightness_temp_p = {
    10: (774.8853, 1321.0789),
    11: (480.8883, 1201.1442),
}

# 遍历，逐一求取标定值并展示。
for i in range(1,12):
    band_file_path = 'D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/'+'LC08_L1TP_118038_20201222_20210310_01_T1_B'+str(i)+'.TIF'
    with rasterio.open(band_file_path) as src:
        band_data = read_band_subset（band_file_path） #原图中的部分图像进行定标
        band_radiance = radiance_calibration(band_data, *radiance_p[i]) 
        plt.imshow(band_radiance, cmap='RdYlGn')
        plt.colorbar(label='辐射率标定值')
        plt.title(f'band{i}辐射标定')
        plt.show()
        
        if i<10:
            band_reflectance = reflectance_calibration(band_data)
            plt.imshow(band_reflectance, cmap='RdYlGn')
            plt.colorbar(label='反射率标定值')
            plt.title(f'band{i}反射率标定值')
            plt.show()     
        else:
            band_K =  brightness_temperature_calibration(band_data, *brightness_temp_p[i])
            plt.imshow(band_K, cmap='RdYlGn')
            plt.colorbar(label='亮度温度标定值')
            plt.title(f'band{i}亮度温度标定')
            plt.show()

3）结果展示

辐射标定

反射率标定

亮度温度标定

（疑问：波段1-9的辐射标定与反射标定，波段10-11的辐射标定与亮度温度标定为什么基本一致）

4. 利用辅助数据对图像数据进行地理定位（编程完成）

1）简介

遥感图像的地理定位：是将遥感图像中的像素与地球表面的实际地理位置相对应的过程。

2）编程实现

找到角点坐标（MTL文件中介绍）>>转换为WGS 坐标格式（ UTM_ZONE = 51）>>建立坐标系>>投射坐标。

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from pyproj import Proj, transform as pyproj_transform
from rasterio.plot import show
# 文件路径
file_path = 'D:/YANJIUSHENG/南航博士申请探索题/Dataset/data/LC08_L1TP_118038_20201222_20210310_01_T1_B1.TIF'

# 定义UTM和WGS84坐标系统
utm_proj = Proj(init='epsg:32651')   # UTM zone 
wgs84_proj = Proj(init='epsg:4326')  # WGS84

# 定义四个角点的UTM坐标
utm_corners = [(282900, 3629700), (514200, 3629700), (282900, 3394200), (514200, 3394200)]

# 转换角点坐标为WGS84
wgs84_corners = [pyproj_transform(utm_proj, wgs84_proj, pt[0], pt[1]) for pt in utm_corners]

# 计算WGS84坐标的边界
left, top = wgs84_corners[0]
right, bottom = wgs84_corners[3]

# 读取图像
with rasterio.open(file_path) as src:
    # 读取图像数据
    data = src.read(1)

    # 为WGS84坐标创建新的仿射变换
    transform = from_bounds(left, bottom, right, top, data.shape[1], data.shape[0])

    # 设置图像显示
    fig, ax = plt.subplots(figsize=(6, 4))
    show(data, ax=ax, transform=transform, cmap='gray', extent=(left, right, bottom, top))
    ax.set_title("Image in WGS84 Coordinate System")
    plt.show()

3) 结果展示

和原文件中四个角点坐标（32.78382，120.68188），（32.80518，123.15168），（30.66079，120.73416），（30.68047，123.14825）是一致的。

5. 通过几何变换，利用该图像数据仿真另一观测几何条件下的图像数据（观测区域不变，观测角度和分辨率发生变化，编程完成）

1）概述

a. 几何变换定义：

几何变换又称空间变换，就是将一组图像数据经过某种数学运算，映射成另外一组图像数据的操作。也就是将一幅图像中的坐标位置映射到另一幅图像中的新坐标位置。或者说是不改变图像的像素值，在图像平面上对像素进行重新安排。

几何变换的关键就是要确定这种空间映射关系。

b. 为什么要进行几何变换：

1）消除图像由于角度、透视关系、拍摄等原因造成的几何失真进而造成的计算机模型或者算法无法正确识别图像的现象，对于遥感图像主要是消除地球曲率、卫星轨道和姿态、大气折射等因素的影响；

2）为了更好地理解或展示地表特征，可能需要从不同的视角查看图像。几何变换可以模拟从不同角度观察地球表面的视图。

3）在多时相或多源的遥感分析中，需要将不同时间或不同来源的图像精确叠加。几何变换用于确保图像之间的对应点在空间上对齐，这对于变化检测、多光谱分析等应用至关重要。

4）在进行多源数据融合、三维建模或景观分析时，需要对来自不同角度和不同传感器的图像进行准确的几何校正，以确保数据的一致性和准确性。

2）几何变换常见方法（部分内容参考了opencv学习笔记（五）：图像几何变换 - 知乎 (zhihu.com)）

为避免大量运算，这里裁剪了原图的一小部分；为了成像清晰，考虑到城市地区，融合了7、6、4波段对应的相同部分（相关代码见代码文档，基本同前）。原图如下：

a. 视觉效果分类

缩放：对图像的大小进行放大和缩小的变换。

def apply_zoom(image, zoom_factor):
      center_x, center_y = image.shape[1] / 2, image.shape[0] / 2 #计算中心点，用作缩放的枢轴
      matrix = cv2.getRotationMatrix2D((center_x, center_y), 0, zoom_factor) #变换矩阵，包括旋转中心、旋转角度、缩放因子三个参数
      return cv2.warpAffine(image, matrix, (image.shape[1], image.shape[0]))

平移

1
2
3

def translation(image, tx, ty):
    matrix = np.float32([[1, 0, tx], [0, 1, ty]])# [1, 0, tx]表示水平方向上不改变像素的位置，但整体向右移动 tx 单位，第二行相反
    return cv2.warpAffine(image, matrix, (image.shape[1], image.shape[0]))

旋转

def rotation(image, angle):
    center_x, center_y = image.shape[1] / 2, image.shape[0] / 2
    matrix = cv2.getRotationMatrix2D((center_x, center_y), angle, 1)#与缩放一样
    return cv2.warpAffine(image, matrix, (image.shape[1], image.shape[0]))

镜像

1
2

def mirror(image):
    return cv2.flip(image, 1)# cv2.flip(img, flipcode),flipcode: flipcode=0 表示绕x轴翻转；flipcode=任意正整数，比如1，2，3，表示绕y轴翻转；flipcode=任意负整数，比如-1，-2，-3，表示绕x轴和y轴同时翻转；

结果展示：

b. 数学原理分类

首先需要了解几种变换矩阵，理解图像的坐标变换：

仿射变换（Affine Transformation）

定义：仿射变换通俗的讲就是就是“线性变换”+“平移”。上述的缩放、平移等操作都可以称为仿射变换。其常用于图像的校正、拼接、配准等操作中。

在二维图像处理中，仿射变换可以理解为一种可以改变图像形状和大小的变换，同时保持图像中平行线的平行性和直线的直线性(点还是点，线还是线，面还是面，点变不成线，线变不成面，同一条直线上的点的位置和长度比例关系不变，变的是向量的夹角、垂直的关系会发生变换）。比如下图：

该图参考了(1 封私信) 如何通俗地讲解「仿射变换」？ - 知乎 (zhihu.com)

仿射变换

如何实现：一般来讲，仿射变换主要基于三种方法，一是在相同坐标空间下，点线面自身移动；二是点线面本身没有移动，而坐标空间发生了变换，比如坐标轴发生了缩放、旋转或者平移导致的点线面坐标值发生变换，这样点线面也就相当于移动了；三是，点线面以及坐标轴本身都发生了变换，但相对的量不等，但也相当于移动。公式如下：

1704355371902

其中，a_11,a_12,a_21,a_22 定义了图像的旋转、缩放和倾斜，t_x,t_y 定义了平移量。(x,y)为原图像中的点。（x’,y’)为变换后图像中的点。

代码：

def affine_transform(image):
    rows, cols = image.shape[:2]
    src_points = np.float32([[0, 0], [cols-1, 0], [0, rows-1]])#原图的三个点
    dst_points = np.float32([[0,image.shape[1]*0.1], [image.shape[1]*0.85,image.shape[0]*0.25], [image.shape[1]*0.15, image.shape[0]*0.7]])
    #仿射变换矩阵，可以是原图的几倍也可以是下面这种固定的坐标dst_points = np.float32([[10,10],[200,120],[50,150]])
    matrix = cv2.getAffineTransform(src_points, dst_points) #仿射变换矩阵
    return cv2.warpAffine(image, matrix, (cols, rows))

结果展示：

仿射变换结果图

透视变换（Perspective Transformation）

定义：把一个图像投影到一个新的视平面的过程。其常用于畸变校正、3D视图的模拟等操作中。

透视变换能够改变图像视点，模拟从不同角度和不同距离观察场景的效果，这是通过改变图像中对象的透视来实现的。透视变换是一个非线性变换过程直线仍然映射为直线，但平行线可能不再保持平行。

该图参考了（十四）透视变换_淡定的炮仗的博客-CSDN博客

透视变换

透视变换包含了仿射变换。

透视变换区别

此外还有逆透视变换，即将一个具有透视效果的图像转换回其原始的视图。常用于摄像头捕获的图像中校正透视畸变，或者将倾斜拍摄的图像矫正为正面视图。

如何实现：该过程包括把一个二维坐标系转换为三维坐标系，然后把三维坐标系投影到新的二维坐标系。与仿射变换不同（仿射变换只需要三个点，因为平行线仍然平行，线的比例相同，确定三条线后即可确定最后一条线），透视变换需要至少四对对应点（原图像中的点和目标图像中的点）来计算变换矩阵。公式如下：

1704355356464

其中（x,y)是原始图像中的点，(x*′,*y′) 是变换后的点，而 w′ 是齐次坐标，用于保持变换的线性表达。在实际应用中，我们通常对结果进行归一化，即 (x’/w’,y’/w’)，以得到最终的图像坐标。变换矩阵中h_00,h_01,h_10,h_11用于线性变换，h_02,h_23用于产生图像透视变换，h_20,h_21用于图像的平移。

代码：

def perspective_transform(image):
    rows, cols = image.shape[:2]
    src_points = np.float32([[0, 0], [cols-1, 0], [0, rows-1], [cols-1, rows-1]])
    dst_points = np.float32([[0, 0], [int(0.6*cols), rows//3], [int(0.4*cols), rows], [cols-1, rows-2*rows//3]])
    matrix = cv2.getPerspectiveTransform(src_points, dst_points)
    return cv2.warpPerspective(image, matrix, (cols, rows))

结果展示：

透视变换结果

**重映射变换(Remap Transformation) **

定义：把一张图像内的像素点放置到另外一幅图像内指定的位置。其常用于图像扭曲、畸变矫正等操作中。

上述的两周变换都是通过变换矩阵来实现，重映射变换则是通过自定义的方式进行映射。

如何实现： 在重映射操作中，输出图像中的每个像素位置都通过映射规则从输入图像的相应位置上取得像素值。这个映射规则由两个独立的映射函数定义，分别是 x*′=f(x,y) 和 y′=g(x,*y)，其中 *(x,*y) 是输入图像中的原始坐标，而(x’,y’) 是输出图像中的新坐标。

代码：

def remap(image):
    rows, cols = image.shape[:2]
    map_x, map_y = np.float32(np.indices((rows, cols)))
    map_x = np.sin(map_x/15)*15 + map_x#定义变换函数
    map_y = np.cos(map_y/20)*10 + map_y
    return cv2.remap(image, map_x, map_y, cv2.INTER_LINEAR)

结果展示：

重射变换

3）文献扩展

弹性变换（elastic_transform)

本文旨在阐述神经网络在文档视觉输入分类中的有效应用（总的任务）。论文指出，尽管神经网络方法繁多且令人困惑，但仍有一些具体的最佳实践可供文档分析研究者使用以获得良好结果。论文强调了扩大训练集的重要性，并提出了一种新的数据扭曲方式——弹性变换（方法1）。此外，论文还提出卷积神经网络比全连接网络更适合视觉文档任务，并描述了一个简单但通用的卷积神经网络架构（方法2）。论文得出结论，使用弹性变形和卷积神经网络在MNIST数据集上取得了最高性能。这些结果反映了两个重要问题：训练集的大小和质量是学习系统质量的主要决定因素，卷积神经网络利用输入数据的空间结构，不同于其他依赖独立元素的分类技术（结论）。

围绕本课题主题，这里只阐述方法1的实现，即弹性变换：

创建一个随机的位移场使图像变形，如下式所示：

其中，rand(-1,+1)是生成一个在（-1，1）之间的随机数。

利用高斯分布的N（0,σ）生成的随机位移场进行卷积操作(和CNN中的卷积操作一样，说白了就是滤波操作)。σ越大产生的图像越平滑。下图是论文中的不同σ值对随机位移场的影响。

用一个控制因子α与随机位移场相乘，用以控制其变形强度，得到一个弹性形变的位移场。
弹性形变的位移场运用到原图的仿射变换的图像上，得到最终结果。

代码：

def elastic_transform(image, alpha, sigma, random_state=None): #四个参数：原图片、两个控制强度的参数（ alpha 越大、sigma越小强度越大）、随机变换的参数（可设置）
    if random_state is None:
        random_state = np.random.RandomState(None) #随机变换

    shape = image.shape #变换后的图像
    dx = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma, mode="constant", cval=0) * alpha #水平方向上的扭曲函数
    dy = gaussian_filter((random_state.rand(*shape) * 2 - 1), sigma, mode="constant", cval=0) * alpha #垂直方向上的扭曲函数

    x, y = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))    #定义新的一个网格数组
    indices = np.reshape(y + dy, (-1, 1)), np.reshape(x + dx, (-1, 1)) #将变换后的图像加到原像素上

    distored_image = map_coordinates(image, indices, order=1).reshape(shape) #创建图像
    return distored_image

结果：

参考文献：

[1] P. Y. Simard, S. David, J. C. Platt, “Best practices for convolutional neural networks applied to visual document analysis”. Icdar. Vol. 3. No. 2003. 2003.

6.根据英文科技论文的写作要求，用1小节对第5题的方法进行描述

以透视变换为例：

Perspective transformation is a method of image geometric transformation that projects an image onto a new viewing plane while ensuring that lines within the image remain mapped as lines [1]. As illustrated in Fig. 1, each pixel point of the original two-dimensional image undergoes projection through a 3x3 transformation matrix into a three-dimensional space. Subsequently, a normalization factor is applied to re-project this three-dimensional image back into a two-dimensional coordinate system, resulting in the final perspective transformed image.

Fig. 1. Perspective transformation schematic diagram

The entire transformation process can be represented by Eq. (1),

Where, (x, y) represents a point in the original image, (x’, y’) denotes the transformed point, w’ is a factor used for normalization. The elements of the perspective transformation matrix are denoted as h_ij. The components h₀₀, h₀₁, h₁₀, h₁₁ are involved in the linear transformation of the image, the elements h₀₂, h₁₂ are responsible for generating the perspective transformation of the image, the elements h₂₀, h₂₁ are used for the translation of the image.

Reference:

[1] F.S.Hill，JR. Computer Graphics Using OpenGL, Second Edition, Science Press and Pearson Education, 2004.

7. 根据中文科技论文的写作要求，用1小节对第5题的方法进行描述（注意中英文语言的写作区别）

透射变换（Perspective Transformation）是一种将图像投影到一个新的视平面，并保证图像中直线仍然映射为直线的图像几何变换方法，计算简便、应用灵活，能够根据预设的参数调整图像的视角和比例实现对图像的校正和增强[1]。在透视变换中，图像中的每一个像素点都会通过一个透视变换矩阵从原始坐标系统中转换到新的坐标系统下。透视变换的公式可以表示为：

式中，(x,y) 是原始图像中的点，(x‘,y’）是变换后的点，w’是用于归一化的因子；h_ij为透视变换矩阵中的元素，h₀₀，h₀₁，h₁₀，h₁₁ 用于图像的线性变换，h₀₂ ，h₁₂用于产生图像透视变换，h₂₀，h₂₁用于图像的平移。

通过适当选择变换矩阵中的元素，透视变换实现了从一种视角到另一种视角的平滑过渡，使得观看者能够从不同的视角观察同一个场景或对象。

参考文献：

[1] F.S.Hill，JR. Computer Graphics Using OpenGL, Second Edition, Science Press and Pearson Education, 2004.

中英文对比：

英文可能更侧重方法的复现性，需要对方法进行详细的描述，语言更直接。中文可能更侧重行文的逻辑性，比如这是一个在怎样的方法，主要怎么实现，实现的结果如何。

博士探索题

题1 通过网络等手段找到任意某颗遥感卫星的一幅**图像数据和辅助数据**，完成以下工作：

1.对图像数据和辅助数据进行解译，解释说明图像的各个参数含义（编程完成）

1).数据来源

2).数据解读（下述内容部分参考【Landsat 8】遥感影像文件内容及命名规则_landsat8影像命名规则_Uyoin的博客-CSDN博客

a. 文件构成

b. 文件解释

3).解译过程 （下述部分内容参考Python 处理遥感影像入门 - 知乎 (zhihu.com)）

a. 导入必须库

b. 打印文件信息

c. 显示图像

d. 波段组合

e. 计算NDVI（归一化植被制数）

f. 更多指数

g. 保存结果(NDVI)到本地

h. 空间索引及提取像素值

2. 利用STK仿真载荷成像时刻轨道位置与视场

1）联合python与STK创建应用获取接口

2）创建场景

3）设置场景时间

4) 创建卫星，设置参数

5）结果展示

3. 利用辅助数据对图像数据进行辐射定标（编程完成）

1）简介

a. 辐射定标定义

b. 为什么要进行辐射定位

c. 还可以进行哪些定标

2) 编程实现

a. 在TXT文件中找到相关参数

b. 代码

3）结果展示

4. 利用辅助数据对图像数据进行地理定位（编程完成）

1）简介

2）编程实现

3) 结果展示

5. 通过几何变换，利用该图像数据仿真另一观测几何条件下的图像数据（观测区域不变，观测角度和分辨率发生变化，编程完成）

1）概述

a. 几何变换定义：

b. 为什么要进行几何变换：

2）几何变换常见方法（部分内容参考了opencv学习笔记（五）：图像几何变换 - 知乎 (zhihu.com)）

a. 视觉效果分类

b. 数学原理分类

3）文献扩展

6.根据英文科技论文的写作要求，用1小节对第5题的方法进行描述

7. 根据中文科技论文的写作要求，用1小节对第5题的方法进行描述（注意中英文语言的写作区别）

题1 通过网络等手段找到任意某颗遥感卫星的一幅图像数据和辅助数据，完成以下工作：

3).解译过程（下述部分内容参考Python 处理遥感影像入门 - 知乎 (zhihu.com)）