1. 实验内容：
批量拼接Modis Swath的逐日数据并输出为Geotiff格式
另外这里需要说一下：
有一个新的函数或许你没有学过: uniq()函数——》字面意思就是获取唯一值的索引，但是这个唯一值有点不一样。
函数解释：首先函数会返回一个数组arr
传入一个数组array，依次检索数组array的元素，每检索一个元素，获取该元素在array的索引。若检索下一个元素时，下一个元素与之前的元素不一样，那么函数会将上一个元素的索引放入到数组arr中；若检索下一个元素时，下一个元素与之前的元素一样，那么函数会继续检索下一个元素。怎么说呢？我举几个例子吧？
如果最后你还是没有理解uniq函数的意思，那么我只能说我解释不够好，你可以去官网查看该函数的解释。
最后，如果你已经理解了uniq函数，那么你可能会困惑这个函数的鸡肋，更会觉得代码里面使用uniq函数是否是不对的，这里我一开始也是比较纳闷的，但是需要多理解。最后对于代码uniq函数我只解释一点，代码可以使用uniq函数完全是因为重复值都是相邻的，如果不相邻，那么也是不可以用uniq函数的。
2. 编程
pro week_five_homework ; 本程序主要解决如何实现Modis Swath产品的每天数据的拼接(如果数据有重叠就均值处理) ; 路径 in_path = 'D:/IDL_program/experiment_data/chapter_3/modis_swath/geo_out' out_path = 'D:/IDL_program/experiment_data/chapter_3/modis_swath/mosaic' if file_test(out_path, /DIRECTORY) eq 0 then begin ; 没有该目录就创建 file_mkdir, out_path endif ; 如何寻找一天的数据不一定只有一张tiff文件，可能有多张，如何找到 ; 获取in_path下所有tiff文件的绝对路径 file_path_array = file_search(in_path, '*.tiff') ; 在in_path目录下寻找tiff文件,以数组形式返回找到的所有文件的绝对路径 ; 获取所有tiff文件的名称 file_name_array = file_basename(file_path_array) ; 获取所有tiff文件中名称的日期 file_date_array = strmid(file_name_array, 10, 7) ; 获取不重复日期的最后一个索引(也许你听的有点别扭，但你可以仔细阅读uniq()函数就明白了) file_date_uniq = uniq(file_date_array) ; 这里如果你了解了uniq()函数的作用，那么你会疑惑，why?这是因为我们获取的时候，相同日期就是连在一起的， 所以可以得到唯一值 ; 获取不重复日期的总个数 file_date_count = n_elements(file_date_uniq) ; 输出tiff文件的分辨率(x、y均是一样的) out_resolution = 0.03 ; 循环进行每天数据的拼接 for file_date_i = 0, file_date_count - 1 do begin ; 获取in_path中该循环日期下的所有文件的路径(以数组形式返回) today_file_path_array = file_search(in_path, '*A' + file_date_array[file_date_uniq[file_date_i]] + '*.tiff') ; 该日期下的所有文件的个数 today_file_count = n_elements(today_file_path_array) ; 由于file_search()函数若只有一个匹配结果，那么不会返回数组而是返回str，所以上面的n_elemen()函数对于str返回0 ; 所以这里做一点改变 if today_file_count eq 0 then begin today_file_count = 1 today_file_path_array = [today_file_path_array] endif ; 获取这一天的所有文件的最大最小经纬度 ; 预设值 lon_min = 9999.0 lon_max = -9999.0 lat_min = 9999.0 lat_max = -9999.0 for today_file_i = 0, today_file_count - 1 do begin ; 获取该文件的路径 today_file_path = today_file_path_array[today_file_i] ; 获取该文件的数据 today_file_data = read_tiff(today_file_path, geotiff=geo_info) ; 同时获取该tiff文件的地理信息(返回给了变量geo_info) ; 该数据的size today_file_data_size = size(today_file_data) ; 获取该数据的空间分辨率 today_file_resolution = geo_info.(0) ; 获取该数据左上角点的地理信息 today_file_loc = geo_info.(1) ; 获取该数据左上角点的经纬度坐标 temp_lon_min = today_file_loc[3] temp_lat_max = today_file_loc[4] temp_lon_max = temp_lon_min + today_file_resolution[0] * today_file_data_size[1] temp_lat_min = temp_lat_max - today_file_resolution[1] * today_file_data_size[2] ; 如果有更大或更小的经纬度，那么替换 if temp_lon_min lt lon_min then lon_min = temp_lon_min if temp_lat_max gt lat_max then lat_max = temp_lat_max 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 endfor ; today所有的文件需要的行列数 today_col = ceil((lon_max - lon_min) / out_resolution) today_row = ceil((lat_max - lat_min) / out_resolution) ; 下面这些是不需要的 ; ; 经纬度box初始化 ; today_lon_box = fltarr(today-col, today_row) ; today_lat_box = fltarr(today-col, today_row) ; ; 初始化值设定 ; today_lon_box[*, *] = -9999.0 ; today_lat_box[*, *] = -9999.0 ; 装数据的sum_box初始化 today_data_box = fltarr(today_col, today_row) ; num_box初始化 today_data_num_box = fltarr(today_col, today_row) ; 循环获取该日期下的所有文件的数据并进行拼接 for today_file_i = 0, today_file_count - 1 do begin ; 获取该文件路径 today_file_path = today_file_path_array[today_file_i] ; 获取该文件的数据 today_file_data = read_tiff(today_file_path, geotiff=geo_info) ; 该数据的size today_file_data_size = size(today_file_data) ; 获取该数据的空间分辨率 today_file_resolution = geo_info.(0) ; 获取该数据左上角点的地理信息 today_file_loc = geo_info.(1) ; 获取该数据左上角点的经纬度坐标 temp_lon_min = today_file_loc[3] temp_lat_max = today_file_loc[4] ; 初始化行列数 today_data_col_box = fltarr(today_file_data_size[1], today_file_data_size[2]) today_data_row_box = fltarr(today_file_data_size[1], today_file_data_size[2]) ; 该数据的在大box中的行列数 for x_i = 0, today_file_data_size[1] - 1 do begin today_data_col_box[x_i, *] = ((temp_lon_min + today_file_resolution[0] * x_i) - lon_min) / out_resolution endfor for y_i = 0, today_file_data_size[2] - 1 do begin today_data_row_box[*, y_i] = (lat_max - (temp_lat_max - today_file_resolution[1] * y_i)) / out_resolution endfor ; 求和 和 频次 计入 today_data_box[today_data_col_box, today_data_row_box] += today_file_data today_data_num_box[today_data_col_box, today_data_row_box] += (today_file_data gt 0) endfor ; 均值运算 ; today_data_num_box优化 today_data_num_box = (today_data_num_box gt 0) * today_data_num_box + (today_data_num_box eq 0) ; 运算 today_data_aver_box = today_data_box / today_data_num_box ; 地理信息填写 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} ; 输出 write_tiff, out_path + '/aver_' + file_date_array[file_date_uniq[file_date_i]] + '.tiff', today_data_aver_box, geotiff=geo_info, /float endfor end
运行结果：
某一张拼接效果：
最后我想说的是，由于这次想睡觉，时间比较紧，所以代码肯定有很多冗余的部分，譬如，如果某一天tiff文件只有一个，那么我没有输出也是直接进入循环。最后的图片虽然大致是一样的，但是还是有很小的偏移，我就懒得改动了，相信你们可以自己直接输出！
