1. 实验内容：
如题所示，就是对Modis Swath产品进行均值运算
这里需要注意几点：
第一：
由于每一个Modis Swath数据的经纬度范围不一样(它不像之前的OMI数据一样是全球数据，之前求和在平均就好了)，这里需要现在知道我们现在手上拿到的Modis Swath数据的经纬度范围有多大，那么就需要遍历所有的MOdis Swath数据找到所有的文件里面中的最大、最小经纬度
第二：
正常求均值就可以(讲的有些模糊，看代码或许会更清楚)
就像之前求借OMI数据差不多，其他的操作的也是常规操作。
2. 编程
pro week_five_test ; 本程序主要解决如何对Modis Swath数据进行均值运算 ; 思路 ; 1. 遍历所有的文件找到经纬度的最大值最小值 ; 路径 in_path = 'D:\IDL_program\experiment_data\chapter_3\modis_swath\geo_out' out_path = in_path ; 输出的分辨率 out_resolution = 0.03 ; 获取in_path下的所有tiff文件的路径以及文件数量 file_path_array = file_search(in_path, '*.tiff', count=file_count) ; 循环每一幅tiff文件，找到最大最小的经纬度 ; 预设值 lon_min = 9999.0 lon_max = -9999.0 lat_min = 9999.0 lat_max = -9999.0 ; 进入循环 for file_i = 0, file_count - 1 do begin ; 获取该循环下的tiff文件的数据以及geotiff data = read_tiff(file_path_array[file_i], geotiff=geo_info) ; 传入tiff文件的路径，geotiff=返回该tiff文件的投影信息，这里用geo_info变量接收 ; 从geo_info中获取空间分辨率 geo_resolution = geo_info.(0) ; 第0行是空间分辨率 geo_resolution_x = geo_resolution[0] geo_resolution_y = geo_resolution[1] ; 从geo_info中获取左上角点的经纬度坐标 geo_tag = geo_info.(1) ; 获取第一行(索引为0开始)即是关于左上角点的x、y、z以及lon、lat、dem数值 ; 获取数据的行列数 data_size = size(data) data_column = data_size[1] data_row = data_size[2] ; 获取经纬度的最大最小值 temp_lon_min = geo_tag[3] temp_lat_max = geo_tag[4] temp_lon_max = temp_lon_min + data_column * geo_resolution_x temp_lat_min = temp_lat_max - data_row * geo_resolution_y ; 判断是否是最大值，如果是，那么替换 if temp_lon_min lt lon_min then lon_min = temp_lon_min if temp_lon_max gt lon_max then lon_max = temp_lon_max if temp_lat_min lt lat_min then lat_min = temp_lat_min if temp_lat_max gt lat_max then lat_max = temp_lat_max endfor ; 计算行列数 data_box_column = ceil((lon_max - lon_min) / out_resolution) data_box_row = ceil((lat_max - lat_min) / out_resolution) ; 均值运算求和 和 频次的存储数组 data_box_sum = fltarr(data_box_column, data_box_row) data_box_num = fltarr(data_box_column, data_box_row) ; 循环每一个文件进行均值运算 for file_i = 0, file_count - 1 do begin ; 获取tiff文件数据 data = read_tiff(file_path_array[file_i], geotiff=geo_info) ; 获取该文件的空间分辨率 data_resolution = geo_info.(0) data_resolution_x = data_resolution[0] data_resolution_y = data_resolution[1] ; 获取文件的size data_size = size(data) ; 获取文件数据的行列数 data_column = data_size[1] data_row = data_size[2] ; 获取geo_info的左上角点x、y、z数据以及对应的lon、lat、dem数据 geo_tag = geo_info.(1) ; 获取文件数据的左上角点的经纬度坐标 data_ul_lon = geo_tag[3] data_ul_lat = geo_tag[4] ; 获取tiff文件中每一个像元的经纬度 data_lon_array = fltarr(data_column, data_row) ; 用来存储经度信息的数组 data_lat_array = fltarr(data_column, data_row) ; 用来存储纬度信息的数组 for data_column_i = 0, data_column - 1 do begin data_lon_array[data_column_i, *] = data_ul_lon + data_resolution_x * data_column_i endfor for data_row_i = 0, data_row - 1 do begin data_lat_array[*, data_row_i] = data_ul_lat - data_resolution_y * data_row_i endfor ; 获取tiff文件的每个像元在box中的行列数、 data_box_column_array = floor((data_lon_array - lon_min) / out_resolution) ; 错误示范 ; 这里犯了重大错误，我想当然的类似的照写(对于纬度而言),仔细想想这显然是有问题的！！！！！！ ; data_box_row_array = floor((data_lat_array - lat_min) / out_resolution) ; 正确示范 data_box_row_array = floor((lat_max - data_lat_array) / out_resolution) ; 均值运算的前提准备 data_box_sum[data_box_column_array, data_box_row_array] += data data_box_num[data_box_column_array, data_box_row_array] += (data gt 0.0) endfor ; 将data_box_num中频次为0的更改为1，避免等会儿出现0/0的情况，当然如果出现那么也会得出结果为NaN(看你是要出现0，还是NaN) data_box_num = (data_box_num gt 0.0) * data_box_num + (data_box_num eq 0.0) ; 均值运算 data_box_aver = data_box_sum / data_box_num ; 地理信息结构体 geo_info={$ MODELPIXELSCALETAG:[out_resolution,out_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} ; data_box_aver = rotate(data_box_aver, 7) ; 输出 write_tiff, out_path + '_aver.tiff', data_box_aver, geotiff=geo_info, /float end
运行结果如此下：
使用ENVI打开处理好的aver.tiff文件
