;VIIRS SDR simplest reader
;Purpose: Use IDL to read and display a VIIRS SDR granule in HDF5 format with minimal coding.  Display the map on screen.
;Author: Changyong Cao, NOAA/NESDIS/STAR/NCC
;Date: December 1, 2011
;      Feb. 9, 2012: replaced the dialate with the missing values approach to handle the bow-tie effect better.  Thanks to Hongqing Liu for suggestions.

 


PRO ReadVIIRSSDR3
;file_name='D:\Home\ccao\Work2011\VIIRS\VIIRS1STImage\SVM02_npp_d20111121_t1610296_e1611538_b00344_c20111121224330078798_noaa_ops.h5'
file_name = 'D:\Home\ccao\Work2011\VIIRS\VIIRS_anomaly\GMODO-SVM01_npp_d20111125_t1631310_e1637124_b00401_c20111126184122718562_noaa_ops.h5'
file_id = H5F_OPEN(file_name)

;get one band image
VIIRSBAND=H5G_GET_MEMBER_NAME(file_id,'All_Data',0)
Data='/All_Data/'+VIIRSBAND+'/Radiance'
print,data
dataset_id= H5D_OPEN(file_id,Data)
img = H5D_Read(dataset_id)
H5D_CLOSE, dataset_id

;get latitude
Geo=H5G_GET_MEMBER_NAME(file_id,'All_Data',1)
Data='/All_Data/'+Geo+'/Latitude'
dataset_id= H5D_OPEN(file_id,Data)
latBuf = H5D_Read(dataset_id)
H5D_CLOSE, dataset_id

;get longitude
Geo=H5G_GET_MEMBER_NAME(file_id,'All_Data',1)
Data='/All_Data/'+Geo+'/Longitude'
dataset_id= H5D_OPEN(file_id,Data)
lonBuf = H5D_Read(dataset_id)
H5D_CLOSE, dataset_id
H5F_CLOSE, file_id

img1=img(0:2000,0:2000)
lat1=latbuf(0:2000,0:2000)
lon1=lonbuf(0:2000,0:2000)

window,0,xsize=1000,ysize=1000;,/PIXMAP
idx=where (img1(*,*) eq 65533)
img1(idx)=!Values.F_NAN
img1 = BYTSCL(img1,TOP=253,/NAN)

latmin = MIN( Lat1, MAX = latmax )
lonmin = MIN( Lon1, MAX = lonmax )

;Hope that the map_set routine will determine the aspect ratio correctly!
ysz=1000 ;note: it appears that the largest PIXMAP I can create is 2000x2000!
xsz=Ysz
print, 'Image size: ',Xsz,' columns, and ',Ysz,' rows. '
;WINDOW, 0, XSize = Xsz, Ysize=Ysz, Title = Swathlist + " : " + FieldUsed,/PIXMAP;note: pixmap makes invisible window
;note: Plat has to be 0 in order to get a geographic projection!!! Wish they documented this!!!
;/GRID didn't do anything

targetlat=0.0
targetlon=-45.
MAP_SET, targetlat,targetlon, /Satellite,limit = [ latmin, lonmin, latmax, lonmax ],/isotropic;,/GRID
img1(idx)=255
img1=MAP_PATCH(img1, Lon1, Lat1,MAX_VALUE=254,XSTART=x0, YSTART=y0)

TV,img1,X0,Y0
MAP_CONTINENTS, /LoRes
;MAP_GRID, /Box_Axes
END