ENVI IDL：如何将txt文本文件转化为GeoTIFF文件？

01 前言

此处的文本文件形式如下：

里面包含了众多点位信息（不是站点数据），我们需要依据上述点的经纬度信息放到对应位置的像素点位置，放置完后如下：

可以发现，还存在部分缺失值，我们还需要进行缺失值的填补。

02 文本文件的读取

IDL读取文本文件还是不够方便，稍微封装了一下。

;+
;   函数用途:
;       用于读取文本文件
;   函数参数:
;       txt_path: 文本文件的路径
;       ds: 读取的输出数据集(不含表头)
;       header(关键字参数): 读取的输出表头
;       separator(关键字参数): 分隔符,默认空白符 
;-
pro read_txt, txt_path, ds, header=header, separator=separatorif ~keyword_set(separator) then separator = " "; 读取和检索openr, 1, txt_path  ; 打开文本文件; 是否指定输出的headerif ~arg_present(header) then beginskip_lun, 1, 0endif else beginheader = ''readf, 1, header  ; 读取表头header = strsplit(header, separator, /extract)endelse; 读取和处理数据集ds = strarr(file_lines(txt_path) - 1)readf, 1, dsds = list(ds, /extract)  ; 字符串数组转化为列表ds = ds.map(lambda(e, separator: double(strsplit(e, separator, /extract))), separator)ds = ds.toarray()  ; 列表转化为数组free_lun, 1
end

还好，算是可以用的程度了。

03 栅格矩阵的放置

于是乎，我们开始进行文件的读取和转数组。

pro txt_to_tiff; 准备in_path = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week6\2013_year_aop.txt'out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week6\out_tif_by_txt_me\'if ~file_test(out_dir, /directory) then file_mkdir, out_dirout_res = 0.18dout_res_half = out_res / 2.0d; 读取和检索read_txt, in_path, ds, header=headerlon = ds[*, 0]lat = ds[*, 1]ds = ds[*, 2:*]header = header[2:*]lon_min = min(lon) - out_res_halflon_max = max(lon) + out_res_halflat_min = min(lat) - out_res_halflat_max = max(lat) + out_res_halfcols = ceil((lon_max - lon_min) / out_res)rows = ceil((lat_max - lat_min) / out_res)lon_cols = floor((lon - lon_min) / out_res)lat_rows = floor((lat_max - lat) / out_res)foreach header_ele, header, header_ix do begintarget = make_array(cols, rows, value=!values.F_NAN)target[lon_cols, lat_rows] = ds[*, header_ix]; 填充window_interp, target, target_interp, interp=2; 输出out_path = out_dir + header_ele + '.tiff'write_img, out_path, target_interp, out_res, lon_min, lat_maxendforeach
end

在循环中，可以发现，使用了自定义的window_interp函数对target栅格矩阵进行缺失值的填补。

关于window_interp函数的定义由于封装的比较多，叠的比较层数比较多，阅读稍微有困难。学python的时候我是真的讨厌那些一个简单的功能的非要定义一个类，类又叠类，方法重组，来回找实现方法，来回折腾，本身功能不算特别复杂，但是被这么一折腾反而给阅读代码带来困难。

然而，我现在还是成为了他们。But 我将尽量让代码的逻辑清晰可见，不过分抽象，如有必要我抽出部分功能进行整合，避免使用过于复杂的代码框架搭建简单的功能。

以下是关于填补缺失值的封装函数，主要基于滑动窗口实现，包括滑动窗口均值填补和最近邻填补。
涉及自定义函数：window_interp、padding、interp_nearest、meshgrid。

