前段時(shí)間需要用到SMAP L4級(jí)的土壤濕度產(chǎn)品，但是它是H5格式的并且沒(méi)有地理坐標(biāo)，找了好多資料都沒(méi)有找到合適的直接可用的教程，在師兄的幫助下經(jīng)過(guò)了幾天的嘗試，最終可算是成功轉(zhuǎn)成了GeoTiff，不知道是不是全網(wǎng)首例hhhhh，在女朋友的鼓勵(lì)之下決定寫(xiě)下人生中第一篇技術(shù)貼分享下。

文章目錄

• • 探索
  • 實(shí)現(xiàn)

---

探索

我找了很久，對(duì)于我來(lái)說(shuō)比較有參考價(jià)值的只有：SMOS、AMSR２以及SMAP三種土壤水分遙感產(chǎn)品的下載和預(yù)處理（ps：這篇文章中提到的兩種方法都不適用）、How to import and geolocate SMAP Level-3 and Level-4 data in ENVI和Perform Time-Series Analysis of soil moisture from SMAP L3 Product，他們處理的都是L3級(jí)的數(shù)據(jù)，但是對(duì)于我的啟發(fā)還是很大的，感謝。

接下來(lái)，我發(fā)現(xiàn)：

SMAP L4產(chǎn)品使用的是EASE-Grid 2.0投影，需要從ftp://sidads.colorado.edu/pub/tools/easegrid2/下載與SMAP數(shù)據(jù)分辨率相對(duì)應(yīng)的格網(wǎng)；大家也可以百度云自取提取碼：e1oh；

上述幾篇教程提到的bulid GLT參數(shù)或者是格網(wǎng)頭文件對(duì)于L4產(chǎn)品來(lái)說(shuō)都沒(méi)有可用的。

既然用手頭現(xiàn)成的軟件不行，就只能嘗試用Python來(lái)試試了~

實(shí)現(xiàn)

在實(shí)現(xiàn)的過(guò)程中我發(fā)現(xiàn)數(shù)據(jù)是全球范圍的，EASE-Grid 2.0格網(wǎng)文件中每?jī)蓚€(gè)坐標(biāo)間隔是不一樣的，但我的研究區(qū)范圍很小，十幾個(gè)像元范圍內(nèi)的影響應(yīng)該也很小，所以我將格網(wǎng)、數(shù)據(jù)都按照我的研究范圍裁剪了，并將中間像元與上一個(gè)像元的坐標(biāo)差作為最終輸出的GeoTiff的經(jīng)緯度分辨率。

import h5py import numpy as np import os import glob from osgeo import gdal,osr import sys import math# get list of all files work_dir = ###your working directory### os.chdir(work_dir)flist = glob.glob('*.h5')# Read binary files and reshape to correct size lats = np.fromfile('EASE2_M09km.lats.3856x1624x1.double', dtype=np.float64).reshape((1624,3856))#< reshape to dimensions above lons = np.fromfile('EASE2_M09km.lons.3856x1624x1.double', dtype=np.float64).reshape((1624,3856))# Set the ROI extent sel_col_start = 3037 sel_row_start = 437 sel_col_end = 3050 sel_row_end = 449num_row = sel_row_end - sel_row_start num_col = sel_col_end - sel_col_start if num_row >= 2:mid_row = math.ceil(num_row/2) if num_col >= 2:mid_col = math.ceil(num_col/2)# select the coordinates of ROI lats = lats[ sel_row_start:sel_row_end, : ][ :, sel_col_start:sel_col_end ] lons = lons[ sel_row_start:sel_row_end, : ][ :, sel_col_start:sel_col_end ]# Calculate the upperleft corner coordinates with the central cell size xlons_per_cell = abs(lons[mid_row][mid_col] - lons[mid_row-1][mid_col-1]) ylats_per_cell = abs(lats[mid_row][mid_col] - lats[mid_row-1][mid_col-1])group_id = 'Geophysical_Data' var_id = [ 'sm_rootzone', 'sm_rootzone_pctl', 'sm_rootzone_wetness' ] for i in range(len(flist)):try:with h5py.File(flist[i], 'r') as f:for j in range(len(var_id)):sm_data = f[group_id][var_id[j]][:,:]sm_data = sm_data.astype(np.float32).reshape((1624,3856))sm_data[sm_data==-9999.000000] = np.nan# select the ranks of ROIsm_data = sm_data[ sel_row_start:sel_row_end, : ][ :, sel_col_start:sel_col_end ] #Write GeoTiff with GDALdriver = gdal.GetDriverByName('GTiff')dataset = driver.Create( var_id[j] + '_' + os.path.splitext(flist[i])[0] + '.tif', num_col, num_row, 1, gdal.GDT_Float32 )dataset.SetGeoTransform([lons[0][0], xlons_per_cell, 0, lats[0][0], 0, -ylats_per_cell])sr = osr.SpatialReference()sr.SetWellKnownGeogCS('WGS84')dataset.SetProjection(sr.ExportToWkt()) dataset.GetRasterBand(1).WriteArray(sm_data)del datasetexcept OSError:print('ERROR:Cannot open file, please check file ', flist[i], '!')with open('Badfiles.txt', 'a') as bf:bf.write(str(flist[i]) + '\n')continue

總結(jié)

以上是生活随笔為你收集整理的SMAP L4级土壤湿度产品的预处理的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。