How to fix the reprojection from EASE-2 grid product SMAP to geographic coordinates?

Question

I'm have been working with SMAP data satellite, specially for moisture and soil proporties.

I follow the idea of use GDAL solve everything, and make something similar to this published in Link to first approach to download SMAP data

Modifing the code and testing:

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

So, using this approach, the result is a global map with many problems of projections, as for example the image below, produced by the python code above.

To compare with a correct data, the same image was extract from h5, using HEG nasa software.

Answer 1

If the data is really in the EASE2 Global grid, you shouldn't assign EPSG:4326 as a the coordinate system with lat/lon degrees in the geotransform.

If you convert the lat/lon coordinates to the EASE2 Grid at 9km, your geotransform should be something like:

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

and the srs:

srs.ImportFromEPSG(6933)