本文只讲干货，不讲客套话，不讲没有由头的技术细节。目的只为将6S模型用于国产卫星大气校正。
6S模型介绍
首先来讲为什么要进行大气校正，讲讲必要性。首先，只有定量化需求才会用到大气校正，当然为了满足项目要求，或者满足老板、甲方要求，强制做大气校正也没啥不可。我还见过对GF系列全色波段进行大气校正呢，技术可行，合理不合理另说。
接下来讲讲为什么要用6S模型，而不是Flaash,或者modtran，因为前者常见于ENVI软件自带，集成有一定约束性。Modtran3以前笔者的印象是精度稍低（也可能我不会用），毕竟看名字嘛modtran（MODerate resolution atmospheric TRANsmission）更适用于中尺度，Modtran4+后来没再用。6S模型使用比较熟练了，所以后续的开发和集成都是基于6S了,甚至笔者用C++重写了一遍6S模型，可以定制各种相机传感器。
6S模型都需要输入哪些参数？具体6S背景介绍就不提了，可以自行搜索一下，这里直接介绍技术细节。6S 的输入参数主要有 9 个部分组成，如下表所示：
| 组成部分 | 内容 |
| 几何参数igeom | 6S 两种输入方法： ⑴ 太阳和卫星的天顶角和方位角以及观测时间(月,日)。 ⑵卫星的接收时间(月,日,年)、像素点数、升交点时间，由程序计算太阳和卫星的天顶角和方位角。特别注意的是这里的时间采用世界时且要精确到 1/6 秒。 |
| 大气模式idatm | 6S 给出几种可供选择的大气模式： Idatm =0：无气体吸收 idatm =1：热带大气 idatm =2：中纬度夏大气 idatm =3：中纬度冬季 idatm =4：亚北极区夏季 idatm =5：亚北极区冬季 idatm =6：美国标准大气（62年） idatm =7：用户定义大气廓线（34层无线电探空数据） 也可自定义大气模式。这个通过简单的时间和经纬度信息可以实现自动判断。 |
| 气溶胶模式iaer | 取值范围：0-12 iaer=0: 无气溶胶 iaer=1: 大陆型气溶胶 iaer=2: 海洋型气溶胶 iaer=3: 城市气溶胶 iaer=5: 沙漠型气溶胶 iaer=6: 生物质燃烧型 iaer=7: 平流层模式 iaer=4: 用户自己输入以下四种粒子所占体积百分比（0-1） c(1) ： 灰尘 c(2) ：水溶型 c(3) ：海洋型 c(4) ：烟灰 iaer=8-10：用户自己按照尺度分布类型定义气溶胶模型 iaer=8：多峰对数正态分布 iaer=9：改进的gamma分布 iaer=10：Junge幂指数律分布 iaer=11：按太阳光度计测量结果定义气溶胶模型 需要输入参数有：粒子半径（µm） 粒径分布（d V / d (logr)，cm3/cm2/micron）和复折射指数的实部和虚部谱。 |
| 气溶胶含量参数 | 两种选择：⑴ 在 550nm 处的光学厚度 ⑵ 气象能见度（km）。故它也提供了两者的相互关系。 取值范围： v=能见度（公里，能见度必须大于5公里） v=0：输入550纳米气溶胶光学厚度 taer55=550纳米气溶胶光学厚度 v=-1：没有气溶胶 |
| 目标高度 | 以千米为单位的地面海拔高度（设为负值）。这里要注意单位一定要换算到KM，并且注意正负号 |
| 传感器高度 | -1000 代表卫星测量，0 为地基观测，飞机航测输入以千米为单位的负值。 xpp= -1000：卫星观测 xpp= 0：地面观测 -100< xpp <0：飞机观测，绝对值代表飞机相对于目标的高度（公里） |
| 光谱参数iwave | 给出了常见卫星 Meteosat，Goes，NOAA/AVHRR 和 HRV，Landsat TM 和 MSS，Modis Polder 的每个通道的光谱响应函数，也可选择自定义。 iwave=-2 – +1， 用户自己定义光谱条件 iwave=-2：用户输入光谱范围的下限和上限（微米），滤光片函数为1，输出文件中给出单色结果。 iwave=-1：单色计算，用户给出单色波长（微米） iwave=0：用户输入光谱范围的下限和上限（微米），滤光片函数为1 iwave=1：用户输入光谱范围的下限和上限（微米）并输入滤光片函数（间隔为0.0025微米。（自定义输入波段范围和光谱响应率） 从iwave=2开始，是一些预设的传感器参数。笔者就是在这一部分增加的国产卫星载荷的标记。 |
| 地表特性inhomo | 可以选择地表均一或不均一，也可选择地表为郎伯体或双向反射。6S 给出了九种比较成熟的 BRDF 模式供用户选择，也可自定义BRDF函数（输入个角度的反射率及入射强度） inhomo=0：均匀表面 所需参数： idirec=0： 无方向效应 输入均匀朗伯表面的反射率igroun（roc=roe） idirec=1： 有方向效应 ibrdf=0： 输入太阳天顶角为thetas时10个观测天顶角（0-80度间隔10度和85度）和30个观测方位角（0-360度间隔30度）下的反射率； 同样，输入观测天顶角为thetav时各太阳入射角度下的反射率； 地表半球反射率； 在所选的观测条件下（太阳天顶角，观测天顶角和相对方位角）的反射率； ibrdf=1-9： 选择模式中储存的模式 ibrdf=1： hapke model ibrdf=2： verstraete et al. model ibrdf=3： Roujean et al. model ibrdf=4： walthall et al. model ibrdf=5： minnaert model ibrdf=6： Ocean ibrdf=7： Iaquinta and Pinty model ibrdf=8： Rahman et al. model ibrdf=9： Kuusk's multispectral CR model |
| 大气校正方式 | 输入反射率或辐射亮度，同时也决定模式是正向还是反向工作。当 RAPP<-1 时是正向。RAPP>0（辐射亮度）或-1<RAPP <0（反射率）均决定是反向过程， 即要进行大气订正过程。>rapp<-1: 不激活大气订正方式 rapp> 0：反演地面反射率，大气层顶的辐射亮度=rapp（w/m2/str/mic） -1.<rapp<0：反演地面反射率，表观反射率= rapp |
6S模型最初基于fortran语言编写实现。大家常用的方法是，直接运行命令行程序，将输入参数构造成6S所需的输入文件。然后6S模型作为一个黑盒子，根据给定的参数，计算出结果并写入到输出文本文件。笔者使用的版本是6Sv4.1版本。传统的6S输入是这样子的（输入举例）：
| 行数 | 内容 | 释意 |
|---|---|---|
| 1 | 0 | 几何参数（igeom）为用户自定义输入 |
| 2 | 88 180 40 0 3 18 | 太阳天顶角（度） 太阳方位角（度） 卫星天顶角（度） 卫星方位角（度） 月（1-12） 日（1-31） |
| 3 | 2 | 大气模式：2中纬度夏大气 |
| 4 | 1 | 气溶胶类型参数：1大陆型气溶胶 |
| 5 | 0 | 气溶胶含量参数：v=0：输入550纳米气溶胶光学厚度 |
| 6 | 1.950000 | 550纳米气溶胶光学厚度 |
| 7 | 0.000000 | 目标高度参数 |
| 8 | -1000 | 传感器高度参数：xpp= -1000：卫星观测 |
| 9 | 1 | 光谱参数，自定义输入波段范围和光谱响应率 |
| 10 | 0.630000 0.690000 | 波段范围 |
| 11 | 0.232000 0.323000 0.413000 0.506000 0.589000 0.682000 0.766000 0.851000 0.936000 0.953000 0.969000 0.985000 1.000000 0.988000 90976.000000 0.965000 0.954000 0.926000 0.898000 0.871000 0.843000 0.757000 0.671000 0.584000 0.497000 | 光谱响应率，可以理解为在波段范围内间隔步长（2.5um）x内的y值 |
| 12 | 0 | 地表反射率类型, inhomo=0：均匀表面 |
| 13 | 0 | idirec=0：不考虑方向效应 |
| 14 | 1 | igroun = 1; 绿色植被 |
| 15 | -2 | rapp = -2; % 无大气校正，正向运算 |
将上述参数保存成TXT作为输入传递给6S模型，然后计算出结果保存成TXT文件，我们需要从结果TXT文件中获取xa,xb,xc三个系数。上述思路的前提是将6s模型算法编译成命令行程序，换成动态库形式原理也类似。
注意，如果都能获取到上述表格中的参数，那么不论是国产卫星还是火星产卫星，都是支持的。但是国产卫星通常有个风格就是参数不便获取，或者即使能获取，位置也比较隐晦，需要有经验和关系才能获取到。获取到参数只是第一步，第二部如何将参数信息（比如传感器定标系数，波谱响应函数等）转化成6S的输入规范也是一个比较麻烦的工作，虽然技术上没有啥难度，但是用户的精力不在此，时间用在刀刃上，这种造轮子技术细节应该有人来做，所以我索性把卫星类型做成了6S模型的一个参数项，用户只要指定了卫星传感器类型，相应的定标系数和光谱响应函数就都准备好了。此外，设定了参数，6S模型也能知道该类型数据该如何从哪里（通常是图像数据的辅助XML文件）获取观测几何等其他元数据信息，极大简化了大气模型调用技术复杂度，节省了用户宝贵科研时间（摸鱼时间）。
通常国产卫星的幅宽都不是特别大，因此认为整个图像具备相同的成像几何条件也说得过去（也就是说一个波段只需要运行一次6S模型）。但是如果向精细计算，那么就要考虑分块（甚至按像素）获取成像几何参数（主要是卫星天顶角，卫星方位角，太阳天顶角，太阳方位角），意味着一张图像要多次运行6S模型，来计算不同成像几何下的6S计算结果3参数。
6S模型C++算法实现
笔者根据6Sv4.1版本，将fortran语言进行C++重构，形成了C++版本的6S模型库，为了照顾fortran版本的使用习惯，C++版本也沿用了TXT文件输入形式，并且格式完全一致。
参数化配置
用户在大气校正之前需要进行辐射定标，确保大气校正模块输入数据的数据具有物理含义。本文辐射定标功能的定标系数来自卫星官网。本文已经处理成JSON格式。中国的资源系列、高分系列、环境系列等系列卫星定标系数原始文件参见（好像已经失效20230411）：http://www.cresda.com/CN/Downloads/dbcs/index.shtml
为了便于用户扩展，定标系数制作成了配置文件形式，便于用户扩展，文件格式为json格式。以GF2多光谱数据为例，定标系数文件如下所示：
上述JSON文件中，数组长度等于图像波段数，日期年份代表地面定标的时期（注意，这是笔者自定义的格式，不是必须的！）。很多传感器在使用寿命内会不定期检测辐射标定情况，以保持精度。上述JSON内容的输入来自资源卫星应用中心官网，目前（20230411）已经失效，笔者将搜集到的原始文件直接放到了CSDN，供大家免费下载使用。
https://download.csdn.net/download/hanbing6174/87676199
大气校正采用6S模型，国产高分系列、资源系列等卫星的光谱响应函数一般不包含在图像中，像这类传感器，我们把光谱响应函数也放到了配置文件夹中，文件格式为txt格式。还是以GF2多光谱数据为例，里面只有4行，每行1501个数据。其中4行代表图像波段数，每行1501个代表从波长250nm开始，步长2.5nm，终止到4000nm波长的光谱响应。（（注意，这是笔者自定义的格式，不是必须的！，你完全可以自己造个形式，然后只要满足6S模型输入即可）
上述TXT内容的输入来自资源卫星应用中心官网，目前（20230411）已经失效，笔者将搜集到的原始文件直接放到了CSDN，供大家免费下载使用。
https://download.csdn.net/download/hanbing6174/87676199
如果用户需要同时做几何校正和辐射校正，用户应该首先进行辐射校正，然后在进行几何校正。因为辐射校正结果会保留图像的RPC等元数据文件。而几何校正不会保留。
本平台支持的辐射定标和大气校正支持的传感器类型列表如下表：
| 卫星型号 | 传感器型号 | 辐射定标 | 辐射定标输入 | 大气校正 | 大气校正输入 | 大气校正数据源 |
| 高分一号 | 全色相机（PMS1，PMS2） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 |
| 多光谱相机（PMS1，PMS2） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 | |
| B、C、D星全色 | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 | |
| B、C、D星多光谱 | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 | |
| 多光谱（WFV） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 | |
| 高分二号 | 全色相机（PMS1，PMS2） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 |
| 多光谱相机（PMS1，PMS2） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 | |
| 高分四号 | 全色多光谱相机 | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 |
| 中波红外相机 | 支持 | tiff图像 | 否 | - | - | |
| 高分五号 | 可见光近红外VN | 支持 | Xml文件 | 支持 | Xml文件 | 辐射定标结果 |
| 短波红外SW | 支持 | tiff图像 | 支持 | |||
| 高分六号 | 全色多光谱（PMS） | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 |
| 多光谱宽幅 | 支持 | tiff图像 | 支持 | Xml文件 | 辐射定标结果 | |
| 资源1号02C | PMS相机 | 支持 | tiff图像 | 支持 | tiff图像 | 辐射定标结果 |
| HRC相机（宽幅） | 否 | - | 否 | - | - | |
| 环境星 | HJ1A-CCD1 | 支持 | Xml文件 | 支持 | Xml文件 | 辐射定标结果 |
| HJ1A-CCD2 | 支持 | Xml文件 | 支持 | Xml文件 | 辐射定标结果 | |
| HJ1B-CCD1 | 支持 | Xml文件 | 支持 | Xml文件 | 辐射定标结果 | |
| HJ1B-CCD2 | 支持 | Xml文件 | 支持 | Xml文件 | 辐射定标结果 | |
| 资源三号 | ZY3-01 多光谱相机 | 支持 | tiff图像 | 否 | tiff图像 | 辐射定标结果 |
| ZY3-02 多光谱相机 | 支持 | tiff图像 | 否 | tiff图像 | 辐射定标结果 | |
| Landsat5 | TM | 支持 | Txt文件(*_MTL.txt) | 支持 | Txt文件(*_MTL.txt) | 原始数据 |
| Landsat7 | ETM | 否 | - | 否 | - | - |
| Landsat8 | L1C | 支持 | Txt文件(*_MTL.txt) | 支持 | Txt文件(*_MTL.txt) | 原始数据 |
| 哨兵2A | MSI | 否 | - | 否 | - | - |
| 哨兵2B | MSI | 否 | - | 否 | - | - |
其中，高分四号中波红外相机因为波长(3300nm-4350nm)超出了6S处理范围(250nm-4000nm)因此不被支持。Landsat7因为传感器故障原因，暂时不被支持。更多国产卫星介绍参见：
6S校正流程
添加新的传感器类型，主要涉及三大部分，分别是MetadataHandle ，RadiationCalcHandle，AtmosCorr6SHandle。这是三个回调函数。分别实现参数读取，辐射定标和大气校正。
RadiationCalcHandle实现
/**
* @brief 辐射定标函数(自动识别传感器类型)。
* @param data_path 定标前的输入数据路径。 注意，根据卫星类型不同，输入的影像形式有所不同。通常直接输入TIFF结尾的影像全路径，但是以下类型
* 需要输入图像的元数据文件：
* LANDSAT5_TM：输入*_MTL.txt文件
* SEN_LANDSAT8：输入*_MTL.txt文件
* @param save_dir 定标后的辐射亮度影像输出文件夹。
* @param progress 进度条。
*/
LIB_WIMATMOS Error AfxRadiationCalibration(
std::string &dnfile,
std::string &save_dir,
Progress *progress = NULL)
{
CalibSensorType st = GetCalibSensorTypeByFile(dnfile);
if (st == SEN_Undefine) {
RETURN_ERROR(_tr("unrecognized satellite image : ") << dnfile);
}
//获取辐射校正处理方法
auto radiation_handle = RadiationCalcEntry(st);
if (!radiation_handle)
return ErrNone;
//获取元数据提取方法
auto meta_handle = MetadataEntry(st);
if (!meta_handle)
return ErrNone;
AtmosMatedata mdata(st);
meta_handle(dnfile, mdata);
wim::CreateAllPaths(save_dir);
radiation_handle(dnfile, save_dir, mdata.gain_offs, progress);
return ErrNone;
}
其中，MetadataEntry来自以下函数：
MetadataHandle MetadataEntry(CalibSensorType st)
{
switch (st)
{
case SEN_Undefine :nullptr;
case SEN_GF1_WFV1 :return Metadata_GF1_WFV1 ;
case SEN_GF1_WFV2 :return Metadata_GF1_WFV2 ;
case SEN_GF1_WFV3 :return Metadata_GF1_WFV3 ;
case SEN_GF1_WFV4 :return Metadata_GF1_WFV4 ;
case SEN_GF1_PAN1 :return Metadata_GF1_PAN1 ;
case SEN_GF1_MSS1 :return Metadata_GF1_MSS1 ;
case SEN_GF1_PAN2 :return Metadata_GF1_PAN2 ;
case SEN_GF1_MSS2 :return Metadata_GF1_MSS2 ;
case SEN_GF1B_PAN :return Metadata_GF1B_PAN ;
case SEN_GF1B_MUX :return Metadata_GF1B_MUX ;
case SEN_GF1C_PAN :return Metadata_GF1C_PAN ;
case SEN_GF1C_MUX :return Metadata_GF1C_MUX ;
case SEN_GF1D_PAN :return Metadata_GF1D_PAN ;
case SEN_GF1D_MUX :return Metadata_GF1D_MUX ;
case SEN_GF2_PAN1 :return Metadata_GF2_PAN1 ;
case SEN_GF2_MSS1 :return Metadata_GF2_MSS1 ;
case SEN_GF2_PAN2 :return Metadata_GF2_PAN2 ;
case SEN_GF2_MSS2 :return Metadata_GF2_MSS2 ;
case SEN_GF4_PMS :return Metadata_GF4_PMS ;
case SEN_GF4_IRS :return Metadata_GF4_IRS;
case SEN_GF5_VN :return Metadata_GF5_VN;
case SEN_GF5_SW :return Metadata_GF5_SW;
case SEN_GF6_PAN :return Metadata_GF6_PAN;
case SEN_GF6_MUX :return Metadata_GF6_MUX;
case SEN_GF6_WFV :return Metadata_GF6_WFV;
case SEN_HJ1A_CCD1 :return Metadata_HJ1A_CCD1 ;
case SEN_HJ1A_CCD2 :return Metadata_HJ1A_CCD2 ;
case SEN_HJ1B_CCD1 :return Metadata_HJ1B_CCD1 ;
case SEN_HJ1B_CCD2 :return Metadata_HJ1B_CCD2 ;
case SEN_ZY02C_PAN :return Metadata_ZY02C_PAN;
case SEN_ZY02C_MUX :return Metadata_ZY02C_MSS;
case SEN_ZY02C_HRC :return Metadata_ZY02C_HRC;
case SEN_ZY301_MUX :return Metadata_ZY301_MSS ;
case SEN_ZY302_MUX :return Metadata_ZY302_MSS ;
case SEN_LANDSAT5_TM :return Metadata_LANDSAT5_TM ;
case SEN_LANDSAT7_TM :return Metadata_LANDSAT7_TM ;
case SEN_LANDSAT8 :return Metadata_LANDSAT8 ;
default: break;
}
return nullptr;
}
以Metadata_GF1_WFV1为例，来看下实现：
Error Metadata_GF1_WFV1(const std::string& data_path, AtmosMatedata& mdata)
{
std::string fxml = wim::Path(data_path).ReplaceExtension(".xml");
if (!wim::Exists(fxml)) {
RETURN_ERROR(_tr("xml file not exist!") << fxml);
}
TiXmlDocument xml_doc;
xml_doc.LoadFile(fxml.c_str());
if (xml_doc.Error())
return Error(xml_doc.ErrorDesc());
TiXmlNode* root = xml_doc.RootElement();
if (NULL == root)
return Error(xml_doc.ErrorDesc());
//基本信息.波段数bands，高分系列并不统一，比如GF1B的xml中bands指的是波段数，二GF2中确是波段序号。因此基本信息干脆直接从图像里面读取\n
RasterDatasetPtr ds = OpenRasterDataset(data_path);
if (!ds){
RETURN_ERROR("failed to open image: " << data_path);
}
mdata.bands = ds->GetBandCount();
mdata.lines = ds->GetYSize();
mdata.samples = ds->GetXSize();
//时间信息
get_gf_xml_time(root, mdata);
//坐标信息
mdata.srs = ds->GetSpatialRef();
if (!mdata.srs){//不存在空间参考信息，采用XML信息
get_gf_xml_envelope(root, mdata);
mdata.srs = SpatialRef::CreateFromWkt(GCS_WGS_1984);
}
else {
mdata.envelope = ds->GetEnvelope();
}
//观测几何参数
get_gf_xml_geom(root, mdata);
//大气模式
mdata.idatm = CalcIdatmByTimeGeo(mdata.month, mdata.day, (mdata.envelope.getMaxY() + mdata.envelope.getMinY()));
//气溶胶模式（iaer）,必须用户填写。一般选择大陆型气溶胶。
//气溶胶数值（visibility或者aod550二选一）,必须用户填写。
//高程数值（height）必须用户填写。后续可考虑改用全球DEM。
//获取定标系数文件
GetGainOffset(mdata.gain_offs, mdata.sensor, mdata.year);
return ErrNone;
}
上述函数主要工作就是构造AtmosMatedata这个类。这个类包含了辐射定标和大气校正所需要的所有参数。定义如下：
/**
* @brief 记录每个波段的波谱信息。
*/
struct WaveInfo
{
double wavelen; //波长(单位纳米)
double fwhm; //半高宽(单位纳米)
};
typedef boost::shared_ptr<WaveInfo> WaveInfoPtr;
typedef std::vector<WaveInfoPtr> WaveInfoList;
typedef wim::StringList BandNameList;
/**
* @brief 参与大气校正、定量计算的图像所需要具备的信息。
*/
class LIB_WIMATMOS AtmosMatedata
{
public:
AtmosMatedata(CalibSensorType st);
//检测输入参数是否完备
Error CheckParamOK();
public:
//传感器参数
CalibSensorType sensor; //传感器类型
//基本信息
int bands;
int lines;
int samples;
//时间信息
Time start_time; //起始观测时间
Time end_time; //截止观测时间
int year;
int month; //观测月
int day; //观测日
//坐标信息
wim::Envelope envelope;
SpatialRefPtr srs;
//观测几何参数
bool use_vertical; //是否垂直观测
double solar_zenith; //太阳天顶角
double solar_azimuth; //太阳方位角
double satellite_zenith; //卫星天顶角
double satellite_azimuth; //卫星方位角
/** 大气模式
* idatm =0：无气体吸收
* idatm = 1：热带大气
* idatm = 2：中纬度夏大气
* idatm = 3：中纬度冬季
* idatm = 4：亚北极区夏季
* idatm = 5：亚北极区冬季
*/
//bool use_idatm; //是否使用大气模式
int idatm; //大气模式
/** 气溶胶模式
* iaer=0: 无气溶胶
* iaer=1: 大陆型气溶胶
* iaer=2: 海洋型气溶胶
* iaer=3: 城市气溶胶
* iaer=5: 沙漠型气溶胶
* iaer=6: 生物质燃烧型
* iaer=7: 平流层模式
*/
int iaer;
//气溶胶参数
bool use_visibility;
double visibility; //能见度，KM
double aod550; //550纳米气溶胶厚度
//高程
double height; //地面海拔,单位km
//定标系数文件
std::vector<CalibFactor> gain_offs;
};
之所以每个传感器类型都对应一个元数据提取函数，是因为每个传感器的元数据格式，读取方式都存在差别。当然高分系列卫星一致性比较好，因此对于读取方式相通的，调用相同的读取函数即可。
RadiationCalcEntry来自以下函数：
RadiationCalcHandle RadiationCalcEntry(CalibSensorType st)
{
switch (st)
{
case SEN_Undefine :nullptr;
case SEN_GF1_WFV1 :return RadiationCalc_GF1_WFV1 ;
case SEN_GF1_WFV2 :return RadiationCalc_GF1_WFV2 ;
case SEN_GF1_WFV3 :return RadiationCalc_GF1_WFV3 ;
case SEN_GF1_WFV4 :return RadiationCalc_GF1_WFV4 ;
case SEN_GF1_PAN1 :return RadiationCalc_GF1_PAN1 ;
case SEN_GF1_MSS1 :return RadiationCalc_GF1_MSS1 ;
case SEN_GF1_PAN2 :return RadiationCalc_GF1_PAN2 ;
case SEN_GF1_MSS2 :return RadiationCalc_GF1_MSS2 ;
case SEN_GF1B_PAN :return RadiationCalc_GF1B_PAN ;
case SEN_GF1B_MUX :return RadiationCalc_GF1B_MUX ;
case SEN_GF1C_PAN :return RadiationCalc_GF1C_PAN ;
case SEN_GF1C_MUX :return RadiationCalc_GF1C_MUX ;
case SEN_GF1D_PAN :return RadiationCalc_GF1D_PAN ;
case SEN_GF1D_MUX :return RadiationCalc_GF1D_MUX ;
case SEN_GF2_PAN1 :return RadiationCalc_GF2_PAN1 ;
case SEN_GF2_MSS1 :return RadiationCalc_GF2_MSS1 ;
case SEN_GF2_PAN2 :return RadiationCalc_GF2_PAN2 ;
case SEN_GF2_MSS2 :return RadiationCalc_GF2_MSS2 ;
case SEN_GF4_PMS :return RadiationCalc_GF4_PMS ;
case SEN_GF4_IRS :return RadiationCalc_GF4_IRS ;
case SEN_GF5_VN :return RadiationCalc_GF5_VN ;
case SEN_GF5_SW :return RadiationCalc_GF5_SW ;
case SEN_GF6_PAN :return RadiationCalc_GF6_PAN ;
case SEN_GF6_MUX :return RadiationCalc_GF6_MUX ;
case SEN_GF6_WFV :return RadiationCalc_GF6_WFV ;
case SEN_HJ1A_CCD1 :return RadiationCalc_HJ1A_CCD1 ;
case SEN_HJ1A_CCD2 :return RadiationCalc_HJ1A_CCD2 ;
case SEN_HJ1B_CCD1 :return RadiationCalc_HJ1B_CCD1 ;
case SEN_HJ1B_CCD2 :return RadiationCalc_HJ1B_CCD2 ;
case SEN_ZY02C_PAN :return RadiationCalc_ZY02C_PAN ;
case SEN_ZY02C_MUX :return RadiationCalc_ZY02C_MSS ;
case SEN_ZY301_MUX :return RadiationCalc_ZY301_MUX ;
case SEN_ZY302_MUX :return RadiationCalc_ZY302_MUX ;
case SEN_LANDSAT5_TM :return RadiationCalc_LANDSAT5_TM;
case SEN_LANDSAT7_TM :return RadiationCalc_LANDSAT7_TM;
case SEN_LANDSAT8 :return RadiationCalc_LANDSAT8 ;
default: break;
}
return nullptr;
}
来看下RadiationCalc_GF1_WFV1的实现
Error RadiationCalc_GF1_WFV1 (const std::string& data_path, const std::string& save_dir, std::vector<CalibFactor>& gain_offs, Progress *progress)
{
std::string base_name = wim::Path(data_path).FileName();
std::string califile = wim::Path(save_dir).Append(base_name);
//对于这类定标系数是单独公布的(不包含在图像文件中)，除了计算图像本身，需要拷贝RPB文件、xml文件、jpg缩略图。
CalcRadiationCalibration(data_path, califile, gain_offs, progress);
//拷贝图像
MakeFileCopy(wim::Path(data_path).ReplaceExtension(".rpb"), wim::Path(califile).ReplaceExtension(".rpb"));
MakeFileCopy(wim::Path(data_path).ReplaceExtension(".xml"), wim::Path(califile).ReplaceExtension(".xml"));
MakeFileCopy(wim::Path(data_path).ReplaceExtension(".jpg"), wim::Path(califile).ReplaceExtension(".jpg"));
std::string name = wim::Path(data_path).Stem();
MakeFileCopy(wim::Path(data_path).ReplaceFileName(name + "_thumb").ReplaceExtension(".jpg"), wim::Path(califile).ReplaceFileName(name + "_thumb").ReplaceExtension(".jpg"));
return ErrNone;
}
AtmosCorr6SHandle实现
/**
* @brief 6s大气校正
* @param src_img 输入辐射亮度文件。注意，根据卫星类型不同，输入的影像形式有所不同。通常直接输入TIFF结尾的影像全路径，但是以下类型
* 需要输入图像的元数据文件：
* LANDSAT5_TM：输入*_MTL.txt文件
* SEN_LANDSAT8：输入*_MTL.txt文件
* @param save_dir 输出大气校正文件所在的文件夹。
* @param param 输入参数
* @param zoom10k 是否放大10000倍结果，以便保存为无符号短整型（u short），否则保存为浮点型（float）
* @param progress 进度条
* @return 错误对象
* @note 不支持 SEN_LANDSAT7_TM
*/
LIB_WIMATMOS Error AtmosCorrection6S(
std::string src_img,
std::string save_dir,
const AtmosMatedata ¶m,
bool zoom10k = true,
Progress* progress = NULL)
{
CalibSensorType st = param.sensor;
Error errt = ErrNone;
if (st == SEN_Undefine) {
RETURN_ERROR(_tr("unrecognized satellite image : ") << src_img);
}
//获取辐射校正处理方法
auto atmoscorr_handle = AtmosCorr6SEntry(st);
if (!atmoscorr_handle)
return ErrNone;
return atmoscorr_handle(src_img, save_dir, param, zoom10k, progress);
}
其中，AtmosCorr6SHandle来自以下函数
AtmosCorr6SHandle AtmosCorr6SEntry(CalibSensorType st)
{
switch (st)
{
case SEN_Undefine :nullptr;
case SEN_GF1_WFV1 :return AtmosCorr6S_GF1_WFV1;
case SEN_GF1_WFV2 :return AtmosCorr6S_GF1_WFV2;
case SEN_GF1_WFV3 :return AtmosCorr6S_GF1_WFV3;
case SEN_GF1_WFV4 :return AtmosCorr6S_GF1_WFV4;
case SEN_GF1_PAN1 :return AtmosCorr6S_GF1_PAN1;
case SEN_GF1_MSS1 :return AtmosCorr6S_GF1_MSS1;
case SEN_GF1_PAN2 :return AtmosCorr6S_GF1_PAN2;
case SEN_GF1_MSS2 :return AtmosCorr6S_GF1_MSS2;
case SEN_GF1B_PAN :return AtmosCorr6S_GF1B_PAN;
case SEN_GF1B_MUX :return AtmosCorr6S_GF1B_MUX;
case SEN_GF1C_PAN :return AtmosCorr6S_GF1C_PAN;
case SEN_GF1C_MUX :return AtmosCorr6S_GF1C_MUX;
case SEN_GF1D_PAN :return AtmosCorr6S_GF1D_PAN;
case SEN_GF1D_MUX :return AtmosCorr6S_GF1D_MUX;
case SEN_GF2_PAN1 :return AtmosCorr6S_GF2_PAN1;
case SEN_GF2_MSS1 :return AtmosCorr6S_GF2_MSS1;
case SEN_GF2_PAN2 :return AtmosCorr6S_GF2_PAN2;
case SEN_GF2_MSS2 :return AtmosCorr6S_GF2_MSS2;
case SEN_GF4_PMS :return AtmosCorr6S_GF4_PMS ;
case SEN_GF4_IRS :return AtmosCorr6S_GF4_IRS ;
case SEN_GF5_VN :return AtmosCorr6S_GF5_VN ;
case SEN_GF5_SW :return AtmosCorr6S_GF5_SW ;
case SEN_GF6_PAN :return AtmosCorr6S_GF6_PAN ;
case SEN_GF6_MUX :return AtmosCorr6S_GF6_MUX ;
case SEN_GF6_WFV :return AtmosCorr6S_GF6_WFV ;
case SEN_HJ1A_CCD1 :return AtmosCorr6S_HJ1A_CCD1;
case SEN_HJ1A_CCD2 :return AtmosCorr6S_HJ1A_CCD2;
case SEN_HJ1B_CCD1 :return AtmosCorr6S_HJ1B_CCD1;
case SEN_HJ1B_CCD2 :return AtmosCorr6S_HJ1B_CCD2;
case SEN_ZY02C_MUX :return AtmosCorr6S_ZY02C_MUX;
case SEN_ZY301_MUX :return AtmosCorr6S_ZY301_MUX;
case SEN_ZY302_MUX :return AtmosCorr6S_ZY302_MUX;
case SEN_LANDSAT5_TM :return AtmosCorr6S_LANDSAT5_TM;
case SEN_LANDSAT7_TM :return AtmosCorr6S_LANDSAT7_TM;
case SEN_LANDSAT8 :return AtmosCorr6S_LANDSAT8;
default: break;
}
return nullptr;
}
以AtmosCorr6S_GF1_WFV1函数为例，介绍详细实现
Error AtmosCorr6S_GF1_WFV1 (const std::string& src_img, const std::string& save_dir, const AtmosMatedata ¶m6s, bool zoom10k, Progress *progress)
{
Progress::Wrapper prog(progress);
prog.Create("6s Correction...");
prog.SetTotalLoops(2);
std::string base_name = wim::Path(src_img).FileName();
std::string dst_img = wim::Path(save_dir).Append(base_name);
//计算大气校正系数
std::vector< SixS3Param > param3;
Error er = Calc6SParam(param3, src_img, dst_img, param6s, progress);
if (er != ErrNone) {
return er;
}
//进行大气校正
Error errt = AtmosphericCorrection_Radiation(src_img, dst_img, param3, zoom10k, progress);
if (errt != ErrNone)
return errt;
//拷贝图像
MakeFileCopy(wim::Path(src_img).ReplaceExtension(".rpb"), wim::Path(dst_img).ReplaceExtension(".rpb"));
MakeFileCopy(wim::Path(src_img).ReplaceExtension(".xml"), wim::Path(dst_img).ReplaceExtension(".xml"));
MakeFileCopy(wim::Path(src_img).ReplaceExtension(".jpg"), wim::Path(dst_img).ReplaceExtension(".jpg"));
std::string name = wim::Path(src_img).Stem();
MakeFileCopy(wim::Path(src_img).ReplaceFileName(name + "_thumb").ReplaceExtension(".jpg"), wim::Path(dst_img).ReplaceFileName(name + "_thumb").ReplaceExtension(".jpg"));
return ErrNone;
}
其中，比较关键的函数是Calc6SParam该函数主要是将用户输入的AtmosMatedata，按照本文一开始介绍的形式，构建TXT，然后通过6S算法返回三参数，用SixS3Param去记录反演结果。还是以GF1_WFV1为例，第一个波段，构建的txt名称为GF1_WFV1_E117.9_N24.6_20200415_L1A0004742543_b0_6sparam.txt，内容如下：
这里面的内容参见第一章节6S介绍部分。其中比较重要的参数是81这个参数并不是通用的，你用你的6S模型处理GF1_WFV1数据时，填写81是得不到结果的，因为这个81是本文自定义的传感器类型。运行完6S模型，会生成GF1_WFV1_E117.9_N24.6_20200415_L1A0004742543_6slog.txt，里面记录了输出结果（也就是三个系数），这里的格式和fortran语言的6S有些差别。但是原理一样：
拿到这系数，根据提示公式，逐像素遍历计算即可得到大气校正后的反射率。
上文说到，81这个参数并不通用，这里是笔者自定义的，定义方法如下：
/* ********************************************************************c
c 2 vis band of meteosat ( 0.350-1.110 ) c
c 3 vis band of goes east ( 0.490-0.900 ) c
c 4 vis band of goes west ( 0.490-0.900 ) c
c 5 1st band of avhrr(noaa6) ( 0.550-0.750 ) c
c 6 2nd " ( 0.690-1.120 ) c
c 7 1st band of avhrr(noaa7) ( 0.500-0.800 ) c
c 8 2nd " ( 0.640-1.170 ) c
c 9 1st band of avhrr(noaa8) ( 0.540-1.010 ) c
c 10 2nd " ( 0.680-1.120 ) c
c 11 1st band of avhrr(noaa9) ( 0.530-0.810 ) c
c 12 2nd " ( 0.680-1.170 ) c
c 13 1st band of avhrr(noaa10 ( 0.530-0.780 ) c
c 14 2nd " ( 0.600-1.190 ) c
c 15 1st band of avhrr(noaa11 ( 0.540-0.820 ) c
c 16 2nd " ( 0.600-1.120 ) c
c 17 1st band of hrv1(spot1) ( 0.470-0.650 ) c
c 18 2nd " ( 0.600-0.720 ) c
c 19 3rd " ( 0.730-0.930 ) c
c 20 pan " ( 0.470-0.790 ) c
c 21 1st band of hrv2(spot1) ( 0.470-0.650 ) c
c 22 2nd " ( 0.590-0.730 ) c
c 23 3rd " ( 0.740-0.940 ) c
c 24 pan " ( 0.470-0.790 ) c
c 25 1st band of tm(landsat5) ( 0.430-0.560 ) c
c 26 2nd " ( 0.500-0.650 ) c
c 27 3rd " ( 0.580-0.740 ) c
c 28 4th " ( 0.730-0.950 ) c
c 29 5th " ( 1.5025-1.890 ) c
c 30 7th " ( 1.950-2.410 ) c
c 31 1st band of mss(landsat5)( 0.475-0.640 ) c
c 32 2nd " ( 0.580-0.750 ) c
c 33 3rd " ( 0.655-0.855 ) c
c 34 4th " ( 0.785-1.100 ) c
c 35 1st band of MAS (ER2) ( 0.5025-0.5875) c
c 36 2nd " ( 0.6075-0.7000) c
c 37 3rd " ( 0.8300-0.9125) c
c 38 4th " ( 0.9000-0.9975) c
c 39 5th " ( 1.8200-1.9575) c
c 40 6th " ( 2.0950-2.1925) c
c 41 7th " ( 3.5800-3.8700) c
c 42 MODIS band 1 ( 0.6100-0.6850) c
c 43 MODIS band 2 ( 0.8200-0.9025) c
c 44 MODIS band 3 ( 0.4500-0.4825) c
c 45 MODIS band 4 ( 0.5400-0.5700) c
c 46 MODIS band 5 ( 1.2150-1.2700) c
c 47 MODIS band 6 ( 1.6000-1.6650) c
c 48 MODIS band 7 ( 2.0575-2.1825) c
c 49 1st band of avhrr(noaa12 ( 0.500-1.000 ) c
c 50 2nd " ( 0.650-1.120 ) c
c 51 1st band of avhrr(noaa14 ( 0.500-1.110 ) c
c 52 2nd " ( 0.680-1.100 ) c
c 53 POLDER band 1 ( 0.4125-0.4775) c
c 54 POLDER band 2 (non polar( 0.4100-0.5225) c
c 55 POLDER band 3 (non polar( 0.5325-0.5950) c
c 56 POLDER band 4 P1 ( 0.6300-0.7025) c
c 57 POLDER band 5 (non polar( 0.7450-0.7800) c
c 58 POLDER band 6 (non polar( 0.7000-0.8300) c
c 59 POLDER band 7 P1 ( 0.8100-0.9200) c
c 60 POLDER band 8 (non polar( 0.8650-0.9400) c
c 61 1st band of etm+(landsat7( 0.435-0.520 ) c
c 62 2nd " ( 0.506-0.621 ) c
c 63 3rd " ( 0.622-0.702 ) c
c 64 4th " ( 0.751-0.911 ) c
c 65 5th " ( 1.512-1.792 ) c
c 66 7th " ( 2.020-2.380 ) c
c 67 8th " ( 0.504-0.909 ) c
c 68 2nd band of liss (IRC 1C)( 0.502-0.620 ) c
c 69 3rd " ( 0.612-0.700 ) c
c 70 4th " ( 0.752-0.880 ) c
c 71 5th " ( 1.452-1.760 ) c
c 72 1st band of aster ( 0.480-0.645 ) c
c 73 2nd " ( 0.588-0.733 ) c
c 74 3N " ( 0.723-0.913 ) c
c 75 4th " ( 1.530-1.750 ) c
c 76 5th " ( 2.103-2.285 ) c
c 77 6th " ( 2.105-2.298 ) c
c 78 7th " ( 2.200-2.393 ) c
c 79 8th " ( 2.248-2.475 ) c
c 80 9th " ( 2.295-2.538 ) c
c 81 1st band of GF1_WFV1 " ( ) c
c 82 2nd " ( ) c
c 83 3 " ( ) c
c 84 4 " ( ) c
c 85 1st band of GF1_WFV2 ( ) c
c 86 2 " ( ) c
c 87 3 " ( ) c
c 88 4 " ( ) c
c 89 1st band of GF1_WFV3 ( ) c
c 90 2 " ( ) c
c 91 3 " ( ) c
c 92 4 " ( ) c
c 93 1st band of GF1_WFV4 ( ) c
c 94 2 " ( ) c
c 95 3 " ( ) c
c 96 4 " ( ) c
c 97 pan band of GF1_PAN1 ( ) c
c 98 1st band of GF1_MSS1 ( ) c
c 99 2 " ( ) c
c 100 3 ( ) c
c 101 4 " ( ) c
c 102 pan band of GF1_PAN2 ( ) c
c 103 1st band of GF1_MSS2 ( ) c
c 104 2 " ( ) c
c 105 3 ( ) c
c 106 4 " ( ) c
c 107 pan band of GF1B_PAN ( ) c
c 108 1st band of GF1B_MUX ( ) c
c 109 2 ( ) c
c 110 3 ( ) c
c 111 4 ( ) c
c 112 pan band of GF1C_PAN ( ) c
c 113 1st band of GF1C_MUX ( ) c
c 114 2 ( ) c
c 115 3 ( ) c
c 116 4 ( ) c
c 117 pan band of GF1D_PAN ( ) c
c 118 1st band of GF1D_MUX ( ) c
c 119 2 ( ) c
c 120 3 ( ) c
c 121 4 ( ) c
c 122 pan band of GF2_PAN1 ( ) c
c 123 1st band of GF2_MSS1 ( ) c
c 124 2 ( ) c
c 125 3 ( ) c
c 126 4 ( ) c
c 127 pan band of GF2_PAN2 ( ) c
c 128 1st band of GF2_MSS2 ( ) c
c 129 2 ( ) c
c 130 3 ( ) c
c 131 4 ( ) c
c 132 1st band of GF4_PMS ( ) c
c 133 2 ( ) c
c 134 3 ( ) c
c 135 4 ( ) c
c 136 5 ( ) c
c 137 1st band of GF5_VN ( ) c
c 286 150 band of GF5_VN ( ) c
c 287 1st band of GF5_SW ( ) c
c 466 180 band of GF5_SW ( ) c
c 467 pan band of GF6_PAN ( ) c
c 468 1st band of GF6_MUX ( ) c
c 469 2 ( ) c
c 470 3 ( ) c
c 471 4 ( ) c
c 472 1st band of GF6_WFV ( ) c
c 473 2 ( ) c
c 474 3 ( ) c
c 475 4 ( ) c
c 475 5 ( ) c
c 477 6 ( ) c
c 478 7 ( ) c
c 479 8 ( ) c
c 480 1st band of HJ1A_CCD1 ( ) c
c 481 2 ( ) c
c 482 3 ( ) c
c 483 4 ( ) c
c 484 1st band of HJ1A_CCD2 ( ) c
c 485 2 ( ) c
c 486 3 ( ) c
c 487 4 ( ) c
c 488 1st band of HJ1B_CCD1 ( ) c
c 489 2 ( ) c
c 490 3 ( ) c
c 491 4 ( ) c
c 492 1st band of HJ1B_CCD2 ( ) c
c 493 2 ( ) c
c 494 3 ( ) c
c 495 4 ( ) c
c 496 1st band of ZY02C_MUX ( ) c
c 497 2 ( ) c
c 498 3 ( ) c
c 499 1st band of ZY301_MUX ( ) c
c 500 2 ( ) c
c 501 3 ( ) c
c 502 4 ( ) c
c note: wl has to be in micrometer c
c**********************************************************************/
也就是从81开始，笔者自己定义了传感器处理代码索引，接下来介绍如何处理自定义的代码索引：
void IWave::parse()
{
iinf = 0;
isup = 1500;
int i;
for(i = 0; i <= isup; i++)
ffu.s[i] = 1;
cin >> iwave;
cin.ignore(numeric_limits<int>::max(),'\n');
if(iwave == 0 || iwave == -2) {
cin >> ffu.wlinf;
cin >> ffu.wlsup;
cin.ignore(numeric_limits<int>::max(),'\n');
} else if(iwave < 0) { /* excludes -2 */
cin >> wl;
cin.ignore(numeric_limits<int>::max(),'\n');
ffu.wlinf = wl;
ffu.wlsup = wl;
} else {
/* da big switch */
if (iwave == 1) {
cin >> ffu.wlinf;
cin >> ffu.wlsup;
cin.ignore(numeric_limits<int>::max(), '\n');
/* moved rest further on */
}
else if (iwave == 2)
meteo();
else if (iwave == 3)
goes_east();
else if (iwave == 4)
goes_west();
else if (iwave <= 16)
avhrr(iwave - 4);
else if (iwave <= 24)
hrv(iwave - 16);
else if (iwave <= 30)
tm(iwave - 24);
else if (iwave <= 34)
mss(iwave - 30);
else if (iwave <= 41)
mas(iwave - 34);
else if (iwave <= 48)
modis(iwave - 41);
else if (iwave <= 52)
avhrr(iwave - 36);
else if (iwave <= 60)
polder(iwave - 52);
else if (iwave <= 67)
etmplus(iwave - 60);
else if (iwave <= 71)
irs_1c_liss(iwave - 67);
else if (iwave <= 80)
aster(iwave - 71);
else if (iwave <= 84)//81 1st band of GF1_WFV1
GF1_WFV1(iwave - 80);
else if (iwave <= 88)//85 1st band of GF1_WFV2
GF1_WFV2(iwave - 84);
else if (iwave <= 92)//89 1st band of GF1_WFV3
GF1_WFV3(iwave - 88);
else if (iwave <= 96)//93 1st band of GF1_WFV4
GF1_WFV4(iwave - 92);
else if (iwave <= 97)//97 pan band of GF1_PAN1
GF1_PAN1(iwave - 96);
else if (iwave <= 101)//98 1st band of GF1_MSS1
GF1_MSS1(iwave - 97);
else if (iwave <= 102)//102 pan band of GF1_PAN2
GF1_PAN2(iwave - 101);
else if (iwave <= 106)//103 1st band of GF1_MSS2
GF1_MSS2(iwave - 102);
else if (iwave <= 107)//107 pan band of GF1B_PAN
GF1B_PAN(iwave - 106);
else if (iwave <= 111)//108 pan band of GF1B_MUX
GF1B_MUX(iwave - 107);
else if (iwave <= 112)//112 1st band of GF1C_PAN
GF1C_PAN(iwave - 111);
else if (iwave <= 116)//113 1st band of GF1C_MUX
GF1C_MUX(iwave - 112);
else if (iwave <= 117)//117 1st band of GF1D_PAN
GF1D_PAN(iwave - 116);
else if (iwave <= 121)//118 1st band of GF1D_MUX
GF1D_MUX(iwave - 117);
else if (iwave <= 122)//122 1st band of GF2_PAN1
GF2_PAN1(iwave - 121);
else if (iwave <= 126)//123 1st band of GF2_MSS1
GF2_MSS1(iwave - 122);
else if (iwave <= 127)//127 1st band of GF2_PAN2
GF2_PAN2(iwave - 126);
else if (iwave <= 131)//128 1st band of GF2_MSS2
GF2_MSS2(iwave - 127);
else if (iwave <= 136)//1st band of GF4_PMS
GF4_PMS(iwave - 131);
else if (iwave <= 128)//125 1st band of GF4_MMS
printf("Unsupported iwave value: %d", iwave);
else if (iwave <= 286)//1st band of GF5_VN
GF5_VN(iwave - 136);
else if (iwave <= 466)//1st band of GF5_SW
GF5_SW(iwave - 286);
else if (iwave <= 467)//129 1st band of GF6_PAN
GF6_PAN(iwave - 466);
else if (iwave <= 471)//129 1st band of GF6_MSS
GF6_MSS(iwave - 467);
else if (iwave <= 479)//133 1st band of GF6_WFV
GF6_WFV(iwave - 471);
else if (iwave <= 483)//141 1st band of HJ1A_CCD1
HJ1A_CCD1(iwave - 479);
else if (iwave <= 487)//145 1st band of HJ1A_CCD2
HJ1A_CCD2(iwave - 483);
else if (iwave <= 491)//149 1st band of HJ1B_CCD1
HJ1B_CCD1(iwave - 487);
else if (iwave <= 495)//153 1st band of HJ1B_CCD1
HJ1B_CCD2(iwave - 491);
else if (iwave <= 498)//157 1st band of ZY1_02C_MSS
ZY02C_MUX(iwave - 495);
else if (iwave <= 502)//160 1st band of ZY3_MSS
ZY301_MUX(iwave - 498);
else
printf("Unsupported iwave value: %d",iwave);
}
iinf = (int)((ffu.wlinf - 0.25f) / 0.0025f + 1.5f) - 1; /* remember indexing*/
isup = (int)((ffu.wlsup - 0.25f) / 0.0025f + 1.5f) - 1; /* "*/
if(iwave == 1) { /* moved here to avoid unnecessery gotos */
for(int i = iinf; i <= isup; i++)
cin >> ffu.s[i];
cin.ignore(numeric_limits<int>::max(),'\n');
}
}
还是以GF1_WFV1为例，介绍实现细节：
void IWave::GF1_WFV1(int iwa)
{
static const float wli[4] = { 0.45, 0.52, 0.63, 0.77 };
static const float wls[4] = { 0.52, 0.59, 0.69, 0.89 };
ffu.wlinf = wli[iwa - 1];
ffu.wlsup = wls[iwa - 1];
for (int i = 0; i < 1501; i++)
ffu.s[i] = 0.0;
GetSpectrumFileLine(SEN_GF1_WFV1, iwa, ffu.s);
}
其中，GetSpectrumFileLine作用是根据传感器类型，找到配置参数文件GF1_WFV1.txt也就是光谱响应函数GF1_WFV1.txt的内容介绍见上文“参数化配置”。
全文完。有问题欢迎留言交流讨论，共同学习进步。也可加我V，hanbing6174，注明来意。
1961
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
