01 前言
此处的文本文件形式如下：
里面包含了众多点位信息（不是站点数据），我们需要依据上述点的经纬度信息放到对应位置的像素点位置，放置完后如下：
可以发现，还存在部分缺失值，我们还需要进行缺失值的填补。
02 文本文件的读取
IDL读取文本文件还是不够方便，稍微封装了一下。
;+
; 函数用途:
; 用于读取文本文件
; 函数参数:
; txt_path: 文本文件的路径
; ds: 读取的输出数据集(不含表头)
; header(关键字参数): 读取的输出表头
; separator(关键字参数): 分隔符,默认空白符
;-
pro read_txt, txt_path, ds, header=header, separator=separator
if ~keyword_set(separator) then separator = " "
; 读取和检索
openr, 1, txt_path ; 打开文本文件
; 是否指定输出的header
if ~arg_present(header) then begin
skip_lun, 1, 0
endif else begin
header = ''
readf, 1, header ; 读取表头
header = strsplit(header, separator, /extract)
endelse
; 读取和处理数据集
ds = strarr(file_lines(txt_path) - 1)
readf, 1, ds
ds = 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_dir
out_res = 0.18d
out_res_half = out_res / 2.0d
; 读取和检索
read_txt, in_path, ds, header=header
lon = ds[*, 0]
lat = ds[*, 1]
ds = ds[*, 2:*]
header = header[2:*]
lon_min = min(lon) - out_res_half
lon_max = max(lon) + out_res_half
lat_min = min(lat) - out_res_half
lat_max = max(lat) + out_res_half
cols = 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 begin
target = 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_max
endforeach
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_value
if ~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
; pad
pad_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)] = array
return, pad_array
end
;+
; 函数用途:
; 用于生成行列号格网矩阵
; 函数参数:
; cols_n: 列数
; rows_n: 行数
;-
function meshgrid, cols_n, rows_n
window_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)
end
function 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 / 2
center_row = rows_n / 2
distance = sqrt((cols - center_col) ^ 2.0 +(rows - center_row) ^ 2.0)
invalid_pos = where(finite(window_ds, /nan))
distance[invalid_pos] = !values.F_NAN
interp_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 = interp
ds_size = size(dataset, /dimensions)
ds_type = size(dataset, /type)
; 边界填充
if ~keyword_set(window_size) then window_size = 3
padding_size = window_size / 2
ds_size += padding_size * 2
dataset_pad = padding(dataset, padding_size)
dataset_interp = padding(dataset, padding_size)
; 插值
for col_ix=padding_size, ds_size[0] - padding_size - 1 do begin
for 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_value
endfor
endfor
; no padding
dataset_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, $ ; 地理坐标系为WGS84
GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $
GEOGANGULARUNITSGEOKEY: 9102} ; 单位为度
; 输出
write_tiff, img_path, img, geotiff=geo_info, /float
end
时间精力有限，不再详细说明，Bye~.
