ENVI IDL：如何基于气象站点数据进行反距离权重插值？

01 前言

仅仅练习，大可使用ArcGIS或者已经封装好的python模块进行插值，此处仅仅从底层理解如何从公式和代码理解反距离权重插值的过程，从而更深刻的理解IDL的使用和插值的理解。

02 函数说明

2.1 Read_CSV()函数

官方语法如下：

Result = READ_CSV( Filename [, COUNT=variable] [, HEADER=variable] [, MISSING_VALUE=value] [, N_TABLE_HEADER=value] [, NUM_RECORDS=value] [, RECORD_START=value] [, TABLE_HEADER=variable] [, TYPES=value] )

Filename表示读取的CSV文件的路径；
COUNT表示读取的CSV文件内表格的行数（不包含标签头即第一行）
HEADER表示读取的CSV文件内表头（以字符串数组存储表头信息，默认第一行记录为表头<如果有>）
MISSING_VALUE表示对于CSV文件内表格中的空值应该赋予何值呢？默认是赋予0。
N_TABLE_HEADER表示表头的行数，或许我们的表头不止一行，那么使用header就很难获取得到所有的表头信息，因此我们需要指定表头到底有多少行。一般与TABLE_HEADER连用，获取的多行表头返回给该参数，且其优先级高于header。
NUM_RECORDS表示读取的总行数，默认是所有行都读取。
RECORD_START表示开始读取的行的索引，默认从0开始(0为表头行)
TYPES传入各个列的数据类型(字符串数组形式，每一列的记录的数据类型）以下是各个数据类型的参数：

""表示该列的数据类型自动确定数据类型。

03 代码

3.1 封装的反距离权重插值函数

;+
;   函数用途:
;       IDW插值相关(私有函数), 用于单个像元值的插值计算
;   函数参数:
;       ···
;-
function _idw, x0, y0, targets_exist, xs_exist, ys_exist,  p = pif ~keyword_set(p) then p = 2.0distances = sqrt((x0 - xs_exist) ^ 2.0 + (y0 - ys_exist) ^ 2.0)distances_coef = total(1.0 / (distances ^ p))interp_target = total(targets_exist / ((distances ^ p) * distances_coef))return, interp_target
end;+
;   函数用途:
;       该函数基于少数点位进行反距离权重插值(IDW)生成指定范围的插值栅格矩阵
;   函数参数:
;       targets_exist: 插值的目标向量(数组形式)
;       xs_exist: 与目标向量对应的X坐标向量集(数组形式)
;       ys_exist: 与目标向量对应的Y坐标向量集(数组形式)
;       out_res: 插值后输出的分辨率大小
;       target_interp: 输出插值后的目标矩阵
;-
pro idw, targets_exist, xs_exist, ys_exist, out_res, target_interp, p=pout_res_half = out_res / 2.0dx_min = min(xs_exist) - out_res_halfx_max = max(xs_exist) + out_res_halfy_min = min(ys_exist) - out_res_halfy_max = max(ys_exist) + out_res_halfcols = ceil((x_max - x_min) / out_res)rows = ceil((y_max - y_min) / out_res)target_interp = make_array(cols, rows, /double, value=!values.F_NAN)existing_cols = floor((xs_exist - x_min) / out_res)existing_rows = floor((y_max - ys_exist) / out_res)target_interp[existing_cols, existing_rows] = targets_existfor col_ix=0, cols - 1 do beginfor row_ix=0, rows - 1 do beginif ~finite(target_interp[col_ix, row_ix], /nan) then continuex0 = x_min + col_ix * out_res + out_res_halfy0 = y_max - row_ix * out_res - out_res_halftarget_interp[col_ix, row_ix] = _idw(x0, y0, targets_exist, xs_exist, ys_exist, p=p)endforendfor
end

3.2 主程序

; @Author	: ChaoQiezi
; @Time		: 2023年11月7日-下午2:17:56
; @Email	: [e-mail]; 该程序用于 对站点(CSV)文件中的空气质量参数(多种污染物浓度)进行指定范围的插值; 主程序
pro idw_interp; 准备in_path = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week7\air_quality_data.csv\'out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week7\out_me\'if ~file_test(out_dir, /directory) then file_mkdir, out_dirout_res = 0.001d  ; 输出分辨率, 度(°)out_res_half = out_res / 2.0d; 读取ds = read_csv(in_path, count=count, header=header, missing_value=!values.F_NAN)lon = ds.(0)lat = ds.(1)targets_name = header[2:*]foreach target_name, targets_name, ix do begintarget = ds.(ix + 2)idw, target, lon, lat, out_res, target_interp; 地理结构体geo_info={$MODELPIXELSCALETAG: [out_res, out_res, 0.0], $  ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, min(lon) - out_res_half, max(lat) + out_res_half, 0.0], $  ; 角点信息GTMODELTYPEGEOKEY: 2, $  ; 设置为地理坐标系GTRASTERTYPEGEOKEY: 1, $  ; 像素的表示类型, 北上图像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $  ; 地理坐标系为WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102}  ; 单位为度; 输出out_path = out_dir + 'IDW_' + target_name + '.tiff'write_tiff, out_path, target_interp, geotiff=geo_info, /doubleprint, target_name, ' IDW插值完成: ', timer_keep(), ' s', format='%s%s%0.2f%s'endforeach
end