TransWikia.com

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

Geographic Information Systems Asked by Isaque Daniel on July 14, 2021

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

One Answer

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

Correct answer by Isaque Daniel on July 14, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP