MCD19A2 MAIAC AOD 数据处理（三）均值+同日镶嵌+重投影（Grid转经纬度）

本节包含以下三方面：

单景均值处理
同日镶嵌
重投影

一、均值
感觉上周均值写的比较乱，重新整理了一下：
使用的数据集：MCD19A2 MAIAC AOD 数据
目的：读取HDF文件AOD550数据集并与对应的QA数据集结合进行质量筛选，并对筛选后的数据（因为卫星观测次数不定，可近似理解为部分数据集为多波段），对其进行均值处理。

;读数据集
function hdf4_data_get,sd_id,sds_namesds_index=hdf_sd_nametoindex(sd_id,sds_name)sds_id=hdf_sd_select(sd_id,sds_index)hdf_sd_getdata,sds_id,datahdf_sd_endaccess,sds_idreturn,data
end;读数据集属性信息
function hdf4_attdata_get,sd_id,sds_name,att_namesds_index=hdf_sd_nametoindex(sd_id,sds_name)sds_id=hdf_sd_select(sd_id,sds_index)att_index=hdf_sd_attrfind(sds_id,att_name)hdf_sd_attrinfo,sds_id,att_index,data=att_datahdf_sd_endaccess,sds_idreturn,att_data
endpro QA_meancompile_opt idl2envi,/restore_base_save_filesenvi_batch_initinput_directory='D:/aerosol_text1/2001001_007/'output_directory='D:/aerosol_text1/2001001_007/average_AOD_out1/'directory_exist=file_test(output_directory,/directory)if(directory_exist eq 0)then beginfile_mkdir,output_directoryendiffile_list=file_search(input_directory,'*.hdf')file_n=n_elements(file_list);*****************************************************************;第一层for循环遍历整个HDF文件，对每个hdf文件进行处理for file_i=0,file_n-1 do beginstart_time=systime(1)sd_id=hdf_sd_start(file_list[file_i],/read)AOD_550_Data=hdf4_data_get(sd_id,'Optical_Depth_055')AOD_Scale_Factor=hdf4_attdata_get(sd_id,'Optical_Depth_055','scale_factor')AOD_QA=hdf4_data_get(sd_id,'AOD_QA')AOD_Size=size(AOD_550_Data,/dimensions)gindex=hdf_sd_attrfind(sd_id,'Orbit_amount')hdf_sd_attrinfo,sd_id,gindex,data=Orbit_amountAOD_average=replicate(-9999.0,AOD_Size[0],AOD_Size[1]);help,AOD_sizeAOD_Data=AOD_550_Data*AOD_Scale_Factor[0];获取AOD有效范围值AOD_Data_valid=(AOD_Data ge -0.1)*(AOD_Data le 5)*AOD_Data+((AOD_Data gt 5)+(AOD_Data lt -0.1))*(-9999.0);符合AOD有效值条件的保留原值，不符合条件的改为背景值hdf_sd_end,sd_id;判断有几个波段if Orbit_amount[0] gt 1 then beginAOD_select=fltarr(AOD_Size[0],AOD_Size[1],AOD_Size[2])for band_i=0,AOD_Size[2]-1 do beginQA_binary=AOD_QA[*,*,band_i].ToBinary()add='000000000000000';自己定义一个15位的字符串，为了填补高位为0add_string=add+QA_binary;把两个字符串连接起来;flag1=strmid(add_string,2,3,/REVERSE_OFFSET);取001无云flag2=strmid(add_string,11,4,/REVERSE_OFFSET);取0000高质量;QA_valid=((flag1 and 1) eq 1 )*((flag2 and 0) eq 0)QA_valid=(flag2 and 0) eq 0help,QA_validQA=reform(QA_valid,AOD_size[1],AOD_size[1])AOD_select[*,*,band_i]=AOD_Data_valid[*,*,band_i]*(QA eq 1)+(QA eq 0)*(-9999.0);符合QA筛选条件的保留原值，不符合条件的改为背景值endforfor i=0,AOD_Size[0]-1 do beginfor j=0,AOD_Size[1]-1 do begin;创建一个二维浮点型数组data_effect=fltarr(AOD_Size[2])for k=0,AOD_Size[2]-1 do begindata_effect[k]=AOD_select[i,j,k]endforpos=where((data_effect ne -9999),count)if count eq 0 then beginAOD_average[i,j]=AOD_average[i,j];endif else beginAOD_average[i,j]=mean(data_effect[pos])endelseendforendforendif else begin;AOD_select=fltarr(AOD_Size[0],AOD_Size[1])QA_binary=AOD_QA.ToBinary()add='000000000000000';自己定义一个15位的字符串，为了填补高位为0add_string=add+QA_binary;把两个字符串连接起来;flag1=strmid(add_string,2,3,/REVERSE_OFFSET);取001flag2=strmid(add_string,11,4,/REVERSE_OFFSET);取0000;QA_valid=((flag1 and 1) eq 1 )*((flag2 and 0) eq 0)QA_valid=(flag2 and 0) eq 0help,QA_validQA=reform(QA_valid,AOD_size[0],AOD_size[1])AOD_average=AOD_Data_valid*(QA eq 1)+(QA eq 0)*(-9999.0)endelseresult_name=output_directory+strmid(file_basename(file_list[file_i],'.hdf'),8,15)+'QA_AOD'+'.tiff'write_tiff,result_name,AOD_average,/floatend_time=systime(1)print,'Reprojection time consuming of file'+strmid(file_basename(file_list[file_i],'.hdf'),8,15)+':'+strcompress(string(end_time-start_time))+'s.'endforend

