注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

面朝大海 春暖花开

一点心得,转载本博客文章,请注明转帖,谢谢!

 
 
 

日志

 
 
关于我

中国科学院博士 主要从事遥感机理、定量反演、数据处理以及GIS应用研究。ArcGIS、Envi 、ERDAS、Ecognition软件、IDL语言、6S、SAIL

网易考拉推荐

6s LUT[转载]  

2013-09-17 14:41:08|  分类: 遥感 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

逐像元大气校正,常预先计算查找表(LUT,LookUp Tabel),6S大气辐射传输模式也可以用来计算LUT。但6S源程序输出信息多,且浮点数输出精度低,不利于提取关键信息生成LUT,本文描述了怎样修改6S源码以生成LUT。

首先确定LUT内容要素。

???????? 阅读MODIS M?D04 气溶胶产品生成算法文档(NASA相关网页),归纳了以下查找表要素:
1)?????? 太阳天顶角 观测天顶角 太阳方位角观测方位角之差(相对方位角) 散射角
2)?????? 大气模式
3)?????? 气溶胶模式
4)?????? 550nm气溶胶光学厚度
5)?????? 波段号
6)?????? 大气透过率
7)?????? atm. intrin. ref.(个人理解,这是大气程辐射换算后的反射率,用于校正计算)
8)?????? total? sca. (总散射,包括rayleigh散射和气溶胶散射,用于校正计算)
9)?????? 表观反射率
10)??? 校正后的地表反射率
11)??? xa xb xc参数(物理意义不明,用于校正计算)

其次,阅读源码,明确LUT各要素在6S源码中的变量名。

6S大气校正计算源码(Excel验证中采用此公式)
???????? rog=rapp/tgasm
???????? rog=(rog-ainr(1,1)/tgasm)/sutott/sdtott
???????? rog=rog/(1.+rog*sast)
? xa=pi*sb/xmus/seb/tgasm/sutott/sdtott
? xb=srotot/sutott/sdtott/tgasm
? xc=sast
由计算源码确定需输出的变量。下面是输出LUT要素的Fortran77代码示例,建议放在源码main.f中stop一句之前。
?? write(*,*) asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55,
1?????? iwave,tgasm,ainr(1,1),sutott*sdtott,rapp,rog,xa,xb,xc
其中,asol是太阳天顶角,phi0是太阳方位角,avis是观测天顶角,phiv是观测方位角,adif是散射角,phi是相对方位角,idatm是大气模式号,iaer是气溶胶模式号,v是水平能见度,taer55是550nm气溶胶光学厚度,iwave表示波段号,tgasm表示大气透过率,ainr(1,1)是大气本身的反射率(姑且这么理解),sutott*sdtott表示总散射,rapp是表观反射率,rog是校正后的地表反射率,xa,xb,xc是校正计算参数。Fortran77每行长度不能超过80个字符,续前行需在特定位置指明(示例中的iwave前的1即表示续行)。
该示例源码没有指定任何输出样式。浮点数会按默认样式输出,小数点后的位数比较多,精度较好。挑选一个查找表文件用Excel验证后表明,excel公式计算值与6S输出值之间最大误差小于1E-7。说明方法是可行的。
再次,编译源码,编写Shell脚本:

编译环境:OpenSUSE操作系统 g95编译器,版本未知。
编译命令:g95 *.f -o sixs(在BRDF相关代码处可能有几个warning,本文不涉及BRDF,故暂不修改调试。在Windows下用f77编译,无warning,编译通过)
生成LUT的bash脚本getLUT.sh:
function LUTCalc(){
#{42,44,48,25,27,30}
IBand=1 echo"Band 1? echo "Band {IBand} is running"
for Iref in {0.1,0.5}
do
fstat=IBand" ′′ {IBand}"_" {Iref}.res???????
echo "asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc" >> {fstat}
for SolarZ in {0,15,30,45,75}
do
? for SolarAz in {0,45,90,135}
? do
?? for ViewZ in {0,15,30,45,75}
?? do
??? for Iaer in {1,2,3,5,6,7}
??? do
???? for Idatm in {1,2,3,4,5,6}
???? do
????? for IAod in {0.1,0.2,0.5,1.0};
????? do
#?? Run sixs and output
./sixs <<EOF | tail -n 1 >> {fstat}? for SolarZ in {0,15,30,45,75}? do?? for SolarAz in {0,45,90,135}?? do??? for ViewZ in {0,15,30,45,75}??? do???? for Iaer in {1,2,3,5,6,7}???? do????? for Idatm in {1,2,3,4,5,6}????? do?????? for IAod in {0.1,0.2,0.5,1.0};?????? do #?? Run sixs and output ./sixs <<EOF | tail -n 1 >> {fstat}
0
SolarZ SolarZ SolarAz ViewZ0315 ViewZ 0 3 15 Idatm
Iaer0 Iaer 0 IAod
-0
-1000
IBand0000.0? {IBand} 0 0 0 0.0 - Iref
EOF
?????? done
????? done
???? done
??? done
?? done
? done

done

}

function getLUT()
{
echo "LUT is Calcing"
for iwave in {42,44,48,25,27,30}
{
? LUTCalc ${iwave}
} &
}
最好晚上计算早晨看结果,如果CPU给力的话,几个小时后就可以得到结果。
下面是生成的LUT示例:
asol,phi0,avis,phiv,adif,phi,idatm,iaer,v,taer55, iwave,tgasm,ainr,tott,rapp,rog,xa,xb,xc
0.0000000?????? 0.0000000?????? 0.0000000?????? 0.0000000?????? 180.00000?????? 0.0000000??? 1?? 1?? 63.664398????? 0.10000000??? 25? 0.98984975????? 6.89081103E-02? 0.80637234????? 0.10000000????? 3.87301818E-02? 1.98730151E-03? 8.71319696E-02? 0.14777245???
?? 0.0000000?????? 0.0000000?????? 0.0000000?????? 0.0000000?????? 180.00000?????? 0.0000000??? 1?? 1?? 26.739149????? 0.20000000??? 25? 0.98984975????? 7.75793791E-02? 0.76532620????? 0.10000000????? 2.94530466E-02? 2.09388486E-03? 0.10336593????? 0.16389242???
?? 0.0000000?????? 0.0000000?????? 0.0000000?????? 0.0000000?????? 180.00000?????? 0.0000000??? 1?? 1?? 8.4940033????? 0.50000000??? 25? 0.98984975????? 0.10173188????? 0.64870018????? 0.10000000???? -2.69861287E-03? 2.47033220E-03? 0.15994170????? 0.20088956???
?? 0.0000000?????? 0.0000000?????? 0.0000000?????? 0.0000000?????? 180.00000?????? 0.0000000??? 1?? 1?? 3.5674956?????? 1.0000000??? 25? 0.98984975????? 0.13688390????? 0.48083964????? 0.10000000???? -7.89646432E-02? 3.33272247E-03? 0.29038188????? 0.24035060???
?????? ……

  评论这张
 
阅读(618)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017