Sunday 10 April 2016

raster - Opening data from NSIDC in QGIS?


i download data from here https://nsidc.org/data/SPL2SMP/versions/3#, but i can not open it in QGIS.


Do you know how to open this data in QGIS?



Answer



The raster can be loaded into QGIS with the HDF5 driver, but it gets displayed as X: 3 Y: 18570 cells. Gdalinfo reports just 2 subdatasets, but no longitudes or latitudes.



The information on the website notes that latitudes and longitudes should be stored in separate subdatasets, along with several other datasets.


Using H5dump from HDF Tools and HDF Compass, the missing tables appear, so this seems to be an unusual HDF5 storage that GDAL imports only partially. The current Windows binary of HDFview does not like HDF5 files at all.


Looking closer, you see that most datasets are one-dimensional arrays, with the EASE Grid row and cell number, latitude and longitude stored in separate datasets. Empty cells are simply discarded, saving much table space.


GDAL does not seem to like those 1-D tables. The landcover_class table that GDAL can read consists of 3 columns with the three dominant landcover classes per EASE grid cell. But without lon/lat information you can not use it.


I finally managed to get the SMAP L2 data tables with Panoply exported as longitude, latitude and landcover_class CSV files, merged them into one CSV table with Libre Office, and imported that as delimited text points into QGIS:


enter image description here


As you can see, this is a dataset from just one half-orbit.


The SMAP L3 dataset works far better. You can extract a subdataset with


gdal_translate -a_ullr -17334194 7356860 17334194 -7356860 -a_srs epsg:3410 -a_nodata -9999 HDF5:"smap_l3.h5"://Soil_Moisture_Retrieval_Data/soil_moisture sm.tif


to get a picture of all orbits of the day:


enter image description here


Unfortunately, this does not work for the landcover_class dataset, which comes as 3-dimensional (the three dominant land cover classes per EASE grid cell). GDAL tries to extract 406 bands of a 3x964 array, instead of 3 bands of a 964x406 array as intended by the data provider.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...