AOD实践，modis数据下载，modis数据处理

modis数据下载-数据读取-重投影-拼接-均值

一、数据下载

1、Cygwin安装

Cygwin安装教程：https://blog.csdn.net/u010356768/article/details/90756742

1.2 数据采集

现提供遥感数据下载服务，主要是NASA数据，数据下载网站包括：

LAADSDAAC: https://ladsweb.modaps.eosdis.nasa.gov/
href="https://ladsweb.modaps.eosdis.nasa.gov/search/

EarthData: https://search.earthdata.nasa.gov/

2、C下载

使用Cygwin批量下载卫星数据，以MODIS数据为例
打开Cygwin64 Terminal，依次输入如下命令，每行命令输完点回车。
1、cd G:\trydata\data\
将路径切换至.sh文件所在的文件夹。斜体改为你自己存放.sh文件的路径。
2、chmod 777 5386246296-download.sh
将.sh文件变为可执行文件。斜体改为你自己的.sh文件名。
3、./5386246296-download.sh
执行.sh文件。./后的斜体改为你自己的.sh文件名。
输入Username和password，即EARTHDATA的用户名和密码，注意输入密码的时候窗口不会显示，输完直接点回车即可。

3、谷歌浏览器扩展下载

在Chrome网上应用商店下载Chrono插件

二、hdf数据读取（MATLAB）

hdr = read_envihdr('****.hdr');
Image = multibandread('****.dat',hdr.size,[hdr.format '=>double'],hdr.header_offset,hdr.interleave,hdr.machine);

如何从气溶胶产品MOD04提取550nm处的气溶胶厚度

如果做数据处理，最好直接写程序去做。matlab中有现成的子程序，hdftool，你可以点击需要输入的参数，下面接着会有输入语句的书写方法，你就照着写到.m文件中做循环就可以了。如果想看数据的内部情况也可以用hdfview
mod04的气溶胶数据都是块状分布，且每天不同时段的都是一个hdf文件，用matlab的hdftool读取某一个hdf文件我可以，但请问如何
将同一天不同时段的数据融合在一起生成逐日的数据？然后如何转换成带坐标系的tiff格式？
假如数据存放的地址为fpath = ‘E:\MODIS’;
file = dir([fpath ‘*.hdf’]); % 得到所有要读文件的名称
for i =1:numel(file)
data = []; %预设的要读取的变量名
fn = [fpath file(i,1).name];
% 下面就可以读取变量了
data = hdfread(fn,……)；
save([outpath matfn],‘data’); % outpath 为输出的文件地址 matfn是你想存成的文件名，此时存储的是mat文件
end

三、MOD04_L2批处理之重投影

MCTK手动应用

批处理

pro modis_mctkcompile_opt idl2e=envi(/headless)fn= dialog_pickfile(title="open the mod04 data",/directory) ;打开数据目录output_location = 'D:\data\data1\' ;输出路径files=file_search(fn,"*.hdf",count=nums)bridges = mctk_create_bridges()for i=0,nums-1 do begin ;每一个影像进行处理modis_swath_file=files[i]basename=file_basename(files[i]) ;获取输入影像文件名output_rootname=strmid(basename,0,32) ;获取指定范围名字swath_name = 'mod04';这个值自己根据自己数据设置，用hdfviewer打开后，会显示sd_names = ['AOD_550_Dark_Target_Deep_Blue_Combined'] ;；这是数据集名称，是一个字符串数组;  以下参数需要自己去看文档，一般保持不变 。官方地址：https://github.com/dawhite/MCTKout_method = 1output_projection = envi_proj_create(/geographic)interpolation_method = 0print,output_rootnameconvert_modis_data, in_file=modis_swath_file, $out_path=output_location, out_root=output_rootname, $swt_name=swath_name, sd_names=sd_names, $out_method=out_method, out_proj=output_projection, $interp_method=interpolation_method, /no_msg, $r_fid_array=r_fid_array, r_fname_array=r_fname_array, $bridges=bridges, msg=msgendformctk_destroy_bridges, bridges
end

四、MOD04_L2批处理之拼接

1、MODIS图像批量镶嵌拼接方法（IDL/ENVI）

（不会，有会得教教我！！！）

2、ENVI mosaic

2.1手动

【ENVI入门系列】09.图像镶嵌：http://blog.sina.com.cn/s/blog_764b1e9d0102v1p9.html
ENVI5.1无缝镶嵌工具（具体功能）
文件选择

2.2二次开发——拼接

