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 = p
    if ~keyword_set(p) then p = 2.0

    distances = 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=p
    out_res_half = out_res / 2.0d
    x_min = min(xs_exist) - out_res_half
    x_max = max(xs_exist) + out_res_half
    y_min = min(ys_exist) - out_res_half
    y_max = max(ys_exist) + out_res_half
    cols = 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_exist

    for col_ix=0, cols - 1 do begin
        for row_ix=0, rows - 1 do begin
            if ~finite(target_interp[col_ix, row_ix], /nan) then continue
            x0 = x_min + col_ix * out_res + out_res_half
            y0 = y_max - row_ix * out_res - out_res_half
            target_interp[col_ix, row_ix] = _idw(x0, y0, targets_exist, xs_exist, ys_exist, p=p)
        endfor
    endfor
end

3.2 主程序

; @Author	: ChaoQiezi
; @Time		: 2023年11月7日-下午2:17:56
; @Email	: chaoqiezi.one@qq.com

; 该程序用于 对站点(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_dir
    out_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 begin
        target = 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, $  ; 地理坐标系为WGS84
            GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $
            GEOGANGULARUNITSGEOKEY: 9102}  ; 单位为度
        
        ; 输出
        out_path = out_dir + 'IDW_' + target_name + '.tiff'
        write_tiff, out_path, target_interp, geotiff=geo_info, /double
        print, target_name, ' IDW插值完成: ', timer_keep(), ' s', format='%s%s%0.2f%s'
    endforeach
end