Как привязать репроекцию из продукта сетки SMAP EASE-2 к географическим координатам?

Question

Я работаю со спутником данных SMAP, специально для влажности и пропорций почвы.

Я следую идее использования GDAL, решаю все и делаю что-то похожее на опубликованное в Link to first подход для загрузки данных SMAP.

Модификация кода и тестирование:

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.

Чтобы сравнить с правильными данными, то же изображение было извлечено из h5 с использованием программного обеспечения HEG nasa.

Answer 1

Если данные действительно находятся в глобальной сетке EASE2, вам не следует назначать EPSG:4326 в качестве системы координат с градусами широты / долготы в геотрансформации.

Если вы преобразуете координаты широты / долготы в сетку EASE2 на 9 км, ваше геотрансформирование должно выглядеть примерно так:

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

и сестры:

srs.ImportFromEPSG(6933)