SMAP L4级土壤湿度产品的预处理

前段时间需要用到SMAP L4级的土壤湿度产品，但是它是H5格式的并且没有地理坐标，找了好多资料都没有找到合适的直接可用的教程，在师兄的帮助下经过了几天的尝试，最终可算是成功转成了GeoTiff，不知道是不是全网首例hhhhh，在女朋友的鼓励之下决定写下人生中第一篇技术贴分享下。

---

探索

我找了很久，对于我来说比较有参考价值的只有：SMOS、AMSR２以及SMAP三种土壤水分遥感产品的下载和预处理（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级的数据，但是对于我的启发还是很大的，感谢。

接下来，我发现：

1. SMAP L4产品使用的是EASE-Grid 2.0投影，需要从ftp://sidads.colorado.edu/pub/tools/easegrid2/下载与SMAP数据分辨率相对应的格网；大家也可以百度云自取提取码：e1oh；
2. 上述几篇教程提到的bulid GLT参数或者是格网头文件对于L4产品来说都没有可用的。

既然用手头现成的软件不行，就只能尝试用Python来试试了~

实现

在实现的过程中我发现数据是全球范围的，EASE-Grid 2.0格网文件中每两个坐标间隔是不一样的，但我的研究区范围很小，十几个像元范围内的影响应该也很小，所以我将格网、数据都按照我的研究范围裁剪了，并将中间像元与上一个像元的坐标差作为最终输出的GeoTiff的经纬度分辨率。

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 = 449

num_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 ROI
				sm_data = sm_data[ sel_row_start:sel_row_end, : ][ :, sel_col_start:sel_col_end ]			
				
				#Write GeoTiff with GDAL
				driver = 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 dataset
	except 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