用c++编写一个遥感图像大气校正代码

遥感图像大气校正是利用大气模型对遥感图像进行校正，以消除大气散射对图像的影响。以下是一个简单的C++代码示例，用于对Landsat 8 OLI图像进行大气校正：

#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

// 定义常量
const double pi = 3.14159265358979323846;
const double esun[] = { 1903.93, 1962.41, 1559.96, 1044.05, 225.71, 82.07 };
const double d2r = pi / 180.0;
const double r2d = 180.0 / pi;
const double L8_W[] = { 0.00002, 0.00002, 0.00002, 0.00002, 0.00002, 0.00002 };

// 声明函数
double deg2rad(double deg);
double rad2deg(double rad);
double get_theta(double sun_zenith, double view_zenith, double relative_azimuth);
double get_airmass(double theta, double altitude);
double get_transmittance(double airmass, double wavelength);
double get_reflectance(double TOA, double transmittance, double sun_zenith, double theta, double altitude);
double get_radiance(double reflectance, double transmittance, double sun_zenith, double theta, double altitude, double wavelength);

int main()
{
    // 读取输入文件
    ifstream infile("Landsat8_OLI.txt");
    double meta_data[10];
    for (int i = 0; i < 10; i++) {
        infile >> meta_data[i];
    }
    double height = meta_data[0];
    double sun_zenith = meta_data[1];
    double sun_azimuth = meta_data[2];
    double view_zenith = meta_data[3];
    double view_azimuth = meta_data[4];
    double relative_azimuth = meta_data[5];
    double TOA[] = { meta_data[6], meta_data[7], meta_data[8], meta_data[9] };

    // 计算大气校正系数
    double theta = get_theta(sun_zenith, view_zenith, relative_azimuth);
    double altitude = height / 1000.0;
    double airmass = get_airmass(theta, altitude);
    double transmittance[6];
    for (int i = 0; i < 6; i++) {
        transmittance[i] = get_transmittance(airmass, L8_W[i]);
    }
    double reflectance[4];
    for (int i = 0; i < 4; i++) {
        reflectance[i] = get_reflectance(TOA[i], transmittance[i], sun_zenith, theta, altitude);
    }
    double radiance[4];
    for (int i = 0; i < 4; i++) {
        radiance[i] = get_radiance(reflectance[i], transmittance[i], sun_zenith, theta, altitude, L8_W[i]);
    }

    // 输出结果
    cout << "大气校正系数：" << endl;
    cout << "波段1：" << reflectance[0] << endl;
    cout << "波段2：" << reflectance[1] << endl;
    cout << "波段3：" << reflectance[2] << endl;
    cout << "波段4：" << reflectance[3] << endl;
    cout << "辐亮度：" << endl;
    cout << "波段1：" << radiance[0] << endl;
    cout << "波段2：" << radiance[1] << endl;
    cout << "波段3：" << radiance[2] << endl;
    cout << "波段4：" << radiance[3] << endl;

    return 0;
}

// 角度与弧度转换
double deg2rad(double deg)
{
    return deg * d2r;
}
double rad2deg(double rad)
{
    return rad * r2d;
}

// 计算观测角
double get_theta(double sun_zenith, double view_zenith, double relative_azimuth)
{
    double cos_theta = cos(deg2rad(sun_zenith)) * cos(deg2rad(view_zenith))
        + sin(deg2rad(sun_zenith)) * sin(deg2rad(view_zenith)) * cos(deg2rad(relative_azimuth));
    return rad2deg(acos(cos_theta));
}

// 计算大气质量
double get_airmass(double theta, double altitude)
{
    double t1 = 1.0 / (cos(deg2rad(theta)) + 0.50572 * pow(96.07995 - altitude, -1.6364));
    double t2 = 1.0 / (cos(deg2rad(theta)) + 0.48596 * pow(96.07995 - altitude, -1.4494));
    return (t1 + t2) / 2.0;
}

// 计算大气透过率
double get_transmittance(double airmass, double wavelength)
{
    double taumol = exp(-0.008735 * pow(wavelength, -4) * airmass);
    double tauaer = exp(-0.003448 * airmass);
    return taumol * tauaer;
}

// 计算表观反射率
double get_reflectance(double TOA, double transmittance, double sun_zenith, double theta, double altitude)
{
    double cos_theta = cos(deg2rad(theta));
    double sin_theta = sin(deg2rad(theta));
    double cos_sun_zenith = cos(deg2rad(sun_zenith));
    double Rmol = 0.008569 * transmittance * pow(cos_theta, 4) * (1.0 + 0.0113 * exp(-1.5 * (cos_theta - 1.0)));
    double Raer = 0.002276 * transmittance * pow(cos_theta, 2) * (1.0 - 0.666 * exp(-0.00146 * pow(1.0 / cos_theta, 1.3)));
    double Rg = (TOA - Raer - Rmol) / (transmittance * cos_sun_zenith);
    double Tg = (pi * esun[0] * pow(cos_sun_zenith, 2)) / (L8_W[0] * pow(altitude + 1.0, 2));
    return Rg * Tg;
}

// 计算辐亮度
double get_radiance(double reflectance, double transmittance, double sun_zenith, double theta, double altitude, double wavelength)
{
    double cos_theta = cos(deg2rad(theta));
    double cos_sun_zenith = cos(deg2rad(sun_zenith));
    double Lmol = 0.008735 * transmittance * pow(cos_theta, 4) * esun[0] * (1.0 + 0.0113 * exp(-1.5 * (cos_theta - 1.0)));
    double Laer = 0.002276 * transmittance * pow(cos_theta, 2) * esun[0] * (1.0 - 0.666 * exp(-0.00146 * pow(1.0 / cos_theta, 1.3)));
    double Lg = reflectance * transmittance * cos_sun_zenith * (pi * esun[0] * pow(cos_sun_zenith, 2)) / (L8_W[0] * pow(altitude + 1.0, 2));
    return Lg + Laer + Lmol;
}

注意：以上代码仅供参考，实际应用中需要根据不同的遥感卫星、传感器和大气模型进行调整和修改。