目录
一、landsat8数据介绍
1）数据下载
2）Landsat8中*_MTL.txt文件详解
GROUP = LANDSAT_METADATA_FILE
#*************************************************************************************************#
#该部分内容是关于元数据文件（即MTL.tXt文件本身）的数据信息的介绍
GROUP = PRODUCT_CONTENTS
#表明影像数据的来源是美国USGS
ORIGIN = "Image courtesy of the U.S. Geological Survey"
DIGITAL_OBJECT_IDENTIFIER = "https://doi.org/10.5066/P9OGBGM6"
“""
这是数据产品（包含但不仅限于影像本身）的ID号，其中主要包含三个信息；第一是数据产品是L1一级产品数据；
第二是数据中影像拍摄的时间为2021-02-12第三是对影像进行处理后制成的数据产品（即数据产品生产的时间）是2021-03-02；
第三是该数据产品是T1类型（稍后会讲解T1表示什么含义
“""
LANDSAT_PRODUCT_ID = "LC08_L2SP_130039_20210212_20210302_02_T1"
PROCESSING_LEVEL = "L2SP"
“""
“Collection 1”与“Collection 2”是USGS对Landsat数据进行的两次不同的处理所形成的采用分级结构管理的数据集。
1）相较于Collection 1，Collection 2 Level 1数据的辐射定标与几何校正精度有所提高；
（2）Collection 2 处理和分发 了Landsat 4-5 TM、Landsat 7 ETM+ 和 Landsat 8 OLI/TIRS 的 Level-2 地表反射率（SR）和地表温度（ST）产品。
“""
COLLECTION_NUMBER = 02
“""
“T”即 Tier，官方全称叫“Landsat Collection Tier”，可以理解为数据等级，针对Level-1数据，基于不同数据质量和处理级别划定不同等级。
Landsat Collection Tiers are the inventory structure for Level-1 data products and are based on data quality and level of processing.
T1（Tier 1）：
具有最高数据质量的Landsat影像被放入T1，并被认为适合于时间序列分析。
T1存放的是L1TP（Level-1 Precision and Terrain）处理等级的数据。
L1TP等级数据包括辐射定标和使用地面控制点 (GCP) 和数字高程模型 (DEM) 数据进行的正射校正（正射校正是几何校正的最高级别），以校正地形位移，所采用的地面控制点来自 Global Land Survey 2000 (GLS2000) 数据集。
该等级数据具有良好的辐射特征，并且在不同的 Landsat 传感器之间进行了相互校准。T1等级影像的地理配准是一致的，校正的均方根误差（RMSE）在12米之内。
T2（Tier 2）：
在处理过程中不符合T1标准（主要指几何校正精度）的 Landsat 影像被分配到T2。
T2等级数据的辐射定标与T1是一致的，但由于轨道信息不准确，无法达到T1等级的几何校正精度。存放的是L1GT和L1GS处理等级的数据。
L1GT等级数据包括辐射定标和使用航天器星历数据和 DEM 数据应用系统几何校正进行辐射校准，以校正地形位移。
L1GS等级数据包括辐射定标和仅使用航天器星历数据进行系统几何校正。
T2等级数据没有明确的使用场景，官方给出的建议是可以根据 RMSE 等属性来确定。
RT（Real-Time）：
RT是实时数据，目前在役的 Landsat 7 ETM+ 和 Landsat 8 OLI/TIRS 数据，获取但尚未处理的临时数据存放在RT中。数据处理之后就会放到T1或T2中，并从RT中删除。
RT数据存在的目的是能够在需要紧急响应的情况下快速分发经过有限校准的数据。
“""
COLLECTION_CATEGORY = "T1"
#影像的输出文件类型为Geottiff因此你可以看到各个波段的后缀均为.tif，
OUTPUT_FORMAT = "GEOTIFF"
#文件夹种各个文件的命名
FILE_NAME_BAND_1 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B1.TIF"
FILE_NAME_BAND_2 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B2.TIF"
FILE_NAME_BAND_3 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B3.TIF"
FILE_NAME_BAND_4 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B4.TIF"
FILE_NAME_BAND_5 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B5.TIF"
FILE_NAME_BAND_6 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B6.TIF"
FILE_NAME_BAND_7 = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_B7.TIF"
FILE_NAME_BAND_ST_B10 = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_B10.TIF"
FILE_NAME_THERMAL_RADIANCE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_TRAD.TIF"
FILE_NAME_UPWELL_RADIANCE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_URAD.TIF"
FILE_NAME_DOWNWELL_RADIANCE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_DRAD.TIF"
FILE_NAME_ATMOSPHERIC_TRANSMITTANCE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_ATRAN.TIF"
FILE_NAME_EMISSIVITY = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_EMIS.TIF"
FILE_NAME_EMISSIVITY_STDEV = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_EMSD.TIF"
FILE_NAME_CLOUD_DISTANCE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_CDIST.TIF"
FILE_NAME_QUALITY_L2_AEROSOL = "LC08_L2SP_130039_20210212_20210302_02_T1_SR_QA_AEROSOL.TIF"
FILE_NAME_QUALITY_L2_SURFACE_TEMPERATURE = "LC08_L2SP_130039_20210212_20210302_02_T1_ST_QA.TIF"
FILE_NAME_QUALITY_L1_PIXEL = "LC08_L2SP_130039_20210212_20210302_02_T1_QA_PIXEL.TIF"
FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION = "LC08_L2SP_130039_20210212_20210302_02_T1_QA_RADSAT.TIF"
FILE_NAME_ANGLE_COEFFICIENT = "LC08_L2SP_130039_20210212_20210302_02_T1_ANG.txt"
FILE_NAME_METADATA_ODL = "LC08_L2SP_130039_20210212_20210302_02_T1_MTL.txt"
FILE_NAME_METADATA_XML = "LC08_L2SP_130039_20210212_20210302_02_T1_MTL.xml"
#数据类型
DATA_TYPE_BAND_1 = "UINT16"
DATA_TYPE_BAND_2 = "UINT16"
DATA_TYPE_BAND_3 = "UINT16"
DATA_TYPE_BAND_4 = "UINT16"
DATA_TYPE_BAND_5 = "UINT16"
DATA_TYPE_BAND_6 = "UINT16"
DATA_TYPE_BAND_7 = "UINT16"
DATA_TYPE_BAND_ST_B10 = "UINT16"
DATA_TYPE_THERMAL_RADIANCE = "INT16"
DATA_TYPE_UPWELL_RADIANCE = "INT16"
DATA_TYPE_DOWNWELL_RADIANCE = "INT16"
DATA_TYPE_ATMOSPHERIC_TRANSMITTANCE = "INT16"
DATA_TYPE_EMISSIVITY = "INT16"
DATA_TYPE_EMISSIVITY_STDEV = "INT16"
DATA_TYPE_CLOUD_DISTANCE = "INT16"
DATA_TYPE_QUALITY_L2_AEROSOL = "UINT8"
DATA_TYPE_QUALITY_L2_SURFACE_TEMPERATURE = "INT16"
DATA_TYPE_QUALITY_L1_PIXEL = "UINT16"
DATA_TYPE_QUALITY_L1_RADIOMETRIC_SATURATION = "UINT16"
END_GROUP = PRODUCT_CONTENTS
GROUP = IMAGE_ATTRIBUTES
# 卫星为Landsat8
SPACECRAFT_ID = "LANDSAT_8"
#卫星传感器为OLlTlRS；
SENSOR_ID = "OLI_TIRS"
WRS_TYPE = 2
WRS_PATH = 130
WRS_ROW = 39
NADIR_OFFNADIR = "NADIR"
TARGET_WRS_PATH = 130
TARGET_WRS_ROW = 39
DATE_ACQUIRED = 2021-02-12
SCENE_CENTER_TIME = "03:39:39.6572859Z"
STATION_ID = "LGN"
CLOUD_COVER = 2.24
CLOUD_COVER_LAND = 2.24
IMAGE_QUALITY_OLI = 9
IMAGE_QUALITY_TIRS = 9
SATURATION_BAND_1 = "N"
SATURATION_BAND_2 = "N"
SATURATION_BAND_3 = "N"
SATURATION_BAND_4 = "N"
SATURATION_BAND_5 = "N"
SATURATION_BAND_6 = "Y"
SATURATION_BAND_7 = "Y"
SATURATION_BAND_8 = "N"
SATURATION_BAND_9 = "N"
ROLL_ANGLE = -0.001
######################################################
影响拍摄时的太阳方位角；太阳高度角；日地距离
SUN_AZIMUTH = 146.94993267
SUN_ELEVATION = 39.60545988
EARTH_SUN_DISTANCE = 0.9871722
######################################################
TRUNCATION_OLI = "UPPER"
TIRS_SSM_MODEL = "FINAL"
TIRS_SSM_POSITION_STATUS = "ESTIMATED"
END_GROUP = IMAGE_ATTRIBUTES
GROUP = PROJECTION_ATTRIBUTES
MAP_PROJECTION = "UTM"
DATUM = "WGS84"
ELLIPSOID = "WGS84"
UTM_ZONE = 48
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
REFLECTIVE_LINES = 7881
REFLECTIVE_SAMPLES = 7741
THERMAL_LINES = 7881
THERMAL_SAMPLES = 7741
ORIENTATION = "NORTH_UP"
# 这是影像的角点信息，UL表示UPLeft表示左上角点，LR表示LowRight右下角点，Lan表示纬度，Lon表示经度；
CORNER_UL_LAT_PRODUCT = 31.34382
CORNER_UL_LON_PRODUCT = 101.77545
CORNER_UR_LAT_PRODUCT = 31.38196
CORNER_UR_LON_PRODUCT = 104.21544
CORNER_LL_LAT_PRODUCT = 29.21381
CORNER_LL_LON_PRODUCT = 101.84442
CORNER_LR_LAT_PRODUCT = 29.24884
CORNER_LR_LON_PRODUCT = 104.23225
#这里类似，也是关于角点信息的内容，唯一不同的是这里记录的不是经纬度信息（地理坐标），而是各个角点的投影坐标；
CORNER_UL_PROJECTION_X_PRODUCT = 193200.000
CORNER_UL_PROJECTION_Y_PRODUCT = 3472200.000
CORNER_UR_PROJECTION_X_PRODUCT = 425400.000
CORNER_UR_PROJECTION_Y_PRODUCT = 3472200.000
CORNER_LL_PROJECTION_X_PRODUCT = 193200.000
CORNER_LL_PROJECTION_Y_PRODUCT = 3235800.000
CORNER_LR_PROJECTION_X_PRODUCT = 425400.000
CORNER_LR_PROJECTION_Y_PRODUCT = 3235800.000
END_GROUP = PROJECTION_ATTRIBUTES
GROUP = LEVEL2_PROCESSING_RECORD
ORIGIN = "Image courtesy of the U.S. Geological Survey"
DIGITAL_OBJECT_IDENTIFIER = "https://doi.org/10.5066/P9OGBGM6"
REQUEST_ID = "P700ohh5vln74_00022"
LANDSAT_PRODUCT_ID = "LC08_L2SP_130039_20210212_20210302_02_T1"
PROCESSING_LEVEL = "L2SP"
OUTPUT_FORMAT = "GEOTIFF"
DATE_PRODUCT_GENERATED = 2021-03-02T13:13:02Z
PROCESSING_SOFTWARE_VERSION = "LPGS_15.4.0"
ALGORITHM_SOURCE_SURFACE_REFLECTANCE = "LaSRC_1.5.0"
DATA_SOURCE_OZONE = "MODIS"
DATA_SOURCE_PRESSURE = "Calculated"
DATA_SOURCE_WATER_VAPOR = "MODIS"
DATA_SOURCE_AIR_TEMPERATURE = "MODIS"
ALGORITHM_SOURCE_SURFACE_TEMPERATURE = "st_1.3.0"
DATA_SOURCE_REANALYSIS = "GEOS-5 FP-IT"
END_GROUP = LEVEL2_PROCESSING_RECORD
GROUP = LEVEL2_SURFACE_REFLECTANCE_PARAMETERS
REFLECTANCE_MAXIMUM_BAND_1 = 1.602213
REFLECTANCE_MINIMUM_BAND_1 = -0.199972
REFLECTANCE_MAXIMUM_BAND_2 = 1.602213
REFLECTANCE_MINIMUM_BAND_2 = -0.199972
REFLECTANCE_MAXIMUM_BAND_3 = 1.602213
REFLECTANCE_MINIMUM_BAND_3 = -0.199972
REFLECTANCE_MAXIMUM_BAND_4 = 1.602213
REFLECTANCE_MINIMUM_BAND_4 = -0.199972
REFLECTANCE_MAXIMUM_BAND_5 = 1.602213
REFLECTANCE_MINIMUM_BAND_5 = -0.199972
REFLECTANCE_MAXIMUM_BAND_6 = 1.602213
REFLECTANCE_MINIMUM_BAND_6 = -0.199972
REFLECTANCE_MAXIMUM_BAND_7 = 1.602213
REFLECTANCE_MINIMUM_BAND_7 = -0.199972
QUANTIZE_CAL_MAX_BAND_1 = 65535
QUANTIZE_CAL_MIN_BAND_1 = 1
QUANTIZE_CAL_MAX_BAND_2 = 65535
QUANTIZE_CAL_MIN_BAND_2 = 1
QUANTIZE_CAL_MAX_BAND_3 = 65535
QUANTIZE_CAL_MIN_BAND_3 = 1
QUANTIZE_CAL_MAX_BAND_4 = 65535
QUANTIZE_CAL_MIN_BAND_4 = 1
QUANTIZE_CAL_MAX_BAND_5 = 65535
QUANTIZE_CAL_MIN_BAND_5 = 1
QUANTIZE_CAL_MAX_BAND_6 = 65535
QUANTIZE_CAL_MIN_BAND_6 = 1
QUANTIZE_CAL_MAX_BAND_7 = 65535
QUANTIZE_CAL_MIN_BAND_7 = 1
REFLECTANCE_MULT_BAND_1 = 2.75e-05
REFLECTANCE_MULT_BAND_2 = 2.75e-05
REFLECTANCE_MULT_BAND_3 = 2.75e-05
REFLECTANCE_MULT_BAND_4 = 2.75e-05
REFLECTANCE_MULT_BAND_5 = 2.75e-05
REFLECTANCE_MULT_BAND_6 = 2.75e-05
REFLECTANCE_MULT_BAND_7 = 2.75e-05
REFLECTANCE_ADD_BAND_1 = -0.2
REFLECTANCE_ADD_BAND_2 = -0.2
REFLECTANCE_ADD_BAND_3 = -0.2
REFLECTANCE_ADD_BAND_4 = -0.2
REFLECTANCE_ADD_BAND_5 = -0.2
REFLECTANCE_ADD_BAND_6 = -0.2
REFLECTANCE_ADD_BAND_7 = -0.2
END_GROUP = LEVEL2_SURFACE_REFLECTANCE_PARAMETERS
GROUP = LEVEL2_SURFACE_TEMPERATURE_PARAMETERS
TEMPERATURE_MAXIMUM_BAND_ST_B10 = 372.999941
TEMPERATURE_MINIMUM_BAND_ST_B10 = 149.003418
QUANTIZE_CAL_MAXIMUM_BAND_ST_B10 = 65535
QUANTIZE_CAL_MINIMUM_BAND_ST_B10 = 1
TEMPERATURE_MULT_BAND_ST_B10 = 0.00341802
TEMPERATURE_ADD_BAND_ST_B10 = 149.0
END_GROUP = LEVEL2_SURFACE_TEMPERATURE_PARAMETERS
GROUP = LEVEL1_PROCESSING_RECORD
ORIGIN = "Image courtesy of the U.S. Geological Survey"
DIGITAL_OBJECT_IDENTIFIER = "https://doi.org/10.5066/P975CC9B"
REQUEST_ID = "P700ohh5vln74_00022"
LANDSAT_SCENE_ID = "LC81300392021043LGN00"
LANDSAT_PRODUCT_ID = "LC08_L1TP_130039_20210212_20210302_02_T1"
PROCESSING_LEVEL = "L1TP"
COLLECTION_CATEGORY = "T1"
OUTPUT_FORMAT = "GEOTIFF"
DATE_PRODUCT_GENERATED = 2021-03-02T12:55:57Z
PROCESSING_SOFTWARE_VERSION = "LPGS_15.4.0"
“""
B1：紫波段，主要用于水体研究B2:蓝波段；
B3：绿波段；
B4红波段；
B5：近红外波段；
B6B7：短波红外；
B8：全色波段；
B9：卷云波段；
B10,B11：热红外波段；
“””
FILE_NAME_BAND_1 = "LC08_L1TP_130039_20210212_20210302_02_T1_B1.TIF"
FILE_NAME_BAND_2 = "LC08_L1TP_130039_20210212_20210302_02_T1_B2.TIF"
FILE_NAME_BAND_3 = "LC08_L1TP_130039_20210212_20210302_02_T1_B3.TIF"
FILE_NAME_BAND_4 = "LC08_L1TP_130039_20210212_20210302_02_T1_B4.TIF"
FILE_NAME_BAND_5 = "LC08_L1TP_130039_20210212_20210302_02_T1_B5.TIF"
FILE_NAME_BAND_6 = "LC08_L1TP_130039_20210212_20210302_02_T1_B6.TIF"
FILE_NAME_BAND_7 = "LC08_L1TP_130039_20210212_20210302_02_T1_B7.TIF"
FILE_NAME_BAND_8 = "LC08_L1TP_130039_20210212_20210302_02_T1_B8.TIF"
FILE_NAME_BAND_9 = "LC08_L1TP_130039_20210212_20210302_02_T1_B9.TIF"
FILE_NAME_BAND_10 = "LC08_L1TP_130039_20210212_20210302_02_T1_B10.TIF"
FILE_NAME_BAND_11 = "LC08_L1TP_130039_20210212_20210302_02_T1_B11.TIF"
FILE_NAME_QUALITY_L1_PIXEL = "LC08_L1TP_130039_20210212_20210302_02_T1_QA_PIXEL.TIF"
FILE_NAME_QUALITY_L1_RADIOMETRIC_SATURATION = "LC08_L1TP_130039_20210212_20210302_02_T1_QA_RADSAT.TIF"
FILE_NAME_ANGLE_COEFFICIENT = "LC08_L1TP_130039_20210212_20210302_02_T1_ANG.txt"
FILE_NAME_ANGLE_SENSOR_AZIMUTH_BAND_4 = "LC08_L1TP_130039_20210212_20210302_02_T1_VAA.TIF"
FILE_NAME_ANGLE_SENSOR_ZENITH_BAND_4 = "LC08_L1TP_130039_20210212_20210302_02_T1_VZA.TIF"
FILE_NAME_ANGLE_SOLAR_AZIMUTH_BAND_4 = "LC08_L1TP_130039_20210212_20210302_02_T1_SAA.TIF"
FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4 = "LC08_L1TP_130039_20210212_20210302_02_T1_SZA.TIF"
FILE_NAME_METADATA_ODL = "LC08_L1TP_130039_20210212_20210302_02_T1_MTL.txt"
FILE_NAME_METADATA_XML = "LC08_L1TP_130039_20210212_20210302_02_T1_MTL.xml"
FILE_NAME_CPF = "LC08CPF_20210201_20210228_02.04"
FILE_NAME_BPF_OLI = "LO8BPF20210212032704_20210212050558.01"
FILE_NAME_BPF_TIRS = "LT8BPF20210205131939_20210212112539.01"
FILE_NAME_RLUT = "LC08RLUT_20150303_20431231_02_01.h5"
DATA_SOURCE_TIRS_STRAY_LIGHT_CORRECTION = "TIRS"
DATA_SOURCE_ELEVATION = "GLS2000"
GROUND_CONTROL_POINTS_VERSION = 5
GROUND_CONTROL_POINTS_MODEL = 782
GEOMETRIC_RMSE_MODEL = 6.002
GEOMETRIC_RMSE_MODEL_Y = 4.252
GEOMETRIC_RMSE_MODEL_X = 4.236
GROUND_CONTROL_POINTS_VERIFY = 107
GEOMETRIC_RMSE_VERIFY = 3.323
END_GROUP = LEVEL1_PROCESSING_RECORD
##这里是关于各个波段最大最小DN值（像元值）所对应的辐射亮度；其主要用于辐射定标
GROUP = LEVEL1_MIN_MAX_RADIANCE
RADIANCE_MAXIMUM_BAND_1 = 779.94421
RADIANCE_MINIMUM_BAND_1 = -64.40805
RADIANCE_MAXIMUM_BAND_2 = 798.67242
RADIANCE_MINIMUM_BAND_2 = -65.95463
RADIANCE_MAXIMUM_BAND_3 = 735.96979
RADIANCE_MINIMUM_BAND_3 = -60.77663
RADIANCE_MAXIMUM_BAND_4 = 620.61121
RADIANCE_MINIMUM_BAND_4 = -51.25027
RADIANCE_MAXIMUM_BAND_5 = 379.78311
RADIANCE_MINIMUM_BAND_5 = -31.36261
RADIANCE_MAXIMUM_BAND_6 = 94.44860
RADIANCE_MINIMUM_BAND_6 = -7.79960
RADIANCE_MAXIMUM_BAND_7 = 31.83424
RADIANCE_MINIMUM_BAND_7 = -2.62888
RADIANCE_MAXIMUM_BAND_8 = 702.36108
RADIANCE_MINIMUM_BAND_8 = -58.00121
RADIANCE_MAXIMUM_BAND_9 = 148.42784
RADIANCE_MINIMUM_BAND_9 = -12.25722
RADIANCE_MAXIMUM_BAND_10 = 22.00180
RADIANCE_MINIMUM_BAND_10 = 0.10033
RADIANCE_MAXIMUM_BAND_11 = 22.00180
RADIANCE_MINIMUM_BAND_11 = 0.10033
END_GROUP = LEVEL1_MIN_MAX_RADIANCE
"""
这个也是类似，各个波段的最大最小DN值所对应的表观反射率(大气层顶端的反射率）
如果你是进行在进行完辐射定标之后需要进行大气校正，并且打算使用ENVI的Flassh模块进行大气校正，
那么你不应该使用反射率进行辐射定标，因为Flassh模型的运行要求得到的辐射定标文件应是辐射亮度而非表观发射率。
"""
GROUP = LEVEL1_MIN_MAX_REFLECTANCE
REFLECTANCE_MAXIMUM_BAND_1 = 1.210700
REFLECTANCE_MINIMUM_BAND_1 = -0.099980
REFLECTANCE_MAXIMUM_BAND_2 = 1.210700
REFLECTANCE_MINIMUM_BAND_2 = -0.099980
REFLECTANCE_MAXIMUM_BAND_3 = 1.210700
REFLECTANCE_MINIMUM_BAND_3 = -0.099980
REFLECTANCE_MAXIMUM_BAND_4 = 1.210700
REFLECTANCE_MINIMUM_BAND_4 = -0.099980
REFLECTANCE_MAXIMUM_BAND_5 = 1.210700
REFLECTANCE_MINIMUM_BAND_5 = -0.099980
REFLECTANCE_MAXIMUM_BAND_6 = 1.210700
REFLECTANCE_MINIMUM_BAND_6 = -0.099980
REFLECTANCE_MAXIMUM_BAND_7 = 1.210700
REFLECTANCE_MINIMUM_BAND_7 = -0.099980
REFLECTANCE_MAXIMUM_BAND_8 = 1.210700
REFLECTANCE_MINIMUM_BAND_8 = -0.099980
REFLECTANCE_MAXIMUM_BAND_9 = 1.210700
REFLECTANCE_MINIMUM_BAND_9 = -0.099980
END_GROUP = LEVEL1_MIN_MAX_REFLECTANCE
#各个波段的影像的最大最小DN值，就是像元属性值的最大最小值；
GROUP = LEVEL1_MIN_MAX_PIXEL_VALUE
QUANTIZE_CAL_MAX_BAND_1 = 65535
QUANTIZE_CAL_MIN_BAND_1 = 1
QUANTIZE_CAL_MAX_BAND_2 = 65535
QUANTIZE_CAL_MIN_BAND_2 = 1
QUANTIZE_CAL_MAX_BAND_3 = 65535
QUANTIZE_CAL_MIN_BAND_3 = 1
QUANTIZE_CAL_MAX_BAND_4 = 65535
QUANTIZE_CAL_MIN_BAND_4 = 1
QUANTIZE_CAL_MAX_BAND_5 = 65535
QUANTIZE_CAL_MIN_BAND_5 = 1
QUANTIZE_CAL_MAX_BAND_6 = 65535
QUANTIZE_CAL_MIN_BAND_6 = 1
QUANTIZE_CAL_MAX_BAND_7 = 65535
QUANTIZE_CAL_MIN_BAND_7 = 1
QUANTIZE_CAL_MAX_BAND_8 = 65535
QUANTIZE_CAL_MIN_BAND_8 = 1
QUANTIZE_CAL_MAX_BAND_9 = 65535
QUANTIZE_CAL_MIN_BAND_9 = 1
QUANTIZE_CAL_MAX_BAND_10 = 65535
QUANTIZE_CAL_MIN_BAND_10 = 1
QUANTIZE_CAL_MAX_BAND_11 = 65535
QUANTIZE_CAL_MIN_BAND_11 = 1
END_GROUP = LEVEL1_MIN_MAX_PIXEL_VALUE
GROUP = LEVEL1_RADIOMETRIC_RESCALING
#辐射定标参数
#for lines in range(parameter_start_line, parameter_start_line + 11):
#multi_parameters
RADIANCE_MULT_BAND_1 = 1.2884E-02
RADIANCE_MULT_BAND_2 = 1.3194E-02
RADIANCE_MULT_BAND_3 = 1.2158E-02
RADIANCE_MULT_BAND_4 = 1.0252E-02
RADIANCE_MULT_BAND_5 = 6.2738E-03
RADIANCE_MULT_BAND_6 = 1.5602E-03
RADIANCE_MULT_BAND_7 = 5.2588E-04
RADIANCE_MULT_BAND_8 = 1.1603E-02
RADIANCE_MULT_BAND_9 = 2.4519E-03
RADIANCE_MULT_BAND_10 = 3.3420E-04
RADIANCE_MULT_BAND_11 = 3.3420E-04
# 由于RADIANCE_MULT_BAND参数和RADIANCE_ADD_BAND参数是挨着的，所以直接+11个参数即可
# add_parameters
RADIANCE_ADD_BAND_1 = -64.42093
RADIANCE_ADD_BAND_2 = -65.96782
RADIANCE_ADD_BAND_3 = -60.78878
RADIANCE_ADD_BAND_4 = -51.26053
RADIANCE_ADD_BAND_5 = -31.36889
RADIANCE_ADD_BAND_6 = -7.80116
RADIANCE_ADD_BAND_7 = -2.62941
RADIANCE_ADD_BAND_8 = -58.01281
RADIANCE_ADD_BAND_9 = -12.25967
RADIANCE_ADD_BAND_10 = 0.10000
RADIANCE_ADD_BAND_11 = 0.10000
REFLECTANCE_MULT_BAND_1 = 2.0000E-05
REFLECTANCE_MULT_BAND_2 = 2.0000E-05
REFLECTANCE_MULT_BAND_3 = 2.0000E-05
REFLECTANCE_MULT_BAND_4 = 2.0000E-05
REFLECTANCE_MULT_BAND_5 = 2.0000E-05
REFLECTANCE_MULT_BAND_6 = 2.0000E-05
REFLECTANCE_MULT_BAND_7 = 2.0000E-05
REFLECTANCE_MULT_BAND_8 = 2.0000E-05
REFLECTANCE_MULT_BAND_9 = 2.0000E-05
REFLECTANCE_ADD_BAND_1 = -0.100000
REFLECTANCE_ADD_BAND_2 = -0.100000
REFLECTANCE_ADD_BAND_3 = -0.100000
REFLECTANCE_ADD_BAND_4 = -0.100000
REFLECTANCE_ADD_BAND_5 = -0.100000
REFLECTANCE_ADD_BAND_6 = -0.100000
REFLECTANCE_ADD_BAND_7 = -0.100000
REFLECTANCE_ADD_BAND_8 = -0.100000
REFLECTANCE_ADD_BAND_9 = -0.100000
END_GROUP = LEVEL1_RADIOMETRIC_RESCALING
GROUP = LEVEL1_THERMAL_CONSTANTS
K1_CONSTANT_BAND_10 = 774.8853
K2_CONSTANT_BAND_10 = 1321.0789
K1_CONSTANT_BAND_11 = 480.8883
K2_CONSTANT_BAND_11 = 1201.1442
END_GROUP = LEVEL1_THERMAL_CONSTANTS
GROUP = LEVEL1_PROJECTION_PARAMETERS
MAP_PROJECTION = "UTM"
DATUM = "WGS84"
ELLIPSOID = "WGS84"
UTM_ZONE = 48
GRID_CELL_SIZE_PANCHROMATIC = 15.00
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
ORIENTATION = "NORTH_UP"
RESAMPLING_OPTION = "CUBIC_CONVOLUTION"
END_GROUP = LEVEL1_PROJECTION_PARAMETERS
END_GROUP = LANDSAT_METADATA_FILE
END
二、代码
代码地址:项目首页 - GitCode
1）环境安装
conda install gdal
conda install -c conda-forge py6s
2）代码分析
辐射定标：辐射亮度L=DN*Gain+Bias
大气校正：表观反射率ρ=π*L*D²/（ESUN*cosθ）式中，ρ为大气层顶( TOA) 表观反射率( 无量纲)，π为常量( 球面度sr)，L 为大气层顶进人卫星传感器的光谱辐射亮度（ W*1/（m²*sr*um）），D为日地之间距离( 天文单位)， ESUN为大气层顶的平均太阳光谱辐照度（ W*1/（m²*sr*um））θ为太阳的天顶角。
#! usr/bin/env python
# -*- coding:utf-8 -*-
# created by zhaoguanhua 2017/9/25
# AtmosphericCorrection for Landsat8
import glob
import os
import sys
import tarfile
import re
#import gdal
import numpy
from Py6S import *
from osgeo import gdal
import pdb
import shutil
import argparse
from base import MeanDEM
def parse_arguments(argv):
parser = argparse.ArgumentParser()
parser.add_argument('--Input_dir', type=str, help='Input dir', default=r'/media/cyz/gaoguangpu/LC08_L2SP_130039_20210212_20210302_02_T1')
parser.add_argument('--Output_dir', type=str, help='Output dir', default=r'/media/cyz/gaoguangpu/jiaozheng/AtmosphericCorrection/LC08_daqijiaozheng')
return parser.parse_args(argv)
# 逐波段辐射定标
def RadiometricCalibration(BandId):
# LandSat8 TM辐射定标参数
global data2,ImgRasterData
parameter_OLI = numpy.zeros((9,2))
#计算辐射亮度参数
parameter_OLI[0,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,0] = float(''.join(re.findall('RADIANCE_MULT_BAND_9.+',data2)).split("=")[1])
parameter_OLI[0,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_1.+',data2)[0]).split("=")[1])
parameter_OLI[1,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_2.+',data2)).split("=")[1])
parameter_OLI[2,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_3.+',data2)).split("=")[1])
parameter_OLI[3,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_4.+',data2)).split("=")[1])
parameter_OLI[4,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_5.+',data2)).split("=")[1])
parameter_OLI[5,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_6.+',data2)).split("=")[1])
parameter_OLI[6,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_7.+',data2)).split("=")[1])
parameter_OLI[7,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_8.+',data2)).split("=")[1])
parameter_OLI[8,1] = float(''.join(re.findall('RADIANCE_ADD_BAND_9.+',data2)).split("=")[1])
Gain = parameter_OLI[int(BandId) - 1,0]
Bias = parameter_OLI[int(BandId) - 1,1]
RaCal = numpy.where(ImgRasterData>0 ,Gain * ImgRasterData + Bias,-9999)
return (RaCal)
# 6s大气校正
def AtmosphericCorrection(BandId):
global data
# 6S模型
s = SixS( )
s.geometry = Geometry.User()
s.geometry.solar_z = 90-float(''.join(re.findall('SUN_ELEVATION.+',data2)).split("=")[1])
s.geometry.solar_a = float(''.join(re.findall('SUN_AZIMUTH.+',data2)).split("=")[1])
s.geometry.view_z = 0
s.geometry.view_a = 0
# 日期
Dateparm = ''.join(re.findall('DATE_ACQUIRED.+',data2)).split("=")
Date = Dateparm[1].split('-')
s.geometry.month = int(Date[1])
s.geometry.day = int(Date[2])
# 中心经纬度
point1lat = float(''.join(re.findall('CORNER_UL_LAT_PRODUCT.+',data2)).split("=")[1])
point1lon = float(''.join(re.findall('CORNER_UL_LON_PRODUCT.+',data2)).split("=")[1])
point2lat = float(''.join(re.findall('CORNER_UR_LAT_PRODUCT.+',data2)).split("=")[1])
point2lon = float(''.join(re.findall('CORNER_UR_LON_PRODUCT.+',data2)).split("=")[1])
point3lat = float(''.join(re.findall('CORNER_LL_LAT_PRODUCT.+',data2)).split("=")[1])
point3lon = float(''.join(re.findall('CORNER_LL_LON_PRODUCT.+',data2)).split("=")[1])
point4lat = float(''.join(re.findall('CORNER_LR_LAT_PRODUCT.+',data2)).split("=")[1])
point4lon = float(''.join(re.findall('CORNER_LR_LON_PRODUCT.+',data2)).split("=")[1])
#print("**********************************************************",point1lat)
sLongitude = (point1lon + point2lon + point3lon + point4lon) / 4
sLatitude = (point1lat + point2lat + point3lat + point4lat) / 4
# 大气模式类型
if sLatitude > -15 and sLatitude <= 15:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)
if sLatitude > 15 and sLatitude <= 45:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)
if sLatitude > 45 and sLatitude <= 60:
if s.geometry.month > 4 and s.geometry.month <= 9:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
else:
s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)
# 气溶胶类型大陆
s.aero_profile = AtmosProfile.PredefinedType(AeroProfile.Continental)
# 目标地物？？？？？？
s.ground_reflectance = GroundReflectance.HomogeneousLambertian(0.36)
# 550nm气溶胶光学厚度,根据日期从MODIS处获取。
#s.visibility=40.0
s.aot550 = 0.14497
# 通过研究去区的范围去求DEM高度。
pointUL = dict()
pointDR = dict()
pointUL["lat"] = point1lat
pointUL["lon"] = point1lon
pointDR["lat"] = point4lat
pointDR["lon"] = point2lon
meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001
# 研究区海拔、卫星传感器轨道高度
s.altitudes = Altitudes()
s.altitudes.set_target_custom_altitude(meanDEM)
s.altitudes.set_sensor_satellite_level()
# 校正波段（根据波段名称）
if BandId == '1':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B1)
elif BandId == '2':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B2)
elif BandId == '3':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B3)
elif BandId == '4':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B4)
elif BandId == '5':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B5)
elif BandId == '6':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B6)
elif BandId == '7':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B7)
elif BandId == '8':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B8)
elif BandId == '9':
s.wavelength = Wavelength(PredefinedWavelengths.LANDSAT_OLI_B9)
# 下垫面非均一、朗伯体
s.atmos_corr = AtmosCorr.AtmosCorrLambertianFromReflectance(-0.1)
# 运行6s大气模型
s.run()
xa = s.outputs.coef_xa
xb = s.outputs.coef_xb
xc = s.outputs.coef_xc
x = s.outputs.values
return (xa, xb, xc)
if __name__ == '__main__':
#输入数据路径
RootInputPath = parse_arguments(sys.argv[1:]).Input_dir
RootOutName = parse_arguments(sys.argv[2:]).Output_dir
#创建日志文件
LogFile = open(os.path.join(RootOutName,'log.txt'),'w')
for root,dirs,RSFiles in os.walk(RootInputPath):
#判断是否进入最底层
if len(dirs)==0:
#根据输入输出路径建立生成新文件的路径
RootInputPathList = RootInputPath.split(os.path.sep)
RootList = root.split(os.path.sep)
StartList = len(RootInputPathList)
EndList = len(RootList)
outname = RootOutName
for i in range(StartList,EndList):
if os.path.exists(os.path.join(outname,RootList[i]))==False:
os.makedirs(os.path.join(outname,RootList[i]))
outname=os.path.join(outname,RootList[i])
else:
outname=os.path.join(outname,RootList[i])
MeteDatas = glob.glob(os.path.join(root,'*MTL.txt'))
for MeteData in MeteDatas:
pass
f = open(MeteData)
data = f.readlines()
data2 =' '.join(data)
shutil.copyfile(MeteData,os.path.join(outname,os.path.basename(MeteData)))
if len(os.path.basename(MeteData))<10:
RSbands = glob.glob(os.path.join(root,"B0[1-8].tiff"))
else:
RSbands = glob.glob(os.path.join(root,"*B[1-8].TIF"))
print('影像'+root+'开始大气校正')
print(RSbands)
for tifFile in RSbands:
BandId = (os.path.basename(tifFile).split('.')[0])[-1]
#捕捉打开数据出错异常
try:
IDataSet = gdal.Open(tifFile)
except Exception as e:
print("文件%S打开失败" % tifFile)
LogFile.write('\n'+tifFile+'数据打开失败')
if IDataSet == None:
LogFile.write('\n'+tifFile+'数据集读取为空')
continue
else:
#获取行列号
cols = IDataSet.RasterXSize
rows = IDataSet.RasterYSize
ImgBand = IDataSet.GetRasterBand(1)
ImgRasterData = ImgBand.ReadAsArray(0, 0, cols, rows)
if ImgRasterData is None:
LogFile.write('\n'+tifFile+'栅格数据为空')
continue
else:
#设置输出文件路径
outFilename=os.path.join(outname,os.path.basename(tifFile))
#如果文件存在就跳过，进行下一波段操作
if os.path.isfile(outFilename):
print("%s已经完成" % outFilename)
continue
else:
# #辐射校正
RaCalRaster = RadiometricCalibration(BandId)
#大气校正
a, b, c = AtmosphericCorrection(BandId)
y = numpy.where(RaCalRaster!=-9999,a * RaCalRaster - b,-9999)
atc = numpy.where(y!=-9999,(y / (1 + y * c))*10000,-9999)
driver = IDataSet.GetDriver()
#输出栅格数据集
outDataset = driver.Create(outFilename, cols, rows, 1, gdal.GDT_Int16)
# 设置投影信息，与原数据一样
geoTransform = IDataSet.GetGeoTransform()
outDataset.SetGeoTransform(geoTransform)
proj = IDataSet.GetProjection()
outDataset.SetProjection(proj)
outband = outDataset.GetRasterBand(1)
outband.SetNoDataValue(-9999)
outband.WriteArray(atc, 0, 0)
print('第%s波段计算完成'%BandId)
# print(root+'计算完成')
print('\n')
#关闭日志文件
LogFile.close()
三、运行结果
1）python运行结果
2）校正后文件目录
3）结果显示
原始影像
该代码大气校正影像
4216
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
