ENVI_IDL:使用反距离权重法选取最近n个点插值(底层实现)并输出为Geotiff格式(效果等价于Arcgis中反距离权重插值)

发布时间：2022-08-19 13:16

1. 实验内容

反权重插值算法简介：

我实现了什么

2. 编程

代码部分

部分运行结果显示(与Arcgis对比)：

编辑

使用Arcgiss打开PM2.5 tiff文件(原本图像时黑白变化——》现经过了一些颜色处理)

这是Arcgis进行反距离权重插值的参数设置(为确保与代码结果能够比较，这里将输出像元大小更改为0.001，幂为2默认，选取的点数更改为11)：

这是Arcgis的的结果(颜色也经过了一些处理)：

这是ENVI_IDL代码的结果：

1. 实验内容

反权重插值算法简介：

我实现了什么

使用反距离权重插值法选取距离最近的n个点进行插值

另外这里面用到了新函数sort(),但是应该不难，可以看ENVI的help，也可以自己敲代码熟悉一下这个函数。

2. 编程

代码部分

pro week_six_test
  ; 使用反距离权重插值法选取距离最近的n个点进行插值
  
  ; 路径
  in_dir = 'D:/IDL_program/experiment_data/chapter_4/air_quality_data.csv'
  out_dir = 'D:/IDL_program/experiment_data/chapter_4/'
  
  ; 获取数据(csv里面第0、1列是经纬度信息、后面几列是诸如pm2.5、pm10等等数据)
  data = read_csv(in_dir, header=par_name)
  ; 获取列索引的相关信息
  data_index_name = strsplit(par_name, /extract)
  data_index_count = n_elements(data_index_name)
  ; 提取lon、lat数据
  lon = data.(0)
  lat = data.(1)
  ; 假定输出的分辨率
  resolution = 0.001  ; 其实分辨率是多少都无所谓的，取决于你需要多大分辨率，因为值都是预测出来的
  
  ; 寻找最大最小的经纬度
  lon_min = min(lon)
  lon_max = max(lon)
  lat_min = min(lat)
  lat_max = max(lat)
  
  ; 计算构建数组的行列数
  column = ceil((lon_max - lon_min) / resolution)
  row = ceil((lat_max - lat_min) / resolution)
  
  
  
  ; 计算每个像元的行列号
  column_pos = floor((lon - lon_min) / resolution)
  row_pos = floor((lat_max - lat) / resolution)
  
  ; 对data的后面几列数据(pm2.5\NO2\SO2等等等)进行空间插值
  for data_index_i = 2, data_index_count - 1 do begin
    ; 时间开始计入
    start_time = systime(1)
    
    ; 输出路径
    out_path = out_dir + 'air_qulity_' + data_index_name[data_index_i] + '.tiff'
    
    ; 创建存储数据的数组
    data_box = fltarr(column, row)
    
    ; 放入已知数据
    data_box[column_pos, row_pos] = data.(data_index_i)
    ; 已知数据(样本值)的个数
    data_count = n_elements(data.(data_index_i))
    
    ; 创建存储预测数据(输出数据)的数组
    data_box_out = fltarr(column, row)
    
    ; 使用反权重插值公式进行预测数据
    for column_i = 0, column - 1 do begin
      for row_i = 0, row - 1 do begin
        ; 首先判断目前值是为已知值(即不是0.0值)，如果是那么原样输出即可
        if data_box[column_i, row_i] eq 0.0 then begin
          distance_sum = 0.0
          value_sum = 0.0
          
          ; 目前该像元点的经纬度坐标(计算距离需要使用这些)
          temp_lon = lon_min + column_i * resolution
          temp_lat = lat_max - row_i * resolution
          
          ; 对已知样本值进行循环——》计算distance(也就是公式里面的Di的求和)
; 下面的方法是取所有点进行公式的运算
;          for data_i = 0, data_count - 1 do begin
;            Di =  sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
;            distance_sum += 1.0 / (Di ^ 2.0)  ; 这里公式里面的p我们一般取2，这里也是
;          endfor
;          for data_i = 0, data_count - 1 do begin
;            Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
;            value_sum += (1.0 / (Di ^ 2.0)) * data.(data_index_i)[data_i] / distance_sum
;          endfor
          
; 下面的方法是可以选取最近的n个点进行计算
          ; 由于我们可能会说有这么一个需求：取最近的n个点进行计算(这里会使用sort()函数)
          ; 由于我们获取最近的n个点，那么就需要所有点的距离，然后排序才知道
          ; 创建存储每个样本值所在像元与所求像元距离的数组
          distance_array = fltarr(data_count)
          ; 进入循环求得每一个样本值所在像元与所求像元地距离并存储
          for data_i = 0, data_count - 1 do begin  ; 还有就是下面的平方什么的最好都加上.0表示小数，否则有可能最后会以整数进行运算——》导致结果的小数部分被舍去
            Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
            distance_array[data_i] = Di
          endfor
          
          ; 然后对distance_array数组进行排序——》使用sort()函数，传入数组distance_array
          near_distance_index = sort(distance_array)
          ; 第一，sort()函数不会对传入的数组进行更改，第二，返回的也是一个数组，这个数组是传入数组中的元素按从小到大(元素值)排列的索引(实在不理解可以自己看help，也可以自己写几行代码看看是怎么回事(我就是这么干的))
          
          ; 这里最近距离的点的个数n的设置
          near_n = 11  ; 假如我们设置为11(当然你可以修改设置为其它值)
          
          ; 进行取near_n = 11的distance_sum的计算
          for near_i = 0, near_n - 1 do begin
            data_i = near_distance_index[near_i]
            Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
            distance_sum += 1.0 / (Di ^ 2.0)  ; 公式里面的p我这里选取默认2，(Arcgis里面默认是2)
          endfor
          ; 进行最终值(value_sum)的运算
          for near_i = 0, near_n - 1 do begin
            data_i = near_distance_index[near_i]
            Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
            value_sum += (1.0 / (Di ^ 2.0)) * data.(data_index_i)[data_i] / distance_sum
          endfor
          
          ; 输出预测的数据
          data_box_out[column_i, row_i] = value_sum
        endif else begin
          data_box_out[column_i, row_i] = data_box[column_i, row_i]
        endelse
      endfor
    endfor
    
    ; 地理信息geotiff填写
    geo_info={$
      MODELPIXELSCALETAG:[resolution,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中该列的数据输出为geotiff文件
    write_tiff, out_path, data_box_out, geotiff=geo_info, /float
    
    stop_time = systime(1)
    print, data_index_name[data_index_i] + '>>> ' + strcompress(string(stop_time - start_time)) + 's'
  endfor
  print, '>>> 程序运行完毕'
end