Thursday 22 November 2018

coordinate system - gdal for reprojection of VSCMO VIIRS data


I downloaded h5 VIIRS file by CLASS NOOA website.


In particular I am interested in Snow Cover Binary Maps EDR and I need to reproject them in lat/lon or UTM.


I tried to use gdal warp to this aim but, although I get correct latitude, I get completely wrong longitude (for example in my AOI I got 170 degrees instead of 10 degrees E), when I open tif file using qgis or envi


Using gdalinfo I get the names of the subdatasets. In particular one of them is the subdataset:


SUBDATASET_16_NAME=HDF5:"GITCO-VSCMO_npp_d20150102_t1051499_e1057303_b16487_c20160728125230331648_noaa_ops.h5"://All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap
SUBDATASET_16_DESC=[6144x6400]//All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap (8-bit unsigned character)


consequently the commmand line I tried is the following:


gdalwarp -t_srs 'epsg:4326' HDF5:"GITCO-VSCMO_npp_d20141212_t1043482_e1049286_b16189_c20160530182604182662_noaa_ops.h5"://All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap rip.tif 

I also tried to insert -s_srs 'epsg:4326' but the result is the same.


Do you have any suggestion, please?



Answer



Similar to Unable to warp HDF5 files, the hdf geolocation seems to be interpreted wrong by GDAL. If you extract the subdataset as a vrt file, it has GCP points inside, but they are above 180° longitude, and a shift of 180° does not fit. Subdatasets 2 and 3 contain latitudes and longitudes with reasonable values.


So I took the same workaround as given in the linked question, building a vrt file manually around the dataset:





0
1
0
1
1
HDF5:"gitco.h5"://All_Data/VIIRS-IMG-GEO-TC_All/Longitude
1
HDF5:"gitco.h5"://All_Data/VIIRS-IMG-GEO-TC_All/Latitude




HDF5:"gitco.h5"://All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap
1






and running gdalwarp without the use of the wrong GCP:



gdalwarp -geoloc -t_srs EPSG:4326 gitco.vrt gitco_16_warped.tif -overwrite -dstnodata 50

The result covers most of Europe:


enter image description here


I have applied a -dstnodata value of 50, because the data is mostly around 0 and 255.


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