1. 内容
——————————————————————————————————————————
——————————————————————————————————————————
這裡主要沒有用到什麼新的知識點,就是先擷取資料,然後找到經緯度的最大值最小值,以此得到行列數去建構裝資料的數組,然後資料按順序放入到數組中。另外,如果你的資料源明确說明資料的分辨率是準确的而不是近似,那麼你可以不進行均值插值,當然為了確定你可以先出圖看看缺失值是否存在再來判斷是否需要均值插值。做好均值插值後最後輸出即可。
——————————————————————————————————————————————————————————————————————————————————————
2. 程式設計
pro geo_out_tiff, lon, lat, geo_data, resolution, geo_data_out, geo_info
; 本程式用于将經緯度、xx資料集轉化為可輸出為tiff的資料
; 确定最大最小的經緯度
lon_min = min(lon)
lon_max = max(lon)
lat_min = min(lat)
lat_max = max(lat)
; 擷取轉化為tiff資料所需要的行列數
column = ceil((lon_max - lon_min) / resolution) ; ceil()函數表示向上取整
row = ceil((lat_max - lat_min) / resolution)
; 得到每一個點的行列位置
column_array = floor((lon - lon_min) / resolution) ; 向下取整
row_array = floor((lat_max - lat) / resolution)
; 建立裝資料的box
data_box = fltarr(column, row)
; 将資料轉入到box裡面
data_box[column_array, row_array] = geo_data
; 用來裝填充之後的值的box
geo_data_out = fltarr(column, row)
; 但是這個樣子還是不行,由于我們這個txt檔案的資料源表明了這個分辨率隻是大緻的分辨率,是以有一些像元位置的data是沒有的,是以需要插值
; 下面這種方法沒有處理最外邊的行列
; for column_i = 1, column - 2 do begin
; for row_i = 1, row - 2 do begin
; if data_box[column_i, row_i] eq 0.0 then begin
; ; 建立九宮格移動視窗
; temp_windows = data_box[column_i - 1:column_i + 1, row_i - 1:row_i + 1]
; temp_windows = (temp_windows gt 0.0) * temp_windows
; if total(temp_windows gt 0.0) gt 3 then begin ; 門檻值設定
; geo_data_out[column_i, row_i] = total(temp_windows) / total(temp_windows gt 0.0)
; endif
; endif else begin
; geo_data_out[column_i, row_i] = data_box[column_i, row_i]
; endelse
; endfor
; endfor
; 現在寫的方法會處理最外邊的行列
geo_data_out_copy = fltarr(column + 2, row + 2)
geo_data_out_copy[1:column, 1:row] = data_box
for column_i = 1, column do begin
for row_i = 1, row do begin
if geo_data_out_copy[column_i, row_i] eq 0.0 then begin
; 建立九宮格移動視窗
temp_windows = geo_data_out_copy[column_i - 1: column_i + 1, row_i - 1:row_i + 1]
temp_windows = (temp_windows gt 0.0) * temp_windows
if total(temp_windows gt 0.0) gt 3 then begin
geo_data_out[column_i - 1, row_i - 1] = total(temp_windows) / total(temp_windows gt 0.0)
endif
endif else begin
geo_data_out[column_i - 1, row_i - 1] = geo_data_out_copy[column_i, row_i]
endelse
endfor
endfor
; 編寫地理資訊geotiff
geo_info={$
MODELPIXELSCALETAG:[resolution,resolution,0.0],$ ; 這裡需要填分辨率
MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$ ; 這裡需要填寫左上角點的經緯度,然後下面的都是預設不用修改
GTMODELTYPEGEOKEY:2,$
GTRASTERTYPEGEOKEY:1,$
GEOGRAPHICTYPEGEOKEY:4326,$
GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
GEOGANGULARUNITSGEOKEY:9102,$
GEOGSEMIMAJORAXISGEOKEY:6378137.0,$
GEOGINVFLATTENINGGEOKEY:298.25722}
end
pro week_six_study2
; 本程式(主程式)是将txt文本資料轉化為geotiff圖像輸出,并且由于有一些缺失值是以需要插值
; 路徑
in_path = 'D:/IDL_program/experiment_data/chapter_4/2013_year_aop.txt'
out_path = 'D:/IDL_program/experiment_data/chapter_4/'
; 打開txt文本檔案
openr, 1, in_path
; 建立存儲數組
geo_data = fltarr(5, 56304) ; 行列數通過記事本檢視得到
; 初始化值
geo_data[*, *] = -9999.0
; 擷取第一行的索引資料
index = ''
readf, 1, index
index = strsplit(index, /extract)
index_num = n_elements(index)
; 擷取txt檔案的資料
readf, 1, geo_data
; 已經将txt檔案的資料全部放在了變量geo_data裡面,是以txt檔案可以關閉了
free_lun, 1
; 擷取經度資料
lon = geo_data[0, *]
; 擷取緯度資料
lat = geo_data[1, *]
; 已知該txt文本資料的分辨率是0.18
resolution = 0.18
; 擷取後面的資料
for data_i = 2, index_num - 1 do begin
; 輸出
geo_out_tiff, lon, lat, geo_data[data_i, *], resolution, geo_data_out, geo_info
write_tiff, out_path + index[data_i] + '.tiff', geo_data_out, geotiff=geo_info, /float
endfor
end
運作結果:
![](https://img.laitimes.com/img/9ZDMuAjOiMmIsIjOiQnIsIyZuBnLhNTOjJmN5UjZhJmYlFGZiNTOxQjZ5kTZzMmMiJzYiRzLc52YucWbp5GZzNmLn9Gbi1yZtl2Lc9CX6MHc0RHaiojIsJye.png)
部分輸出結果展示:
這是未經過均值插值的結果:
———————————————————————————————————————————
這是經過了均值插值的結果:
——————————————————————————————————————————————————————————————————————————————————————
我是炒茄子,謝謝大家!
——————————————————————————————————————————————————————————————————————————————————————