%======================================================
%  Quick Display tool for HDF file
%                                 H.Murakami 2002.11.15
%=========================================== disp_hdf.m
clear
%
binout='off';
%
filename=input('> Input (path/)File = ','s');
%
%--Search for HDF data sets
fileID = hdfh('open',filename,'read',0);
if fileID<0;
 disp(['No input data = ',filename]);
 return;
end
anID = hdfan('start',fileID);
%
ga_on=input('> Show global attributes? [y or n] = ','s');
if size(ga_on,2)<1;ga_on='n';end
%
%--Scientific Data Set
sdID = hdfsd('start',filename,'read');
[ndatasets,nglobal_attr,status] = hdfsd('fileinfo',sdID);
%
if ga_on(1:1)=='y';
 for i=0:nglobal_attr-1
  [name,data_type,count,status] = hdfsd('attrinfo',sdID,i);
  if data_type(1:4)=='char';
   [data,status] = hdfsd('readattr',sdID,i);
   ll=min(42,count);
   disp(['GA_idx=',int2str(i),', GA= ',name,' = ',data(1:ll)])
  elseif count<3
   [data,status] = hdfsd('readattr',sdID,i);
   gclass=class(data);
   if gclass(1:4)=='doub';
    disp(['GA_idx=',int2str(i),', GA= ',name,' = ',num2str(data)])
   else
    disp(['GA_idx=',int2str(i),', GA= ',name,' = ',num2str(double(data))])
   end
  else
   disp(['GA_idx=',int2str(i),', GA= ',name,' --'])
  end
 end
end
%
dataname=input('> SDS Name = ','s');
if size(dataname,2)>0;
 index = hdfsd('nametoindex',sdID,dataname);
 if index<0;index=str2num(dataname);
  if size(index,2)<1;index=-1;end
 end
else
 index=-1;
end
if index<0;
 for i=0:ndatasets-1
  sdsID = hdfsd('select',sdID,i);
  if sdsID>-1;
   [sdsName,rank,dimSizes,sddataType,nattrs,status] = hdfsd('getinfo',sdsID);
   disp(['index=',int2str(i),', SDS= ',sdsName])
  end
 end
 index=input('> SDS index = ');
 if index<0;
  dataname=input('> SDS Name = ','s');
  index = hdfsd('nametoindex',sdID,dataname);
 end
end
sdsID = hdfsd('select',sdID,index);
%ref = hdfsd('idtoref',sdsID);
[sdsName,rank,dimSizes,sddataType,nattrs,status] = hdfsd('getinfo',sdsID);
%
%--Get SD attribute information. The index for readattr is zero based.
unit='DN';
if nattrs>0
 for i = 1:nattrs
  [arrayAttribute(i).Name,attrDataType,count,status] = hdfsd('attrinfo',sdsID,i-1);
  [arrayAttribute(i).Value, status] = hdfsd('readattr',sdsID,i-1);
  if attrDataType(1:3)=='cha';
   disp([arrayAttribute(i).Name,' = ',arrayAttribute(i).Value])
   if arrayAttribute(i).Name(1:3)=='uni';unit=arrayAttribute(i).Value;end
  else
   if mean(size(arrayAttribute(i).Value))<9;
    value1=double(arrayAttribute(i).Value);
    disp([arrayAttribute(i).Name,' = ',num2str(value1)])
   else;disp([arrayAttribute(i).Name])
   end
  end
 end
else
 arrayAttribute = [];
 sdDataType = '';
end
%
IsScale = logical(hdfsd('iscoordvar',sdsID));
if IsScale == 0
 for i=1:rank
   dimID = hdfsd('getdimid',sdsID,i-1);
   %Use sizes from SDgetinfo because this size may be Inf
   [dimName{i}, sizeDim,DataType{i}, nattrs, status] = hdfsd('diminfo',dimID);
   if strcmp(DataType{i},'none')
     Scale{i} = 'none';
   else
     [Scale{i},status] = hdfsd('getdimscale',dimID);
   end
   Size{i} = dimSizes(i);
   if nattrs>0
     for j=1:nattrs
       [Name{j},dataType,count,status] = hdfsd('attrinfo',dimID,j-1);
       [Value{j}, status] = hdfsd('readattr',dimID,j-1);
     end
     Attributes{i} = struct('Name',Name(:),'Value',Value(:));
   else
     Attributes = [];
   end
 end
 dims = struct('Name',dimName(:),'DataType',DataType(:),'Size',Size(:),'Scale',Scale(:),'Attributes',Attributes(:));
else
 dims = [];
