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级的数据，但是对于我的启发还是很大的，感谢。

接下来，我发现：

SMAP L4产品使用的是EASE-Grid 2.0投影，需要从ftp://sidads.colorado.edu/pub/tools/easegrid2/下载与SMAP数据分辨率相对应的格网；大家也可以百度云自取提取码：e1oh；

上述几篇教程提到的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 = 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

总结

以上是生活随笔为你收集整理的SMAP L4级土壤湿度产品的预处理的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。