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文件 

——————————————————————————————————————————

——————————————————————————————————————————

 我是炒茄子，谢谢大家！