end
%
%--Read SDS
start=zeros(rank,1);
stride=ones(rank,1);
edge=dimSizes(1:rank);
nb=1;
if rank==3
  disp(['number of band = ',int2str(dimSizes(1))])
  nb=input('> band = ');
  start(1)=nb-1;stride(1)=1;edge(1)=1;
end
if size(sdsName,2)>12;
 if sdsName(1:7)=='l1b_ch1';
   nb0=input('> band [0(rgb) or 1]= ');
   if size(nb0,1)<1;nb=1;else;nb=nb0;end
 end
end
if nb>0;
 [data,status] = hdfsd('readdata',sdsID,start,stride,edge);
elseif sdsName(1:6)=='l1b_ch'
 [data,status] = hdfsd('readdata',sdsID,start,stride,edge);
 status = hdfsd('endaccess',sdsID);
 index = hdfsd('nametoindex',sdID,'l1b_ch8_data');  %for GLI
% index = hdfsd('nametoindex',sdID,'l1b_ch9_data');  %for MODIS CNV
 sdsID = hdfsd('select',sdID,index);
 [data2,status] = hdfsd('readdata',sdsID,start,stride,edge);
 status = hdfsd('endaccess',sdsID);
 index = hdfsd('nametoindex',sdID,'l1b_ch5_data');
 sdsID = hdfsd('select',sdID,index);
 [data3,status] = hdfsd('readdata',sdsID,start,stride,edge);
else
 start(1)=6-1;stride(1)=1;edge(1)=1;
 [data,status] = hdfsd('readdata',sdsID,start,stride,edge);
 start(1)=5-1;stride(1)=1;edge(1)=1;
 [data2,status] = hdfsd('readdata',sdsID,start,stride,edge);
 start(1)=3-1;stride(1)=1;edge(1)=1;
 [data3,status] = hdfsd('readdata',sdsID,start,stride,edge);
end
%
%--Close interfaces
status = hdfsd('endaccess',sdsID);
%
%--Populate output structure
sdinfo = struct('Filename',[],'Type',[],'Name',[],'Rank',[],'DataType',[],'Attributes',[],'Dims',[],'Label',[],'Description',[],'Index',[]);
sdinfo.Filename = filename;
sdinfo.Name = sdsName;
sdinfo.Index = index;
sdinfo.Rank = rank;
sdinfo.DataType = sddataType;
if ~isempty(arrayAttribute)
  sdinfo.Attributes = arrayAttribute;
end
if ~isempty(dims)
  sdinfo.Dims = dims;
end
sdinfo.Type = 'Scientific Data Set';
%
slope=1;
if size(sdsName,2)>5 & sum(sdsName(1:6)=='l1b_ch')==6;
 if size(sdsName,2)<14;
  index = hdfsd('nametoindex',sdID,'gsys');
  if index<0;
   index = hdfsd('nametoindex',sdID,'c1');
  end
 else
  index = hdfsd('nametoindex',sdID,'gsys_2km');
 end
 sdsID = hdfsd('select',sdID,index);
 if sdsID>-1;
  [sdsName2,rank,dimSizes,sddataType,nattrs,status] = hdfsd('getinfo',sdsID);
  start=zeros(rank,1);
  stride=ones(rank,1);
  edge=dimSizes(1:rank);
  [gsys,status] = hdfsd('readdata',sdsID,start,stride,edge);
  status = hdfsd('endaccess',sdsID);
  if sdsName(8:8)=='_';gch=str2num(sdsName(7:7));
  else;gch=str2num(sdsName(7:8));
  end
  if gch>19;
   if     size(data,1)>1300;
    if gch<28;slope=gsys(gch-19);  %250m
    else     ;slope=gsys(gch-23);  %250m
    end
   elseif     gch>19 & gch<24;
    if gch<28;slope=gsys(gch-19);  %250m
    else     ;slope=gsys(gch-23);  %250m
    end
   elseif gch<28;slope=gsys(gch-23);
   elseif gch<30;slope=gsys(gch-27);
   elseif gch<37;slope=gsys(gch-29);
   else;slope=1;
   end
  else
   nbr=0;
   for nb1=1:gch-1
    nbr=nbr+1;
    if (nb1==4)|(nb1==5)|(nb1==7)|(nb1==8);
     nbr=nbr+1;
    end
   end
   nbr=nbr+1;
   slope=gsys(nbr);
  end
 end
end
 
%--Close annotation interface
status = hdfan('end',anID);
status = hdfh('close',fileID);
%
%
%---------------------- Display image ---------------------------
mll=size(data,2);nll=size(data,1);dclass='     ';
dclass=class(data);
filename1=filename;sdsName1=sdsName;
nlg=size(filename1,2);
for k=1:nlg;if filename1(k:k)=='_';filename1(k:k)='-';end;end
nlg=size(sdsName1,2);
for k=1:nlg;if sdsName1(k:k)=='_';sdsName1(k:k)='-';end;end
%
if mll==1 & nll==1;
  if dclass(1:4)=='doub';
   disp(['SDS:',sdsName,' = ',num2str(data)])
  else
   disp(['SDS:',sdsName,' = ',num2str(double(data))])
 end
