2

我一直在使用 SMAP 数据卫星,专门用于湿度和土壤比例。

我遵循使用 GDAL 解决所有问题的想法,并制作了类似于Link to first approach to download SMAP data

修改代码和测试:

import os
import h5py
import numpy as np
from osgeo import gdal, gdal_array, osr


  # the file to download 

https://n5eil01u.ecs.nsidc.org/SMAP/SPL4SMAU.003/2017.08.01/SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5

  path = "/path/to/data"
  h5File = h5py.File(path + "SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5", 'r')
 data = h5File.get('Analysis_Data/sm_rootzone_analysis')
 lat = h5File.get("cell_lat")
 lon = h5File.get("cell_lon")

 np_data = np.array(data)
 np_lat = np.array(lat)
 np_lon = np.array(lon)

 num_cols = float(np_data.shape[1])
 num_rows = float(np_data.shape[0])

xmin = np_lon.min()
xmax = np_lon.max()
ymin = np_lat.min()
ymax = np_lat.max()
xres = (xmax - xmin) / num_cols
yres = (ymax - ymin) / num_rows

nrows, ncols = np_data.shape
xres = (xmax - xmin) / float(ncols)
yres = (ymax - ymin) / float(nrows)
geotransform = (xmin, xres, 0, ymax, 0, -xres)

dataFileOutput = path + "sm_rootzone_analysis.tif"
output_raster = gdal.GetDriverByName('GTiff').Create(dataFileOutput, ncols, nrows, 1, gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  
srs = osr.SpatialReference() 
srs.ImportFromEPSG(4326)  

output_raster.SetProjection(srs.ExportToWkt()) 
output_raster.GetRasterBand(1).WriteArray(np_data)  # Writes my array to the raster

del output_raster

因此,使用这种方法,结果是具有许多投影问题的全局地图,例如下面的图像,由上面的 python 代码生成。 使用上述代码生成的地图

为了与正确的数据进行比较,使用 HEG nasa 软件从 h5 中提取了相同的图像。 地图正确

4

1 回答 1

1

如果数据确实在 EASE2 全局网格中,则不应将 EPSG:4326 指定为地理变换中具有纬度/经度的坐标系。

如果您将纬度/经度坐标转换为 9 公里处的 EASE2 网格,您的地理变换应该类似于:

geotransform = (-17367530.44516138, 9000, 0, 7314540.79258289, 0, -9000.0)

和 srs:

srs.ImportFromEPSG(6933)

于 2017-09-01T12:00:13.810 回答