import numpy as np
from osgeo import gdal
def calculate_reflectance(panel_dn_list, panel_reflectance, input_path, output_path, default_bands=(105, 59, 34),
sensor_params={
"FriendlyName": "Basler acA1920-155um (24997755)",
"FrameRate": 128.766418,
"hwFrameRate": 25,
"Integral": 3000,
"Gain": 0.0
}, wavelength_list=None):
# 参数说明
# panel_dn_list: 灰板各波段DN值列表(需按波段顺序)
# panel_reflectance: 灰板反射率(0.5)
# input_path: 输入ENVI格式.hdr文件路径
# output_path: 输出反射率文件路径
# 打开输入文件
ds = gdal.Open(input_path, gdal.GA_ReadOnly)
if ds is None:
raise ValueError("无法打开输入文件!")
# 获取图像元数据
bands = ds.RasterCount
cols = ds.RasterXSize
rows = ds.RasterYSize
# 校验波段数一致性
if len(panel_dn_list) != bands:
print("图像波段数:", bands)
print("灰板DN值数量:", len(panel_dn_list))
raise ValueError("灰板DN值数量与图像波段数不匹配!")
# 创建输出文件
driver = gdal.GetDriverByName('ENVI')
out_ds = driver.Create(output_path, cols, rows, bands, gdal.GDT_Float32)
# 逐波段处理
for b in range(bands):
# 读取当前波段数据
band_data = ds.GetRasterBand(b + 1).ReadAsArray().astype(np.float32)
# 计算增益系数(处理除零错误)
panel_dn = panel_dn_list[b]
if panel_dn <= 0:
gain = 0.0
else:
gain = panel_reflectance / panel_dn
# 计算反射率
reflectance = band_data * gain
# 写入输出
out_band = out_ds.GetRasterBand(b + 1)
out_band.WriteArray(reflectance)
out_band.SetDescription(f"Band {b + 1} Reflectance")
# 生成ENVI头文件
wavelength_str = ', '.join(map(str, wavelength_list)) if wavelength_list else ''
header_content = f"""ENVI
description = {{
File Imported into ENVI.}}
samples = {cols}
lines = {rows}
bands = {bands}
header offset = 0
file type = ENVI Standard
data type = 4 # Float32对应类型码4
interleave = bil
sensor type = Unknown
byte order = 0
x start = 0
y start = 0
default bands = {{{','.join(map(str, default_bands))}}}
FriendlyName:{sensor_params['FriendlyName']}
FrameRate = {sensor_params['FrameRate']}
hwFrameRate = {sensor_params['hwFrameRate']}
Integral = {sensor_params['Integral']}
Gain = {sensor_params['Gain']:.6f}
wavelength units = Nanometers
wavelength = {{
{wavelength_str}
}}"""
with open(output_path + '.hdr', 'w') as f:
f.write(header_content)
# 释放资源
ds = None
out_ds = None
if __name__ == "__main__":
# 使用示例
panel_dn = {75.7568, 72.6486, 75.5135, 77.8649, 78.8108, 78.8649, 82.1081, 84.4865, 89.1081, 98.1081, 107.9459,
124.6486,
146.9189, 172.2162, 208.2973, 241.7568, 279.7838, 321.2432, 370.7027, 410.027, 455.4324, 496.2703,
490.7838,
616.1892, 712.3514, 832.8649, 863.4324, 996.1892, 1127.6486, 1197.2703, 1332.5676, 1479.5676, 1569.6757,
1622.4054,
1759.2703, 1829.5135, 1948.2703, 2063.7838, 2091.5405, 2156.8649, 2248.2973, 2360.4865, 2475.9189,
2531.3514,
2665.3784, 2737.8919, 2814.027, 2725.3514, 2609.9189, 2910.9459, 3024.0541, 3098.027, 3241.6216,
3268.4595,
3257.4865, 3246.8919, 3358.3243, 3506.9459, 3509.4054, 3600.5405, 3584.0541, 3556.6216, 3323.8919,
3553.3243,
3751.027, 3849.3784, 3818.8108, 3830.5946, 4030.3784, 3988.8919, 4013.5676, 4045.4054, 4005.6486,
3923.3243,
4010.973, 4052.1081, 4033.7838, 4048.2703, 4046.7838, 4051.0811, 4030.1081, 3934.973, 3927.1892,
3950.9459,
3916.8919, 3872.1892, 3843.2162, 3766.2162, 3852.7027, 3768.8649, 3742.2703, 3721.3784, 3779.3243,
3775.5405,
3675.9459, 3519.6486, 3414.7297, 3485.0, 3442.3243, 3412.3514, 3376.7838, 3305.3784, 3314.8378,
3290.1351,
3237.7568, 3164.3243, 3128.3784, 3060.3784, 2993.3243, 2999.8378, 2993.1351, 2929.8108, 2843.8378,
2816.7297,
2755.3784, 2732.8919, 2725.7297, 2717.6486, 2713.6757, 2665.3514, 2621.8378, 2595.9459, 2542.1081,
2492.4324,
2477.0541, 2483.8919, 2377.6216, 2220.8919, 2379.2703, 2402.3243, 2375.2703, 2361.5676, 2335.4054,
2310.6216,
2266.1622, 2250.2432, 2216.8919, 2187.5135, 2173.9459, 2142.9459, 2081.2973, 1805.4865, 1756.3243,
1848.5676,
1893.7838, 1888.2973, 1905.5946, 1859.6486, 1821.7838, 1833.7297, 1834.5946, 1825.8378, 1805.0541,
1774.4324,
1735.3514, 1598.9189, 1422.1081, 1488.1351, 1505.9459, 1445.973, 1447.5946, 1436.5135, 1467.8108,
1505.0811,
1485.8649, 1469.2162, 1436.2703, 1432.4324, 1426.4054, 1418.2703, 1384.2432, 1357.027, 1335.5135,
1317.5946,
1279.7568, 862.1892, 556.4865, 676.1892, 899.7297, 1084.5946, 1129.2432, 1120.9189, 1100.9459,
1086.6757,
1063.4865, 1041.2162, 1021.2432, 1002.6216, 978.6486, 949.4324, 919.5135, 894.0541, 881.6486, 870.6757,
849.9189,
828.7027, 807.8649, 798.6486, 773.0541, 746.5946, 704.6216, 636.9459, 598.973, 595.1351, 605.4324,
554.7297,
595.1351, 571.1351, 549.2703, 547.7297, 542.5405, 545.9459, 539.1892, 529.7568, 515.0, 503.0811,
490.027, 471.1351,
448.973, 430.6757, 408.3784, 431.2162, 417.9189, 408.0541, 402.1622, 365.2973, 364.5135, 365.5405,
356.1892,
343.3514, 336.3514, 326.7568, 313.3243, 308.8378, 298.9459, 284.1622, 278.9459, 272.0811, 256.7297,
232.973,
213.6486, 199.6216, 197.7568, 212.2162, 198.973, 176.973, 172.9459, 171.6486, 163.8919, 159.8649, 160.0,
159.0,
155.1351, 146.7568, 141.5676, 119.2432, 98.4054, 82.6216, 80.6757, 83.5135, 97.1622, 87.6486, 80.6216,
78.1892,
81.4324, 79.6216, 78.7027, 77.0, 80.9459, 78.9189, 79.0541, 81.2703, 76.3784, 81.3243, 83.8919, 78.4054,
74.3243,
72.8378, 74.1622, 72.2162, 70.4324, 66.7297, 68.8919, 68.0, 65.3514, 66.0, 63.1622, 63.8378, 59.9459,
60.4595,
58.7568, 54.7838, 54.7568, 54.7297, 55.9459, 51.3243, 52.0811, 48.8649, 49.0, 48.1351}
input_hdr = "D://高光谱数据九三农场//高光谱//九三5航高60//20250414134504.spe"
output_file = "D://高光谱数据九三农场//reflectance"
calculate_reflectance(
panel_dn_list=panel_dn,
panel_reflectance=0.5,
input_path=input_hdr,
output_path=output_file,
wavelength_list=[375.6236, 377.86405, 380.10332, 382.3414, 384.57833, 386.81409, 389.04872, 391.28221,
393.51458, 395.74585, 397.97601, 400.20508, 402.43307, 404.73992, 406.9486, 409.15618,
411.36264, 413.56801, 415.77231, 417.97553, 419.96943, 422.13808, 424.35444, 426.58128,
428.85044, 431.12472,
433.40251, 435.67257, 437.93984, 440.18425, 442.37761, 444.55551, 446.72028, 448.88473,
451.03863, 453.19598, 455.35828, 457.52154, 459.67714, 461.82435, 463.97206, 466.14151,
468.26702, 470.413, 472.54014, 474.69333, 476.85601, 479.04709, 481.29279, 483.54239,
485.78723,
488.03179, 490.27648, 492.47912, 494.68291, 496.8622, 499.02272, 501.18392, 503.32623,
505.46991, 507.6186, 509.76711, 511.90234, 514.04112, 516.17058, 518.30139, 520.42076,
522.52896, 524.65049, 526.77098, 528.92009, 531.08737, 533.29257, 535.50312, 537.74579,
539.98786,
542.22061, 544.4337, 546.62996, 548.80112, 550.95963, 553.10646, 555.24833, 557.3898,
559.51965, 561.66327, 563.8039, 565.9184, 568.01801, 570.13229, 572.23666, 574.32743,
576.43528, 578.55963, 580.69099, 582.83569, 585.01512, 587.20918, 589.43483, 591.6698,
593.89598,
596.09915, 598.28529, 600.4457, 602.59109, 604.70005, 606.82995, 608.95335, 611.08223,
613.20129, 615.32344, 617.46209, 619.58238, 621.69155, 623.79777, 625.90151, 628.00567,
630.11592, 632.24984, 634.40306, 636.56602, 638.73111, 640.9332, 643.1498, 645.3721, 647.57279,
649.74416, 651.90161, 654.04517, 656.16354, 658.26701, 660.37511, 662.4852, 664.58416,
666.68141, 668.77766, 670.86687, 672.96076, 675.04138, 677.12867, 679.21427, 681.32407,
683.45333, 685.5965, 687.76943, 689.98218, 692.17936, 694.40878, 696.60662, 698.8439,
701.01672,
703.16179, 705.3098, 707.44132, 709.57707, 711.71348, 713.83126, 715.95493, 718.0755,
720.18548, 722.28876, 724.38782, 726.48104, 728.57768, 730.67546, 732.7666, 734.90332,
737.05871, 739.24784, 741.47798, 743.7034, 745.90836, 748.11667, 750.29548, 752.43513,
754.5789,
756.6867, 758.79223, 760.91059, 763.02172, 765.13218, 767.23627, 769.33747, 771.4313,
773.51155, 775.60704, 777.69687, 779.79199, 781.89675, 784.02384, 786.17737, 788.37067,
790.58057, 792.80676, 795.02314, 797.23193, 799.40111, 801.55596, 803.68111, 805.7943,
807.87454,
809.98513, 812.09726, 814.21372, 816.31819, 818.40966, 820.50485, 822.58657, 824.65779,
826.74563, 828.84876, 830.9539, 833.09098, 835.25159, 837.45193, 839.66527, 841.90645,
844.1094, 846.29963, 848.44414, 850.58494, 852.71508, 854.82339, 856.94649, 859.06082,
861.19367,
863.27509, 865.38199, 867.46129, 869.54041, 871.60973, 873.68963, 875.79473, 877.91663,
880.05366, 882.21485, 884.4061, 886.61912, 888.85185, 891.06904, 893.24914, 895.39386,
897.52424, 899.65282, 901.75885, 903.86346, 905.97361, 908.08025, 910.19921, 912.29565,
914.37154,
916.44548, 918.51177, 920.5919, 922.68277, 924.78512, 926.91117, 929.08136, 931.29005,
933.50503, 935.70445, 937.90025, 940.07842, 942.20405, 944.33368, 946.44308, 948.52033,
950.61563, 952.71577, 954.8231, 956.93438, 958.98088, 961.04423, 963.11908, 965.17531,
967.26554,
969.34805, 971.46161, 973.62139, 975.81091, 978.03392, 980.23021, 982.4186, 984.58094,
986.71892, 988.838, 990.94424, 993.05166, 995.12993, 997.11134, 999.22406, 1001.33633,
1003.44813, 1005.55947, 1007.67033, 1009.78071, 1011.8906, 1013.99999, 1016.10889, 1018.21728]
)import numpy as np
from osgeo import gdal
def calculate_reflectance(panel_dn_list, panel_reflectance, input_path, output_path, default_bands=(105, 59, 34),
sensor_params={
"FriendlyName": "Basler acA1920-155um (24997755)",
"FrameRate": 128.766418,
"hwFrameRate": 25,
"Integral": 3000,
"Gain": 0.0
}, wavelength_list=None):
# 参数说明
# panel_dn_list: 灰板各波段DN值列表(需按波段顺序)
# panel_reflectance: 灰板反射率(0.5)
# input_path: 输入ENVI格式.hdr文件路径
# output_path: 输出反射率文件路径
# 打开输入文件
ds = gdal.Open(input_path, gdal.GA_ReadOnly)
if ds is None:
raise ValueError("无法打开输入文件!")
# 获取图像元数据
bands = ds.RasterCount
cols = ds.RasterXSize
rows = ds.RasterYSize
# 校验波段数一致性
if len(panel_dn_list) != bands:
print("图像波段数:", bands)
print("灰板DN值数量:", len(panel_dn_list))
raise ValueError("灰板DN值数量与图像波段数不匹配!")
# 创建输出文件
driver = gdal.GetDriverByName('ENVI')
out_ds = driver.Create(output_path, cols, rows, bands, gdal.GDT_Float32)
# 逐波段处理
for b in range(bands):
# 读取当前波段数据
band_data = ds.GetRasterBand(b + 1).ReadAsArray().astype(np.float32)
# 计算增益系数(处理除零错误)
panel_dn = panel_dn_list[b]
if panel_dn <= 0:
gain = 0.0
else:
gain = panel_reflectance / panel_dn
# 计算反射率
reflectance = band_data * gain
# 写入输出
out_band = out_ds.GetRasterBand(b + 1)
out_band.WriteArray(reflectance)
out_band.SetDescription(f"Band {b + 1} Reflectance")
# 生成ENVI头文件
wavelength_str = ', '.join(map(str, wavelength_list)) if wavelength_list else ''
header_content = f"""ENVI
description = {{
File Imported into ENVI.}}
samples = {cols}
lines = {rows}
bands = {bands}
header offset = 0
file type = ENVI Standard
data type = 4 # Float32对应类型码4
interleave = bil
sensor type = Unknown
byte order = 0
x start = 0
y start = 0
default bands = {{{','.join(map(str, default_bands))}}}
FriendlyName:{sensor_params['FriendlyName']}
FrameRate = {sensor_params['FrameRate']}
hwFrameRate = {sensor_params['hwFrameRate']}
Integral = {sensor_params['Integral']}
Gain = {sensor_params['Gain']:.6f}
wavelength units = Nanometers
wavelength = {{
{wavelength_str}
}}"""
with open(output_path + '.hdr', 'w') as f:
f.write(header_content)
# 释放资源
ds = None
out_ds = None
if __name__ == "__main__":
# 使用示例
panel_dn = {75.7568, 72.6486, 75.5135, 77.8649, 78.8108, 78.8649, 82.1081, 84.4865, 89.1081, 98.1081, 107.9459,
124.6486,
146.9189, 172.2162, 208.2973, 241.7568, 279.7838, 321.2432, 370.7027, 410.027, 455.4324, 496.2703,
490.7838,
616.1892, 712.3514, 832.8649, 863.4324, 996.1892, 1127.6486, 1197.2703, 1332.5676, 1479.5676, 1569.6757,
1622.4054,
1759.2703, 1829.5135, 1948.2703, 2063.7838, 2091.5405, 2156.8649, 2248.2973, 2360.4865, 2475.9189,
2531.3514,
2665.3784, 2737.8919, 2814.027, 2725.3514, 2609.9189, 2910.9459, 3024.0541, 3098.027, 3241.6216,
3268.4595,
3257.4865, 3246.8919, 3358.3243, 3506.9459, 3509.4054, 3600.5405, 3584.0541, 3556.6216, 3323.8919,
3553.3243,
3751.027, 3849.3784, 3818.8108, 3830.5946, 4030.3784, 3988.8919, 4013.5676, 4045.4054, 4005.6486,
3923.3243,
4010.973, 4052.1081, 4033.7838, 4048.2703, 4046.7838, 4051.0811, 4030.1081, 3934.973, 3927.1892,
3950.9459,
3916.8919, 3872.1892, 3843.2162, 3766.2162, 3852.7027, 3768.8649, 3742.2703, 3721.3784, 3779.3243,
3775.5405,
3675.9459, 3519.6486, 3414.7297, 3485.0, 3442.3243, 3412.3514, 3376.7838, 3305.3784, 3314.8378,
3290.1351,
3237.7568, 3164.3243, 3128.3784, 3060.3784, 2993.3243, 2999.8378, 2993.1351, 2929.8108, 2843.8378,
2816.7297,
2755.3784, 2732.8919, 2725.7297, 2717.6486, 2713.6757, 2665.3514, 2621.8378, 2595.9459, 2542.1081,
2492.4324,
2477.0541, 2483.8919, 2377.6216, 2220.8919, 2379.2703, 2402.3243, 2375.2703, 2361.5676, 2335.4054,
2310.6216,
2266.1622, 2250.2432, 2216.8919, 2187.5135, 2173.9459, 2142.9459, 2081.2973, 1805.4865, 1756.3243,
1848.5676,
1893.7838, 1888.2973, 1905.5946, 1859.6486, 1821.7838, 1833.7297, 1834.5946, 1825.8378, 1805.0541,
1774.4324,
1735.3514, 1598.9189, 1422.1081, 1488.1351, 1505.9459, 1445.973, 1447.5946, 1436.5135, 1467.8108,
1505.0811,
1485.8649, 1469.2162, 1436.2703, 1432.4324, 1426.4054, 1418.2703, 1384.2432, 1357.027, 1335.5135,
1317.5946,
1279.7568, 862.1892, 556.4865, 676.1892, 899.7297, 1084.5946, 1129.2432, 1120.9189, 1100.9459,
1086.6757,
1063.4865, 1041.2162, 1021.2432, 1002.6216, 978.6486, 949.4324, 919.5135, 894.0541, 881.6486, 870.6757,
849.9189,
828.7027, 807.8649, 798.6486, 773.0541, 746.5946, 704.6216, 636.9459, 598.973, 595.1351, 605.4324,
554.7297,
595.1351, 571.1351, 549.2703, 547.7297, 542.5405, 545.9459, 539.1892, 529.7568, 515.0, 503.0811,
490.027, 471.1351,
448.973, 430.6757, 408.3784, 431.2162, 417.9189, 408.0541, 402.1622, 365.2973, 364.5135, 365.5405,
356.1892,
343.3514, 336.3514, 326.7568, 313.3243, 308.8378, 298.9459, 284.1622, 278.9459, 272.0811, 256.7297,
232.973,
213.6486, 199.6216, 197.7568, 212.2162, 198.973, 176.973, 172.9459, 171.6486, 163.8919, 159.8649, 160.0,
159.0,
155.1351, 146.7568, 141.5676, 119.2432, 98.4054, 82.6216, 80.6757, 83.5135, 97.1622, 87.6486, 80.6216,
78.1892,
81.4324, 79.6216, 78.7027, 77.0, 80.9459, 78.9189, 79.0541, 81.2703, 76.3784, 81.3243, 83.8919, 78.4054,
74.3243,
72.8378, 74.1622, 72.2162, 70.4324, 66.7297, 68.8919, 68.0, 65.3514, 66.0, 63.1622, 63.8378, 59.9459,
60.4595,
58.7568, 54.7838, 54.7568, 54.7297, 55.9459, 51.3243, 52.0811, 48.8649, 49.0, 48.1351}
input_hdr = "D://高光谱数据九三农场//高光谱//九三5航高60//20250414134504.spe"
output_file = "D://高光谱数据九三农场//reflectance"
calculate_reflectance(
panel_dn_list=panel_dn,
panel_reflectance=0.5,
input_path=input_hdr,
output_path=output_file,
wavelength_list=[375.6236, 377.86405, 380.10332, 382.3414, 384.57833, 386.81409, 389.04872, 391.28221,
393.51458, 395.74585, 397.97601, 400.20508, 402.43307, 404.73992, 406.9486, 409.15618,
411.36264, 413.56801, 415.77231, 417.97553, 419.96943, 422.13808, 424.35444, 426.58128,
428.85044, 431.12472,
433.40251, 435.67257, 437.93984, 440.18425, 442.37761, 444.55551, 446.72028, 448.88473,
451.03863, 453.19598, 455.35828, 457.52154, 459.67714, 461.82435, 463.97206, 466.14151,
468.26702, 470.413, 472.54014, 474.69333, 476.85601, 479.04709, 481.29279, 483.54239,
485.78723,
488.03179, 490.27648, 492.47912, 494.68291, 496.8622, 499.02272, 501.18392, 503.32623,
505.46991, 507.6186, 509.76711, 511.90234, 514.04112, 516.17058, 518.30139, 520.42076,
522.52896, 524.65049, 526.77098, 528.92009, 531.08737, 533.29257, 535.50312, 537.74579,
539.98786,
542.22061, 544.4337, 546.62996, 548.80112, 550.95963, 553.10646, 555.24833, 557.3898,
559.51965, 561.66327, 563.8039, 565.9184, 568.01801, 570.13229, 572.23666, 574.32743,
576.43528, 578.55963, 580.69099, 582.83569, 585.01512, 587.20918, 589.43483, 591.6698,
593.89598,
596.09915, 598.28529, 600.4457, 602.59109, 604.70005, 606.82995, 608.95335, 611.08223,
613.20129, 615.32344, 617.46209, 619.58238, 621.69155, 623.79777, 625.90151, 628.00567,
630.11592, 632.24984, 634.40306, 636.56602, 638.73111, 640.9332, 643.1498, 645.3721, 647.57279,
649.74416, 651.90161, 654.04517, 656.16354, 658.26701, 660.37511, 662.4852, 664.58416,
666.68141, 668.77766, 670.86687, 672.96076, 675.04138, 677.12867, 679.21427, 681.32407,
683.45333, 685.5965, 687.76943, 689.98218, 692.17936, 694.40878, 696.60662, 698.8439,
701.01672,
703.16179, 705.3098, 707.44132, 709.57707, 711.71348, 713.83126, 715.95493, 718.0755,
720.18548, 722.28876, 724.38782, 726.48104, 728.57768, 730.67546, 732.7666, 734.90332,
737.05871, 739.24784, 741.47798, 743.7034, 745.90836, 748.11667, 750.29548, 752.43513,
754.5789,
756.6867, 758.79223, 760.91059, 763.02172, 765.13218, 767.23627, 769.33747, 771.4313,
773.51155, 775.60704, 777.69687, 779.79199, 781.89675, 784.02384, 786.17737, 788.37067,
790.58057, 792.80676, 795.02314, 797.23193, 799.40111, 801.55596, 803.68111, 805.7943,
807.87454,
809.98513, 812.09726, 814.21372, 816.31819, 818.40966, 820.50485, 822.58657, 824.65779,
826.74563, 828.84876, 830.9539, 833.09098, 835.25159, 837.45193, 839.66527, 841.90645,
844.1094, 846.29963, 848.44414, 850.58494, 852.71508, 854.82339, 856.94649, 859.06082,
861.19367,
863.27509, 865.38199, 867.46129, 869.54041, 871.60973, 873.68963, 875.79473, 877.91663,
880.05366, 882.21485, 884.4061, 886.61912, 888.85185, 891.06904, 893.24914, 895.39386,
897.52424, 899.65282, 901.75885, 903.86346, 905.97361, 908.08025, 910.19921, 912.29565,
914.37154,
916.44548, 918.51177, 920.5919, 922.68277, 924.78512, 926.91117, 929.08136, 931.29005,
933.50503, 935.70445, 937.90025, 940.07842, 942.20405, 944.33368, 946.44308, 948.52033,
950.61563, 952.71577, 954.8231, 956.93438, 958.98088, 961.04423, 963.11908, 965.17531,
967.26554,
969.34805, 971.46161, 973.62139, 975.81091, 978.03392, 980.23021, 982.4186, 984.58094,
986.71892, 988.838, 990.94424, 993.05166, 995.12993, 997.11134, 999.22406, 1001.33633,
1003.44813, 1005.55947, 1007.67033, 1009.78071, 1011.8906, 1013.99999, 1016.10889, 1018.21728]
)
以上代码出现错误:line 33, in calculate_reflectance
raise ValueError("灰板DN值数量与图像波段数不匹配!")
最新发布