2

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. 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

Iron Banker Of Braavos
  • 353
  • 5
  • 18

1 Answers1

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)

Rutger Kassies
  • 61,630
  • 17
  • 112
  • 97
  • 1
    Can you explain further why `6933`? I made a search and found epsg:3410 seems for the data http://spatialreference.org/ref/epsg/3410/ though I haven't tried it yet. – pateheo Oct 24 '18 at 11:09
  • 1
    That's the EPSG of the original EASE grid, not EASE-2 as referenced in this question. See for example the remark part at: https://epsg.io/3410 – Rutger Kassies Nov 04 '18 at 12:09
  • 2
    @RutgerKassies Your corner coordinates, geotransform[0] and geotransform[3], from above, don't match those defined by the EASE-Grid 2.0 authors. See: https://nsidc.org/ease/ease-grid-projection-gt Your pixel sizes are also slightly different. Any idea why? – Arthur Aug 08 '19 at 21:28
  • I don't recall exactly, I probably converted the lat/lon corner coordinates from the metadata in the file to EASE-2. The old EASE website didn't provide the corner coordinates if i recall correctly, but now it does. It's certainly best to stick to what NSIDC is saying themselves. Thanks for the link! – Rutger Kassies Aug 09 '19 at 06:59