IDL 解析葵花8Himawari-8标准数据（HSD），辐射定标、重投影、裁剪

代码解析全云盘FLDK全波段（16个）的数据并进行定标、重投影和裁剪，重投影用的GLT生成方法参考之前的博客https://blog.csdn.net/qq_33339770/article/details/102957857
输入的是包含全部DAT数据的文件夹
输出的是与NC相同范围的2个影像（1000m和2000m各一景），其中Band3的500m分辨率影像被我采样到1000m

PRO Himawari8COMPILE_OPT idl2start=systime(1)inputfolder='***\HS_H08_20190916_0300_FLDK';获取影像时间basename=file_basename(inputfolder);获取1km和2km定标后影像数据imgdata=H8_Preprocess(inputfolder=inputfolder) R10_data=imgdata['R10']R20_data=imgdata['R20'];销毁哈希表OBJ_DESTROY,imgdata;应用GTL重投影并裁剪outR10=GLTproject(inputdata=R10_data,results=results,timename=timename)outR20=GLTproject(inputdata=R20_data,results=results,timename=timename)print,systime(1)-startEND
;对合成的影像应用GLT几何校正并裁减成NC的范围[80,60,200,-60]
Function GLTproject,inputdata=inputdata,results=results,timename=timenameCOMPILE_OPT idl2e=envi(/headless)inputraster=ENVIRaster(inputdata,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))inputraster.Savebands=inputraster.NBANDScols=inputraster.NCOLUMNSif cols eq 11000 then beginres='1000'gltfile='***\GLT_1000M.dat'endif else if cols eq 5500 then beginres='2000'gltfile='***\GLT_2000M.dat'endifgltraster=e.OpenRaster(gltfile)glt_id=ENVIRastertoFID(gltraster)ENVI_FILE_QUERY,glt_id,dims=glt_dims;应用GLT校正fid=ENVIRasterToFID(inputraster)out_name=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True')ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',FID=fid,GLT_DIMS=glt_dims,GLT_FID=glt_id,OUT_NAME=out_name,pos=lindgen(bands),r_fid=out_fid;裁剪输出outraster=ENVIFIDtoRaster(out_fid)SpatialRef = outraster.SPATIALREFSpatialRef.ConvertLonLatToMap, 80.0, 60.0, MapX, MapYSpatialRef.ConvertLonLatToMap, 200.0, -60.0, MapX2, MapY2subraster = outraster.Subset(SPATIALREF=SpatialRef,SUB_RECT=[MapX, MapY2, MapX2, MapY])subraster.Export,'******.dat','ENVI'gltraster.Closeinputraster.Closeoutraster.Closesubraster.Closereturn,outtiff
END;--解析数据并将500m重采样至1000m,合成两个数据集--
Function H8_Preprocess,inputfolder=inputfolderCOMPILE_OPT idl2R10_data=make_array(11000,11000,4,/FLOAT)R20_data=make_array(5500,5500,12,/FLOAT)bandnames=['01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16'];按波段查找文件for i=0,15 do beginband_files=file_search(inputfolder,'*B'+bandnames[i]+'_FLDK*.DAT',/test_regular);必须找到全部10个部分的数据if N_ELEMENTS(band_files) EQ 10 THEN BEGINband_data=[]seg=['01','02','03','04','05','06','07','08','09','10']for j=0,9 do begininfile=band_files[j]segment=read_hsd(inputfile=infile[0])imgdata=segment['pixels']imgtime=segment['time'];后续制作大气校正需要的太阳方位角、天顶角使用，代码后续补充sun_pos=segment['sun_pos']band_data=[[band_data],[imgdata]];销毁哈希表OBJ_DESTROY,segmentendforif i lt 2 then beginR10_data[*,*,i]=band_dataendif else if i eq 2 then begine=envi(/headless);把B03五百米分辨率重采样band_raster=ENVIRaster(band_data,URI=e.GetTemporaryFilename(CLEANUP_ON_EXIT='True'))band_raster.SavesampleRaster=ENVIResampleRaster(band_raster, DIMENSIONS=[11000,11000]);, METHOD='Bilinear')band_data=sampleRaster.GetData()sampleRaster.Closeband_raster.Closee.CloseR10_data[*,*,i]=band_dataendif else if i lt 4 then beginR10_data[*,*,i]=band_dataendif else beginR20_data[*,*,i-4]=band_dataendelseendifendforoutdata=hash('R10',R10_data,'R20',R20_data)return,outdata
END;读取HSD数据并定标为albedo/BT
Function read_hsd,inputfile=inputfileCOMPILE_OPT idl2spli=strsplit(file_basename(inputfile),'_');获取分辨率resolution=fix(strmid(file_basename(inputfile),spli[6]+1,2))if resolution eq 5 then begincols=22000rows=2200endif else if resolution eq 10 then begincols=11000rows=1100endif else if resolution eq 20 then begincols=5500rows=550endifopenr,h8_lun,inputfile,/get_lun;获取Observation start timepoint_lun,h8_lun,46imgtime=dblarr(1);R8readu,h8_lun,imgtime;获取Total header lengthpoint_lun,h8_lun,70header_length=uintarr(1);I4readu,h8_lun,header_length;获取影像point_lun,h8_lun,header_lengthimgdata=uintarr(cols,rows)readu,h8_lun,imgdata ;获取Sun's positionpoint_lun,h8_lun,510sun_pos=dblarr(3);R8readu,h8_lun,sun_pos;获取Band numberpoint_lun,h8_lun,601band=uintarr(1);I2readu,h8_lun,band  ;获取Gainpoint_lun,h8_lun,617Gain=dblarr(1);R8readu,h8_lun,Gain;获取Offsetpoint_lun,h8_lun,625Offset=dblarr(1);R8readu,h8_lun,Offset;计算radianceradiance=imgdata*Gain[0]+Offset[0]if band le 6 then begin;计算反射率;获取radiance to albedopoint_lun,h8_lun,633cc=dblarr(1);R8readu,h8_lun,cc;计算反射率albedo=radiance*cc[0]   ;返回值outdataoutdata=albedoendif else begin;计算亮温;获取Central wave lengthpoint_lun,h8_lun,603wv=dblarr(1);R8readu,h8_lun,wv;获取radiance to brightness temperature(c0)point_lun,h8_lun,633c0=dblarr(1);R8readu,h8_lun,c0;获取radiance to brightness temperature(c1)point_lun,h8_lun,641c1=dblarr(1);R8readu,h8_lun,c1;获取radiance to brightness temperature(c2)point_lun,h8_lun,649c2=dblarr(1);R8readu,h8_lun,c2;获取Speed of light (c)point_lun,h8_lun,681c=dblarr(1);R8readu,h8_lun,c    ;获取Planck constant (h)point_lun,h8_lun,689h=dblarr(1);R8readu,h8_lun,h;获取Boltzmann constant(k)point_lun,h8_lun,697k=dblarr(1);R8readu,h8_lun,k;计算亮温wv=wv[0]*1e-6rad=radiance*1e6Te=h[0]*c[0]/k[0]/wv[0]/(ALOG(2*h[0]*c[0]^2/(wv[0]^5*rad)+1))BT=c0[0]+c1[0]*Te+c2[0]*Te^2;返回值 outdataoutdata=BTendelsefree_lun,h8_lun ;返回：albedo/BT,时间,太阳坐标out=hash('pixels',outdata,'time',imgtime,'sun_pos',sun_pos)return,out;dataEND

IDL 解析葵花8Himawari-8标准数据（HSD），辐射定标、重投影、裁剪相关推荐

ENVI软件对Landsat-8数据进行辐射定标、大气校正、提取NDVI、估算植被覆盖度等操作【图说GIS】
目录一.前言二.数据下载第一步访问网站第二步登录账户第三步点击高级检索, 下载数据三.软件下载及安装四.辐射定标第一步打开数据第二步辐射定标五.计算区域平均高程第一步 ...
【python】批量实现modis数据的辐射定标，大气校正及地形校正
批量实现modis数据的辐射定标,大气校正及地形校正 1.定义辐射定标函数 def radiance_cal(band, gain, bias, scale_factor):rad = gain * ...
FY-4A图像数据的辐射定标和几何校正
1.1 FY-4A AGRI图像数据(辐射定标.几何校正) (1)辐射定标辐射定标就是把星上扫描辐射仪输出的计数值转换为要探测的物理量的绝对值.对于可见光一般转换为灰度和反照率,对于红外和水汽通道转 ...
利用python进行Landsat8数据的辐射定标
上次用于处理Landsat8数据的类,增加了辐射定标的功能.目前是用于可见光波段的辐射定标. 这段辐射定标的代码,首先从元文件中读取定标参数,然后进行辐射定标.对NaN数据做了处理,使其一直保持为Na ...
中秋节祝福程序源代码分享：土地分类数据阈值筛选和重投影分类
今年的中秋又要到啦,我的中秋节祝福程序源代码分享:过什么节,代码走起! 今天主要给大家介绍的是关于一个如何对影像重投影然后获取特定阈值情况下的影像先看结果: 原始影像重投影后的影像: 这里用到的数 ...
Linux系统下Sen2Cor对Sentinel哨兵2号遥感数据预处理(辐射定标和大气校正）Sen2Cor下载，使用
Sen2Cor插件的作用是对哨兵2号数据进行大气校正和辐射定标,将L1C数据,处理成为L2A数据.网上用window的教程有很多,这是linux的下载安装使用Sen2Sor的教程,我的系统是ubunt ...
利用MRT对MODIS数据进行批量重投影+批量波段合成
写在前面: 官方渠道已经下载不到MRT了,为什么呢?退休了呗. LP DAAC - The downloadable MODIS Reprojection Tool (MRT) and MRTSwat ...
php写的接口返回数据的页面,PHP怎么解析 WEBSERVICES接口返回的数据
PHP如何解析 WEBSERVICES接口返回的数据 WEBSERVICE服务器,他们用.net写了个接口 http://XXX.XXX.XXX.XX/Info.asmx?wsdl 我用PHP调用这个 ...
retrofit 解析百度地图api 返回数据_新版百度地图建筑数据含高度解析
注:本文中所述内容仅作为研究用途,不包括任何技术细节,不提供数据爬取.非法攻击的咨询和支持.另:百度建筑轮廓数据中存在明显签名式错误,下载使用容易被追求侵权. 矢量瓦片:爬取百度地图必须要了解矢量瓦片 ...
retrofit 解析百度地图api 返回数据_阿里二面：关于 Retrofit 你知道多少？看完你的offer稳了
一.整体思路从使用方法出发,首先是怎么使用,其次是我们使用的功能在内部是如何实现的, 实现方案上有什么技巧,有什么范式.全文基本上是对 Retrofit 源码的一个分析与导读,非常建议大家下载 R ...

IDL 解析葵花8Himawari-8标准数据（HSD），辐射定标、重投影、裁剪

IDL 解析葵花8Himawari-8标准数据（HSD），辐射定标、重投影、裁剪相关推荐

最新文章

热门文章