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
部分运行结果显示(与Arcgis对比)：
使用Arcgiss打开PM2.5 tiff文件(原本图像时黑白变化——》现经过了一些颜色处理)
这是Arcgis进行反距离权重插值的参数设置(为确保与代码结果能够比较，这里将输出像元大小更改为0.001，幂为2默认， 选取的点数更改为11)：
这是Arcgis的的结果(颜色也经过了一些处理)：
这是ENVI_IDL代码的结果：
可以发现，两者几乎没有区别，这说明我们的代码是有效的@！