PRO MOSAIC_BATCHCOMPILE_OPT IDL2; 启动ENVI 5.1e = ENVI()fn= dialog_pickfile(title="Select input scenes",/directory) ;打开数据目录files=file_search(fn,"*.dat",count=nums); output_location = 'D:\data\outdata3\' ;输出路径scenes = !NULL; 将每一个Raster放在一个Scenes中FOR i=0, N_ELEMENTS(files)-1 DO BEGINbasename=file_basename(files[i]) ;获取输入影像文件名name=strmid(basename,0,20) ;获取指定范围名字for n=0, 12 do beginm=i+nbasename_1=file_basename(files[m]) ;获取输入影像文件名name_1=strmid(basename_1,0,20) ;获取指定范围名字print,basename_1if name_1 ne name then breakraster = e.OpenRaster(files[m])scenes = [scenes, raster]endfori=m-1; 创建ENVIMosaicRaster对象mosaicRaster = ENVIMosaicRaster(scenes,$background = -999,$color_matching_method = 'histogram matching',$color_matching_stats = 'overlapping area',$feathering_distance = 20,$feathering_method = 'seamline',$resampling = 'bilinear',$seamline_method = 'none')；  newFile = ENVI_PICKFILE(title='Select output file', $ /output)
；  IF FILE_TEST(newFile) THEN FILE_DELETE, newFile; 设置输出路径output_location = 'G:\NDVI\mosaic\' + name + '.dat' ;输出路径IF FILE_TEST(output_location) THEN FILE_DELETE, output_location ; 输出镶嵌结果mosaicRaster.Export, newFile, 'ENVI'ENDFOR
END

五、裁剪

app商店下载安装批处理工具包

六、均值计算

pro modis_swath_averagestart_time= systime(1)input_directory= 'D:\study\AOD\Mosaicdata1\'output_directory=input_directoryfile_list=file_search(input_directory,'*.tif');路径file_n=n_elements(file_list)output_resolution=0.03 ;像元分辨率output_name=output_directory+'avr.tif';结果print,file_n;循环，读图像信息:经纬度lon_min=9999.0lon_max=-9999.0lat_min=9999.0lat_max=-9999.0for file_i=0,file_n-1 do begindata=read_tiff(file_list[file_i],geotiff=geo_info)data_size=size(data)data_col=data_size[1];范围:最大data_line=data_size[2]resolution_tag=geo_info.(0)geo_tag=geo_info.(1)temp_lon_min=geo_tag[3]temp_lon_max=temp_lon_min+data_col*resolution_tag[0]temp_lat_max=geo_tag[4]temp_lat_min=temp_lat_max -data_line*resolution_tag[1]if temp_lon_min lt lon_min then lon_min=temp_lon_minif temp_lon_max gt lon_max then lon_max=temp_lon_maxif temp_lat_min lt lat_min then lat_min=temp_lat_minif temp_lat_max gt lat_max then lat_max=temp_lat_maxendfor;经纬度数据存储data_box_geo_col=ceil((lon_max-lon_min)/output_resolution)data_box_geo_line=ceil((lat_max-lat_min)/output_resolution)data_box_geo_sum=fltarr(data_box_geo_col,data_box_geo_line)data_box_geo_num=fltarr(data_box_geo_col,data_box_geo_line);逐像元循环位置for file_i=0,file_n-1 do beginprint,file_list[file_i]data=read_tiff(file_list[file_i],geotiff=geo_info)data_size=size(data)data_col=data_size[1]data_line=data_size[2]resolution_tag=geo_info.(0)geo_tag=geo_info.(1)temp_lon_min=geo_tag[3]temp_lat_max=geo_tag[4]for data_col_i=0,data_col-1 do beginfor data_line_i=0,data_line-1 do begintemp_lon=temp_lon_min+data_col_i*resolution_tag[0]temp_lat=temp_lat_max-data_line_i*resolution_tag[1]data_box_col_pos=floor((temp_lon-lon_min)/output_resolution)data_box_line_pos=floor((lat_max-temp_lat)/output_resolution)if (data[data_col_i,data_line_i]eq 0.0)then continuedata_box_geo_sum[data_box_col_pos,data_box_line_pos]=data_box_geo_sum[data_box_col_pos,data_box_line_pos]+data[data_col_i,data_line_i]data_box_geo_num[data_box_col_pos,data_box_line_pos]=data_box_geo_num[data_box_col_pos,data_box_line_pos]+1.0endforendfor;col_start=floor((temp_lon_min-lon_min)/output_resolution);  line_start=floor((lat_max-temp_lat_max)/output_resolution);  data_box_geo_sum[col_start:col_start+data_col-1,line_start:line_start+data_line-1]+=data;  data_box_geo_num[col_start:col_start+data_col-1,line_start:line_start+data_line-1]+=(data gt 0.0)endfordata_box_geo_num=(data_box_geo_num gt 0.0)*data_box_geo_num+(data_box_geo_num eq 0.0)data_box_geo_avr=data_box_geo_sum/data_box_geo_numgeo_info={$MODELPIXELSCALETAG:[output_resolution,output_resolution,0.0],$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}write_tiff,output_name,data_box_geo_avr,geotiff=geo_info,/floatend_time=systime(1)print,'Time consuming: '+strcompress(string(end_time-start_time))
end

hdr转TIFF

七、tiff数据读取

[A,R] = geotiffread(filename)