二、同日镶嵌
根据均值输出文件命名中的"h、v"范围确定镶嵌范围，根据日期确定同日镶嵌的单景影像。
使用数据：单景均值影像（TIFF格式）

pro Mosaiccompile_opt idl2envi,/restore_base_save_filesenvi_batch_initinput_directory='D:/aerosol_text1/2001001_007/average_AOD_out1/'output_directory='D:/aerosol_text1/2001001_007/mosaic1_out/'directory_exist=file_test(output_directory,/directory)if(directory_exist eq 0)then beginfile_mkdir,output_directoryendiffile_list=file_search(input_directory,'*.tiff',count=bnum)file_n=n_elements(file_list) dims=size(read_tiff(file_list[0]),/dimensions) size_x=dims[0];数据集本身的列数size_y=dims[1];数据集本身的行数;for file_i=0,file_n-1 do begin;print,file_list[file_i]file_prefix=(strmid(file_basename(file_list),0,8)).uniq() ;取A+年积日，uniq返回'A+年积日'不同的值构成的集合;对不同的天做循环for prefix_i=0,n_elements(file_prefix)-1 do begin;n_elements(file_prefix)-1代表不同天数的个数，即我们要合成多少张图像，又或者说天数start_time=systime(1)print,file_prefix[prefix_i]year=fix(strmid(file_prefix[prefix_i],1,4))doy=fix(strmid(file_prefix[prefix_i],5,3))date_julian=imsl_datetodays(31,12,year-1)+doyimsl_daystodate,date_julian,day,month,yearresult_name=output_directory+string(year,format='(I04)')+string(month,format='(I02)')+string(day,format='(I02)')+'_mosaic.tiff'AOD_file_list=file_search(input_directory+file_prefix[prefix_i]+'*.tiff',count=day_n);寻找同一天的tiff数据h_value=fix(strmid(file_basename(AOD_file_list),10,2))v_value=fix(strmid(file_basename(AOD_file_list),13,2))h_min=fix(min(h_value),type=3)h_max=fix(max(h_value),type=3)v_min=fix(min(v_value),type=3)v_max=fix(max(v_value),type=3)h_num=h_max-h_min+1Lv_num=v_max-v_min+1Lcols=h_num*size_xlines=v_num*size_xgrid_box=replicate(-9999.0,cols,lines);构建能够填入一天内所有数据的矩阵;对同一天的数据进行循环,并赋值给新组好的GRID矩阵for file_i=0,day_n-1 do begin;print,AOD_file_list[file_i]h_value=fix(strmid(file_basename(AOD_file_list[file_i]),10,2))v_value=fix(strmid(file_basename(AOD_file_list[file_i]),13,2))data_temp=read_tiff(AOD_file_list[file_i])col_start_pos=(h_value-h_min)*size_xcol_end_pos=col_start_pos+size_x-1line_start_pos=(v_value-v_min)*size_yline_end_pos=line_start_pos+size_y-1grid_box[col_start_pos:col_end_pos,line_start_pos:line_end_pos]=data_tempendforend_time=systime(1)write_tiff,result_name,grid_box,/floatprint,'mosaic time consuming of file'+string(year,format='(I04)')+string(month,format='(I02)')+string(day,format='(I02)')+':'+strcompress(string(end_time-start_time))+'s.'endforend