elseif dclass(1:4)=='char' & (mll<2 |  nll<2);
 if     mll==1;disp(['SDS:',sdsName,' = ',data']);
 elseif nll==1;disp(['SDS:',sdsName,' = ',data]);
 end
elseif mll==1 | nll==1;
 clf
 axes('position',[0.1 0.1 0.8 0.8])
 ll=max(mll,nll);
 plot([1:ll],data,'-k')
 title(['SDS Name= ',sdsName1])
 grid on;
else
 a0=fopen('/export/emc2/murakami/Scat.pal','r');
 if a0>-1;ba1=fread(a0,[256 3],'uint8');fclose(a0);gg=ba1/255;
 else;gg=jet(256);
 end
%
clf
if nb==0;
 img=uint8(zeros(mll,nll,3));
%-- for GLI --
 if sdsName(1:6)=='l1b_ch';
  data =bitand(data,8191);
  data2=bitand(bitxor(data2,61440),8191);
  data3=bitand(bitxor(data3,61440),8191);
 end
%-- for GLI --
 img(:,:,1)=uint8(bitshift(data ,-2))';
 img(:,:,2)=uint8(bitshift(data2,-4))';
 img(:,:,3)=uint8(bitshift(data3,-4))';
%-- for GLI --
% img(:,:,1)=uint8(bitshift(data ,-2))';
% img(:,:,2)=uint8(bitshift(data2,-2))';
% img(:,:,3)=uint8(bitshift(data3,-3))';
%-- for MODISCNV --
% img(:,:,1)=uint8(bitshift(data ,-3))';
% img(:,:,2)=uint8(bitshift(data2,-3))';
% img(:,:,3)=uint8(bitshift(data3,-2))';
%------------------
 dmax1=max(max(data));
 dmax2=max(max(data2));
 dmax3=max(max(data3));
 dmax=255;
 dmin=0;
elseif sum(dclass(1:4)=='int8')==4 & nb>0;
 img=data';
 dmax=max(max(img));
 dmin=min(min(img));
elseif (sdsName(1:2)=='ch' | sdsName(1:2)=='cd' | (sdsName(1:1)=='s' & size(sdsName,2)==2));
 if sdsName(1:2)=='ch';slope=0.0015;else;slope=0.001;end
 img=double(data');
 msk=(img<1)+(img>65535);
 img=log10((img.*(1-msk)+msk)*slope);
 dmax=max(max(img.*(1-msk)));
 dmin=min(min(img));
 if sdsName(1:2)=='ch';dmin=-2;end
 dclass=class(data);
 if size(dclass,2)==5;
  if dclass(1:5)=='uint8';
   if sdsName(1:2)=='ch';slope=0.015;else;slope=0.01;end
   img=double(data');
   msk=(img<1)+(img>254);
   img=(img.*(1-msk)+msk)*slope-2;
   dmax=log10(100);dmin=log10(0.01);
  end
 end
 gg(1,:)=[1,1,1];gg(256,:)=[0.6,0.6,0.6];
elseif size(sdsName,2)>5 & sum(sdsName(1:6)=='l1b_ch')==6;
 disp('> This is GLI L1B data');
 img=bitand(data,8191)';
 dmax=double(max(max(img)));
 dmin=double(min(min(img)));
 dmin=max(0,dmin-(dmax-dmin)*0.01);
elseif size(sdsName,2)>3 & sum(sdsName(1:4)=='flag')==4;
 disp('>> This is GLI L2 flag data');
 img=uint16(log2(double(data)+1)*1000)';
 dmax=max(max(img));
 dmin=min(min(img));
elseif size(sdsName,2)>7 & sum(sdsName(1:8)=='l2_flags')==8;
 disp('>> This is GLI L2 l2_flags data');
 img=uint16(log2(double(data)+1)*1000)';
 dmax=max(max(img));
 dmin=min(min(img));
elseif dclass(1:4)=='char' & (mll>1 &  nll>1);
 img=uint8(data)';
 dmax=max(max(img));
 dmin=min(min(img));
elseif size(sdsName,2)>7 & sum(sdsName(1:8)=='Mean for')==8;
 img=data';
 dmax=max(max(img));
 dmin=min(min(img));
 dmin=double(dmin);
 if dmin==-9999;
  dmin=0;
  dmax=110;
  img=uint16(double(img)+double(img==-9999)*(9999+110));
  gg=[0:0.01:1;0:0.01:1;0.2:0.008:1]';gg(1,:)=[0,0,0.2];gg(101,:)=[0.5,0.5,0.5];colormap(gg)
 end
else
 img=data';
 dmax=max(max(img));
 dmin=min(min(img));
end
%
if mll>nll;
 set(gcf,'papertype','usletter','PaperOrientation','portrait','PaperPosition',[0.24706 0.25882 8 10.494],'Units','inches');
 scps=get(gcf,'position');
 set(gcf,'position',[scps(1)+scps(3)-7.2,scps(2)+scps(4)-9.4446,7.2,9.4446])
else
 set(gcf,'papertype','usletter','PaperOrientation','landscape','PaperPosition',[0.24706 0.25882 10.494 8],'Units','inches');
 scps=get(gcf,'position');
 set(gcf,'position',[scps(1)+scps(3)-9.4446,scps(2)+scps(4)-7.2,9.4446,7.2])
end
if mll>nll;nsz=min(0.76,0.80*nll/mll*1.3125);msz=nsz/1.3125*mll/nll;
else;nsz=min(0.76,0.80*nll/mll/1.3125);msz=nsz*1.3125*mll/nll;end
%
%----------------------------- Scale bar ------------------------------------
if nb>0;
 axes('position',[0.10+nsz 0.08 0.05 0.50])
 dmax=double(dmax);
 dmin=double(dmin);
 ddd=(dmax-dmin)/100;ba=[dmin:ddd:dmax];
 bb=mat2gray([ba;ba;ba;ba;ba;ba;ba;ba;ba;ba]',[dmin dmax]);
 imshow(bb);axis on;axis('xy')
 set(gca,'xtick',[ ]);
 if size(dclass,2)<5;
  rbar=[dmin:(dmax-dmin)/5:dmax];rbar0=[1:20:101];
  xlabel([unit],'fontsize',10);
  text( 5,110,['DN'],'fontsize',10);
  text(-2,105,['x ',num2str(double(slope))],'fontsize',10);
  set(gca,'ytick',rbar0,'yticklabel',rbar,'fontsize',10,'yaxislocation','right');
 elseif sdsName(1:2)=='ch';
  rbar=[0.01,0.03,0.1,0.3,1,3,10,30,90];rbar0=((log10(rbar)-dmin)/ddd)+1;
  xlabel('[mg m^{-3}]','fontsize',10);
  set(gca,'ytick',rbar0,'yticklabel',rbar,'fontsize',10,'yaxislocation','right');
 elseif dclass(4:5)~='t8' & sdsName(1:2)=='ss';
  rbar=[dmin:(dmax-dmin)/5:dmax];rbar0=[1:20:101];
  xlabel(unit,'fontsize',10);
  set(gca,'ytick',rbar0,'yticklabel',rbar,'fontsize',10,'yaxislocation','right');
 else;
  rbar=[dmin:(dmax-dmin)/5:dmax];rbar0=[1:20:101];
  xlabel([unit],'fontsize',10);
  text( 5,110,['DN'],'fontsize',10);
  text(-2,105,['x ',num2str(double(slope))],'fontsize',10);
  set(gca,'ytick',rbar0,'yticklabel',rbar,'fontsize',10,'yaxislocation','right');
 end
 grid on
end
%------------------------------------------------------------------------------
%
%------------------------image--------------------
h=axes('position',[0.08 0.08 nsz msz]);
imshow(img,[dmin dmax]);
axis on;axis('ij');axis('equal')
grid off;
%grid on;
%set(gca,'xtick',lo1,'xticklabel',lon,'fontsize',11);
%xlabel('Longitude','fontsize',13);
%set(gca,'ytick',la1,'yticklabel',lat,'fontsize',11);
%ylabel('Latitude','fontsize',13);
title([filename1,', ',sdsName1],'fontsize',13)
%
%---------------------------------------------------
colormap(gg)
end
%
if binout(1:2)=='on';
 f1=fopen('dump.dat','w');
 if     dclass(1:2)=='ch';fwrite(f1,data,'char');
 elseif dclass(3:4)=='t8';fwrite(f1,data,'int8');
 elseif dclass(3:4)=='t1';fwrite(f1,data,'int16');
 elseif dclass(4:5)=='t8';fwrite(f1,data,'uint8');
 elseif dclass(4:5)=='t1';fwrite(f1,data,'uint16');
 elseif dclass(1:2)=='lo';fwrite(f1,data,'int32');
 elseif dclass(1:2)=='si';fwrite(f1,data,'single');
 elseif dclass(1:2)=='do';fwrite(f1,data,'double');
 end
 fclose(f1);
end
