pro read_hdf

dir='/misc/mcc31/common/satellite/modis/land_cover/'
fn='MCD12C1.A2023001.061.2024251212901.hdf'
print,fn

lon=findgen(7200)*0.05-180.+0.025
lat=reverse(findgen(3600)*0.05- 90.+0.025)
nlon=n_elements(lon)
nlat=n_elements(lat)
print,nlon,nlat

; Open hdf file
; Frlm hdfdump: There are three landcover datasets, IGBP, UMD, LAI_FPAR,
; each has different landcover classes. See their respective index of classes
; below. Name of each class is listed in the hdf file.

sd_id=hdf_sd_start(dir+fn,/read)
sds_name='Majority_Land_Cover_Type_1'
sds_index=hdf_sd_nametoindex(sd_id,sds_name)
sds_id=hdf_sd_select(sd_id,sds_index)
hdf_sd_getinfo,sds_id,name=sds_name,natts=num_attributes,ndim=num_dims, $
dims=dim_sizes
hdf_sd_getdata,sds_id,landc1
print,max(landc1),min(landc1)
class1=[indgen(8),indgen(8)+10,20]
nclass1=n_elements(class1)
data1='IGBP'

sds_name='Majority_Land_Cover_Type_2'
sds_index=hdf_sd_nametoindex(sd_id,sds_name)
sds_id=hdf_sd_select(sd_id,sds_index)
;hdf_sd_getinfo,sds_id,name=sds_name,natts=num_attributes,ndim=num_dims, $
;dims=dim_sizes
hdf_sd_getdata,sds_id,landc2
print,max(landc2),min(landc2)
class2=[indgen(8),indgen(3)+10,14,15,17]
nclass2=n_elements(class2)
data2='UMD'

sds_name='Majority_Land_Cover_Type_3'
sds_index=hdf_sd_nametoindex(sd_id,sds_name)
sds_id=hdf_sd_select(sd_id,sds_index)
;hdf_sd_getinfo,sds_id,name=sds_name,natts=num_attributes,ndim=num_dims, $
;dims=dim_sizes
hdf_sd_getdata,sds_id,landc3
print,max(landc3),min(landc3)
class3=[indgen(8),indgen(3)+10]
nclass3=n_elements(class3)
data3='IAI_FPAR'

; Below is what I did for a quick look at the global map for selected dataset
;------------------------------------------------------
ict=39
fign='Landcover_'+data1
lvl=class1
clr=indgen(nclass1)*15+1 & clr(0)=255
pmapgif,landc1,lon,lat,levels=lvl,colortable=ict,colorindex=clr,region='glb',$
figname=fign,title='Landcover '+data1
print,fign+'.gif'
;------------------------------------------------------

stop
end