三、重投影(Grid转经纬度)
核心函数：map_proj_init()、map_proj_inverse()
所用数据：镶嵌.TIFF
代码如下：

pro Grid_to_Geocompile_opt idl2envi,/restore_base_save_filesenvi_batch_initinput_directory='D:/aerosol_text1/2001001_007/Terra_out/day_mosaic/';已经镶嵌过的TIFF影像output_directory='D:/aerosol_text1/2001001_007/Terra_out/day_mosaic/day_mosaic_geo/'directory_exist=file_test(output_directory,/directory)if(directory_exist eq 0)then beginfile_mkdir,output_directoryendifstart_time=systime(1)file_list=file_search(input_directory,'*.tiff',count=bnum)file_n=n_elements(file_list)dims=size(read_tiff(file_list[0]),/dimensions)cols=dims[0];数据集本身的列数lines=dims[1];数据集本身的行数corner=[5559752.598333,6671703.118,13343406.236,0.0];角点信息，我感觉利用函数读比较浪费时间，就直接写出来了：左上角点，右下角点Xmin\Ymax\Xmax\Ymingrid_resolution=(corner[2]-corner[0])/cols;计算格网分辨率;创建和源图像等大的数组，分别用来存放像元X和Y坐标信息grid_col_pos=fltarr(cols,lines)grid_line_pos=fltarr(cols,lines)for col_i=0,cols-1 do begingrid_col_pos[col_i,*]=corner[0]+grid_resolution*col_iendforfor line_i=0,lines-1 do begingrid_line_pos[*,line_i]=corner[1]-grid_resolution*line_iendfor;坐标转换参数的获取，里面的数值是固定的，是专门用于Grid转经纬度的sin_prj=map_proj_init('sinusoidal',/gctp,sphere_radius=6371007.181,center_longitude=0.0,false_easting=0.0,false_northing=0.0);坐标转换，map_proj_inverse专门用于向经纬度转换的函数;结果是与格网坐标一一对应的经纬度坐标,第一列是经度，第二列是纬度geo_loc=map_proj_inverse(grid_col_pos,grid_line_pos,map_structure=sin_prj);存储经纬度并获取经纬度图像的经纬度范围geo_x=geo_loc[0,*];经度geo_y=geo_loc[1,*];纬度lon_min=min(abs(geo_x));中国区最东边临界经度是180，可能是因为我没有区分卫星和时刻？我不加绝对值通过map_proj_inverse函数转换出来的经度范围是-180到180。加了绝对值，经度范围就正常了，可以试试不加与我的情况是否一样lon_max=max(geo_x)lat_min=min(geo_y)lat_max=max(geo_y);设置经纬度格网分辨率geo_resolution=0.01;分辨率的设置老师说要合理一些，最好和原图像像元覆盖范围差不多，咱们可以看一看MTCK做完的分辨率进行参考。;获取新图像的范围geo_col=ceil((lon_max-lon_min)/geo_resolution)geo_line=ceil((lat_max-lat_min)/geo_resolution)data_geo=replicate(-9999.0,geo_col,geo_line);设置背景值为-9999for file_i=0,file_n-1 do begin;start_time=systime(1);新图像有效值的填充及平滑处理;建立新旧图像坐标联系，获取像素值geo_col_pos=floor((geo_x-lon_min)/geo_resolution);新图像的列号geo_line_pos=floor((lat_max-geo_y)/geo_resolution);新图像的行号data_geo[geo_col_pos,geo_line_pos]=read_tiff(file_list[file_i]);平滑处理，填补孔隙;geo_out=fltarr(geo_col,geo_line)geo_out=data_geo;采取的方式为3*3窗口均值平滑处理中间的空洞，这个速度特别慢，不知道怎么修改快一点for geo_col_i=1,geo_col-2 do beginfor geo_line_i=1,geo_line-2 do begin      if data_geo[geo_col_i,geo_line_i] eq -9999.0 then begintemp_window=data_geo[geo_col_i-1:geo_col_i+1,geo_line_i-1:geo_line_i+1]temp_window_num=total(temp_window gt -0.1)if(temp_window_num gt 3) then begintemp_window=(temp_window gt -0.1)*temp_windowtemp_window_sum=total(temp_window)interpol_value=temp_window_sum/temp_window_numgeo_out[geo_col_i,geo_line_i]=interpol_valueendif else begincontinueendelseendifendforendforgeo_info={$MODELPIXELSCALETAG:[0.01,0.01,0.0],$;0.01是输出的经纬度分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$;左上角点经纬度信息(自西向东)GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}result_name=output_directory+string(file_basename(file_list[file_i],'.tiff'))+'_geo.tiff'write_tiff,result_name,geo_out,/float,geotiff=geo_infoend_time=systime(1)print,'Reprojection time consuming of file'+':'+strcompress(string(end_time-start_time))+'s.'endfor
end

