01 前言

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

里面包含了眾多點位信息（不是站點數(shù)據(jù)），我們需要依據(jù)上述點的經(jīng)緯度信息放到對應(yīng)位置的像素點位置，放置完后如下：

可以發(fā)現(xiàn)，還存在部分缺失值，我們還需要進行缺失值的填補。

02 文本文件的讀取

IDL讀取文本文件還是不夠方便，稍微封裝了一下。

;+
;   函數(shù)用途:
;       用于讀取文本文件
;   函數(shù)參數(shù):
;       txt_path: 文本文件的路徑
;       ds: 讀取的輸出數(shù)據(jù)集(不含表頭)
;       header(關(guān)鍵字參數(shù)): 讀取的輸出表頭
;       separator(關(guān)鍵字參數(shù)): 分隔符,默認空白符 
;-
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; 讀取和處理數(shù)據(jù)集ds = strarr(file_lines(txt_path) - 1)readf, 1, dsds = list(ds, /extract)  ; 字符串數(shù)組轉(zhuǎn)化為列表ds = ds.map(lambda(e, separator: double(strsplit(e, separator, /extract))), separator)ds = ds.toarray()  ; 列表轉(zhuǎn)化為數(shù)組free_lun, 1
end

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

03 柵格矩陣的放置

于是乎，我們開始進行文件的讀取和轉(zhuǎn)數(shù)組。

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

在循環(huán)中，可以發(fā)現(xiàn)，使用了自定義的window_interp函數(shù)對target柵格矩陣進行缺失值的填補。

關(guān)于window_interp函數(shù)的定義由于封裝的比較多，疊的比較層數(shù)比較多，閱讀稍微有困難。學(xué)python的時候我是真的討厭那些一個簡單的功能的非要定義一個類，類又疊類，方法重組，來回找實現(xiàn)方法，來回折騰，本身功能不算特別復(fù)雜，但是被這么一折騰反而給閱讀代碼帶來困難。

然而，我現(xiàn)在還是成為了他們。But 我將盡量讓代碼的邏輯清晰可見，不過分抽象，如有必要我抽出部分功能進行整合，避免使用過于復(fù)雜的代碼框架搭建簡單的功能。

以下是關(guān)于填補缺失值的封裝函數(shù)，主要基于滑動窗口實現(xiàn)，包括滑動窗口均值填補和最近鄰填補。
涉及自定義函數(shù)：window_interp、padding、interp_nearest、meshgrid。

;+
;   函數(shù)用途:
;       用于對二維數(shù)組進行邊界填充
;   函數(shù)參數(shù):
;       array: 用于邊界填充的數(shù)組
;       pad_size: 單邊(上下左右)填充的大小
;       pad_value(關(guān)鍵字參數(shù): NAN): 填充的數(shù)值
;-
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;+
;   函數(shù)用途:
;       用于生成行列號格網(wǎng)矩陣
;   函數(shù)參數(shù):
;       cols_n: 列數(shù)
;       rows_n: 行數(shù)
;-
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; 獲取數(shù)組尺寸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;+
;   函數(shù)用途:
;       該函數(shù)用于對柵格矩陣中缺失值基于滑動窗口進行填充
;   函數(shù)參數(shù):
;       dataset: 需要進行填補的柵格矩陣
;       dataset_interp: 輸出的經(jīng)填補好的柵格矩陣
;       window_size(默認: 3): 窗口大小(奇數(shù))
;       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; 取窗口數(shù)組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  ; 若窗口內(nèi)均為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每次都得自己寫地理結(jié)構(gòu)體，也稍微封裝了一下。

;+
;   函數(shù)用途:
;       用于輸出tiff文件(封裝write_tiff)
;   函數(shù)參數(shù):
;       img_path: tiff文件的輸出路徑
;       img: 柵格矩陣
;       out_res: 輸出分辨率
;       ul_x: 左上角格點的左上角位置的X坐標
;       ul_y: 左上角格點的左上角位置的Y坐標
;-
pro write_img, img_path, img, out_res, ul_x, ul_y; 地理結(jié)構(gòu)體geo_info={$MODELPIXELSCALETAG: [out_res, out_res, 0.0], $  ; 分辨率MODELTIEPOINTTAG: [0.0, 0.0, 0.0, ul_x, ul_y, 0.0], $  ; 角點信息GTMODELTYPEGEOKEY: 2, $  ; 設(shè)置為地理坐標系GTRASTERTYPEGEOKEY: 1, $  ; 像素的表示類型, 北上圖像(North-Up)GEOGRAPHICTYPEGEOKEY: 4326, $  ; 地理坐標系為WGS84GEOGCITATIONGEOKEY: 'GCS_WGS_1984', $GEOGANGULARUNITSGEOKEY: 9102}  ; 單位為度; 輸出write_tiff, img_path, img, geotiff=geo_info, /float
end

時間精力有限，不再詳細說明，Bye~.