;+
;   函数用途:
;       用于对二维数组进行边界填充
;   函数参数:
;       array: 用于边界填充的数组
;       pad_size: 单边(上下左右)填充的大小
;       pad_value(关键字参数: NAN): 填充的数值
;-
function padding, array, pad_size, pad_value=pad_valueif ~keyword_set(pad_value) then pad_value = !values.F_NAN; 获取基本信息ds_size = size(array, /dimensions)ds_type = size(array, /type)ds_size += pad_size * 2; padpad_array = make_array(ds_size, type=ds_type, value=pad_value)pad_array[pad_size:(ds_size[0] - pad_size - 1), $pad_size:(ds_size[1] - pad_size - 1)] = arrayreturn, pad_array
end;+
;   函数用途:
;       用于生成行列号格网矩阵
;   函数参数:
;       cols_n: 列数
;       rows_n: 行数
;-
function meshgrid, cols_n, rows_nwindow_cols = rebin(findgen(cols_n, 1), cols_n, rows_n)window_rows = rebin(findgen(1, rows_n), cols_n, rows_n)return, list(window_cols, window_rows)
endfunction interp_nearest, window_ds; 获取数组尺寸window_size = size(window_ds, /dimensions)cols_n = window_size[0]rows_n = window_size[1]; 生成行列号矩阵cols_rows = meshgrid(cols_n, rows_n)cols = cols_rows[0]rows = cols_rows[1]; 计算距离矩阵center_col = cols_n / 2center_row = rows_n / 2distance = sqrt((cols - center_col) ^ 2.0 +(rows - center_row) ^ 2.0)invalid_pos = where(finite(window_ds, /nan))distance[invalid_pos] = !values.F_NANinterp_value = (window_ds[where(distance eq min(distance, /nan))])[0]return, interp_value
end;+
;   函数用途:
;       该函数用于对栅格矩阵中缺失值基于滑动窗口进行填充
;   函数参数:
;       dataset: 需要进行填补的栅格矩阵
;       dataset_interp: 输出的经填补好的栅格矩阵
;       window_size(默认: 3): 窗口大小(奇数)
;       interp: 填充的方法(1: 窗口均值; 2: 最近邻)
;-
pro window_interp, dataset, dataset_interp, window_size = window_size, interp = interpds_size = size(dataset, /dimensions)ds_type = size(dataset, /type); 边界填充if ~keyword_set(window_size) then window_size = 3padding_size = window_size / 2ds_size += padding_size * 2dataset_pad = padding(dataset, padding_size)dataset_interp = padding(dataset, padding_size); 插值for col_ix=padding_size, ds_size[0] - padding_size - 1 do beginfor row_ix=padding_size, ds_size[1] - padding_size - 1 do begin; 若不是NAN跳过if ~finite(dataset_pad[col_ix, row_ix], /nan) then continue; 取窗口数组window_ds = dataset_pad[col_ix-padding_size: col_ix+padding_size, $row_ix-padding_size: row_ix+padding_size]if (where(~finite(window_ds, /nan), /null) eq !null) then continue  ; 若窗口内均为NAN则跳过; 插值if interp eq 1 then interp_value = mean(window_ds, /nan)if interp eq 2 then interp_value = interp_nearest(window_ds); 赋值dataset_interp[col_ix, row_ix] = interp_valueendforendfor; no paddingdataset_interp = dataset_interp[padding_size:(ds_size[0] - padding_size - 1), $padding_size:(ds_size[1] - padding_size - 1)]
end

还有write_img熬，自带的write_tiff每次都得自己写地理结构体，也稍微封装了一下。

;+
;   函数用途:
;       用于输出tiff文件(封装write_tiff)
;   函数参数:
;       img_path: tiff文件的输出路径
;       img: 栅格矩阵
;       out_res: 输出分辨率
;       ul_x: 左上角格点的左上角位置的X坐标
;       ul_y: 左上角格点的左上角位置的Y坐标
;-
pro write_img, img_path, img, out_res, ul_x, ul_y; 地理结构体geo_info={$MODELPIXELSCALETAG: [out_res, out_res, 0.0], $  ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, ul_x, ul_y, 0.0], $  ; 角点信息GTMODELTYPEGEOKEY: 2, $  ; 设置为地理坐标系GTRASTERTYPEGEOKEY: 1, $  ; 像素的表示类型, 北上图像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $  ; 地理坐标系为WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102}  ; 单位为度; 输出write_tiff, img_path, img, geotiff=geo_info, /float
end

时间精力有限，不再详细说明，Bye~.