3

I have been working with SMAP satellite data, especially for moisture and soil properties.

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 from   server

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, -yres)

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. Map produced using the code above

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

first suggestion

Following the modification suggested by @Ali:

The same code, with the modification using yres

The solution

Proposed by @AndreJ in Solution that solve the problem

reprojection working fine by gdal_translate

1

Following the solution proposed by @AndreJ proposed in this post I could fix.

My code:

gdal_translate -a_ullr -17367530.45 7314540.11 17367530.45 -7314540.11 -a_srs "+proj=cea +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m" -a_nodata -9999 HDF5:"SMAP_L4_SM_aup_20170801T030000_Vv3030_001.h5"://Analysis_Data/sm_rootzone_analysis sm-rootZone_analysis.tif

And the result was SMAP extraction vs costlines shape: Result of SMAP extraction vs costlines shape

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.