在我学习过程中，特别致谢在各平台无私奉献的前辈老师们。没有前辈们的无私分享观点，我将永远是个菜鸡，虽然现在还是很菜，但是有一些进步了。
不积跬步无以至千里，不积小流无以成江海。
下一个目标：
分时刻单波段分离（已实现，下周再发）
AERONET数据下载（不知道有没有批处理的方法）
精度验证（目前思路尚不清晰）
有没有哪位前辈可以指点一二。

MCD19A2 MAIAC AOD 数据处理（三）均值+同日镶嵌+重投影（Grid转经纬度）相关推荐

MCD19A2 MAIAC AOD 数据处理（二）质量控制+均值计算
2.MCD19A2数据质量控制及日均值计算好几天没更新了,为什么呢?因为突然发现自己上一篇的方向搞错了,不应该用MCTK插件,不是说不能用,而是说不方便,应该利用IDL编程进行相关操作,可以调用其中 ...
MCD19A2 MAIAC AOD数据处理（一）MTCK插件批处理
1.MTCK插件批处理MCD19A2数据 MTCK插件介绍:将数据进行重投影,从正弦投影转化为WGS-84投影系统,且数据格式由.HDF转化为了.dat格式. 缺点: 运行速度超级慢
三本毕业论文查重吗？
现在大部分高校都需要进行论文查重,不管是一本.二本还是三本,我们需要认真对待论文,不要抄袭和复制,今天给大家讲解一下三本毕业论文查重吗? 三本毕业论文查重吗? 如今有许多人对三本大学有些偏见,认为三本 ...
【MaskRCNN】源码系列一：train数据处理三
目录 data_generator(接着数据处理二,data_generator中还有一部分没讲) build_rpn_targets data_generator(接着数据处理二,data_gene ...
python Excel数据处理三兄弟：xlrd/xlwt/xlutils！
常规的Excel数据处理中,就是对Excel数据文件的读/写/文件对象操作. 通过对应的python非标准库xlrd/xlwt/xlutils,来实现具体的数据处理业务逻辑. 在复杂的Excel业务数 ...
Camera - dump 预览帧数据处理(三)
应用层处理CPU.GPU框架下TextureView.GLSurfaceView所对应的Surface获取到的ImageReader数据以及dump相关的方法前面已经讲了,本次主要以展讯平台的HAL层 ...
图像处理(三) 均值滤波与中值滤波的对比
均值滤波与中值滤波的对比均值滤波与中值滤波实验对比 matlab 代码均值滤波与中值滤波均值滤波是典型的线性滤波算法,均值滤波是对目标像素及周边像素取平均值后再填回目标像素来实现滤波目的的方法 ...
视频图像数据处理三：将yuv420视频图像转换为灰度图像
文章目录函数代码测试用例下载本文介绍了将yuv420视频图像转换为灰度图像的方法,附有详细的代码和图像示例.文中yuv420文件需要使用yuv/rgb播放器才能查看,参考播放器可选择雷神推荐的 ...
Android View详解（三）视图状态及重绘流程分析
转载:http://blog.csdn.net/guolin_blog/article/details/17045157 [本文出自郭霖的博客] 在前面一篇文章中,我带着大家一起从源码的层面上分析了视 ...

MCD19A2 MAIAC AOD 数据处理（三）均值+同日镶嵌+重投影（Grid转经纬度）

MCD19A2 MAIAC AOD 数据处理（三）均值+同日镶嵌+重投影（Grid转经纬度）相关推荐

最新文章

热门文章