%===== Display binary data of EORC/TRIC MODIS Near real time data =====%
%  Example;
%  >> dsp_modnrt
%  Input File = ...please type file name with path, and hit return...
%
%           MATLAB m-file  coded by H. Murakami, EORC/NASDA 4 June 2002
%======================================================================%
clear
infile=input('Input File = ','s');
%
f_print='off';     % <- PS-file output 'on' or 'off'
%
%--- lat lon range ---
lats=[-999,-999];lons=[-999,-999];  % if lats or lons ==-999;full scene
%------ if you want extract, please define lat-lon area like follows -----
%lats=[41,46];lons=[139,146];  % for zoom up of Hokkaido
%lats=[30,47];lons=[127,148];  % for zoom up of winter serge
%lats=[32,33.5];lons=[129.5,130.8];  % for zoom up of Ariake Bay
%lats=[25,35];lons=[124,132];  % for East China Sea
%lats=[43,48];lons=[142,146];  % for Okhotsk
%lats=[24,41];lons=[120,136];  % for East China Sea test
%lats=[43.8,44.8];lons=[143.2,144.6];  % for Saroma
%lats=[19,32];lons=[124,132];  % for 20020707
%lats=[34,37];lons=[138,141];  % for 20020711
%-------------------------------------------------------------------------
%
heda=zeros(1,7);
f1=fopen([infile],'r');
if f1>-1;
%
heda=fscanf(f1,'%d %d %g %g %g %g %g');fclose(f1);
f1=fopen([infile],'r');
fseek(f1,55,'bof');para0=fread(f1,[33,1],'char');fclose(f1);
para=char(para0(1:8)');
k1=0;k2=0;
for k=1:8;
 if k1==0 & para(k:k)~=' '; k1=k;end;
 if para(k:k)~=' '; k2=k;end;
end
para(1:k2-k1+1)=para(k1:k2);
if para(2:3)=='Kr';label=char(para0(10:33)');
else;              label=char(para0(10:24)');
end
%
nl=heda(1);ml=heda(2);dd=heda(5);
if lons(1)==-999;slon=heda(3);else;slon=lons(1);end
if lons(2)==-999;elon=slon+(nl-1)*dd;else;elon=lons(2);end
if lats(2)==-999;slat=heda(4);else;slat=lats(2);end
if lats(1)==-999;elat=slat-(ml-1)*dd;else;elat=lats(1);end
%
subdeg=[slon,elon,slat,elat];
%
subset0(1:2)=round((subdeg(1:2)-heda(3))/dd+1);
subset0(3:4)=round((heda(4)-subdeg(3:4))/dd+1);
nl0=subset0(2)-subset0(1)+1;ml0=subset0(4)-subset0(3)+1;
subset(1)=max( 1,subset0(1));if subset0(1)<1;ofsn=-subset0(1)-1;else;ofsn=0;end
subset(2)=min(nl,subset0(2));
subset(3)=max( 1,subset0(3));if subset0(3)<1;ofsm=-subset0(3)+1;else;ofsm=0;end
subset(4)=min(ml,subset0(4));
nll=subset(2)-subset(1)+1;mll=subset(4)-subset(3)+1;
%
if subdeg(2)-subdeg(1)>16;inc=floor((subdeg(2)-subdeg(1))/8);
elseif subdeg(2)-subdeg(1)<3;inc=0.2;
elseif subdeg(2)-subdeg(1)<5;inc=0.5;
else;inc=1;end
lon=[ceil(subdeg(1)/inc-0.0001)*inc:inc:subdeg(2)];
lat=[floor(subdeg(3)/inc+0.0001)*inc:-inc:subdeg(4)];
lo1=round((lon-subdeg(1))/dd+1);
la1=round((subdeg(3)-lat)/dd+1);
%
f1=fopen([infile],'r');
if para(1:3)=='sst';
 c1=uint8(zeros(ml0,nl0));
 fseek(f1,nl*2,'bof');
 a1=fread(f1,[nl ml],'uint16');fclose(f1);
 b1=a1(subset(1):subset(2),subset(3):subset(4))'*heda(6)+heda(7)-273.15;
 lnd=(a1(subset(1):subset(2),subset(3):subset(4))>65534)';
 clear a1;
 dat_max=35;dat_min=-5;npt=256;
% dat_max=32;dat_min=20;npt=256;  % for 20020707
 c1(ofsm+[1:mll],ofsn+[1:nll])=uint8((npt-1)*(b1-dat_min)./(dat_max-dat_min).*(1-lnd)+lnd*(npt-2));
 clear lnd b1;
 f1=fopen('GOCC.pal','r');
 if f1>-1;
  ba1=fread(f1,[256 3],'uint8');fclose(f1);
  gg=zeros(npt,3);gg=ba1(round([34:0.7:213]),:)/255;
 else;
  gg=jet(npt);
 end
 gg(1,1:3)=[1,1,1];gg(npt,1:3)=[0.5,0.5,0.5];gg(npt-1,1:3)=[0.4,0.4,0.4];
elseif para(1:3)=='chl';npt=256;
 c1=uint8(zeros(ml0,nl0));
 fseek(f1,nl*2,'bof');
 a1=fread(f1,[nl ml],'uint16');fclose(f1);
 lnd=(a1(subset(1):subset(2),subset(3):subset(4))>65534)';
 out=(a1(subset(1):subset(2),subset(3):subset(4))>65533)'.*(1-lnd);
 msk=(a1(subset(1):subset(2),subset(3):subset(4))<1)';
 b1=log10(a1(subset(1):subset(2),subset(3):subset(4))'*heda(6).*(1-msk)+0.000001*msk);
 clear a1 msk;
 dat_max=2.02;dat_min=-2; 
% dat_max=1.2;dat_min=-2;  % for 20020707
 c1(ofsm+[1:mll],ofsn+[1:nll])=uint8((npt-1)*(b1-dat_min)./(dat_max-dat_min).*(1-lnd).*(1-out)+lnd*(npt-2)+out*(npt-1));
 clear lnd out b1;
 a0=fopen('GOCC.pal','r');
 if a0>-1;
  ba1=fread(a0,[256 3],'uint8');fclose(a0);
  gg=ba1/255;
 else
  gg=jet(npt);
 end
 gg(1,1:3)=[1,1,1];gg(npt,1:3)=[0.5,0.5,0.5];gg(npt-1,1:3)=[0.4,0.4,0.4];
elseif para(1:3)=='RcR' | para(1:3)=='nLw';
 c1=uint8(zeros(ml0,nl0,3));
 fseek(f1,nl,'bof');
 dat_max=255;dat_min=0;
 for j=1:3;
  a1=uint8(fread(f1,[nl ml],'uint8'));
  c1(ofsm+[1:mll],ofsn+[1:nll],j)=a1(subset(1):subset(2),subset(3):subset(4))';
 end
 fclose(f1);
 clear a1;
elseif para(1:3)=='HKr';
 c1=uint8(zeros(ml0,nl0,3));
 fseek(f1,nl,'bof');
 for j=1:3;
  a1=uint8(fread(f1,[nl ml],'uint8'));
  c1(ofsm+[1:mll],ofsn+[1:nll],j)=a1(subset(1):subset(2),subset(3):subset(4))';
 end
 fclose(f1);
 clear a1;
elseif para(1:3)=='QKr';
 c1=uint8(zeros(ml0,nl0,3));
 fseek(f1,nl,'bof');
 for j=1:2;
  a1=uint8(fread(f1,[nl ml],'uint8'));
  c1(ofsm+[1:mll],ofsn+[1:nll],j)=a1(subset(1):subset(2),subset(3):subset(4))';
 end
 fclose(f1);
 c1(ofsm+[1:mll],ofsn+[1:nll],3)=a1(subset(1):subset(2),subset(3):subset(4))';
 clear a1;
end
%
%---------------------- Display image ---------------------------
clf
if ml0>nl0;
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)-6,scps(2)+scps(4)-7.875,6,7.875])
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 para(1:3)=='chl'|para(1:3)=='sst';
 if mll>nll;nsz=min(0.80,0.80*nl0/ml0*1.3125);msz=nsz/1.3125*ml0/nl0;
 else;nsz=min(0.80,0.80*nl0/ml0/1.3125);msz=nsz*1.3125*ml0/nl0;end
else
 if mll>nll;nsz=min(0.87,0.87*nl0/ml0*1.3125);msz=nsz/1.3125*ml0/nl0;
 else;nsz=min(0.87,0.87*nl0/ml0/1.3125);msz=nsz*1.3125*ml0/nl0;end
end
axes('position',[0.08 0.08 nsz msz])
%[c8,map]=rgb2ind(c1,256,'nodither');imshow(c8,map);  %% for RGB on 8-bit display %%
image(c1);
axis on;grid off;axis('ij');axis('equal')
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(['\it{ }',label,', ',para(k1:k2)],'fontsize',15)
%
%----------------------------- Scale bar ------------------------------------
if para(1:3)=='chl'|para(1:3)=='sst';
 axes('position',[0.10+nsz 0.08 0.05 0.50])
 ddd=(dat_max-dat_min)/100;ba=[dat_min:ddd:dat_max];
 bb=uint8(([ba;ba;ba;ba;ba;ba;ba;ba;ba;ba]'-dat_min)/(dat_max-dat_min)*255);
 image(bb);axis on;axis('xy')
 set(gca,'xtick',[ ]);
 if para(1:3)=='chl';
  rbar=[0.01,0.03,0.1,0.3,1,3,10,30];rbar0=((log10(rbar)-dat_min)/ddd)+1;
  xlabel('[mg m^{-3}]','fontsize',10);
 elseif para(1:3)=='sst';
  rbar=[dat_min:(dat_max-dat_min)/5:dat_max];rbar0=[1:20:101];
  xlabel('[^{o}C]','fontsize',10);
 else;
  rbar=[dat_min:(dat_max-dat_min)/5:dat_max];rbar0=[1:20:101];
  xlabel('[mW cm^{-2} str^{-1} um^{-1}]','fontsize',10);
 end
 set(gca,'ytick',rbar0,'yticklabel',rbar,'fontsize',10,'yaxislocation','right');
 grid on
 colormap(gg)
end
%------------------------------------------------------------------------------
%-------------------- PS file ---------------------------------
if f_print(1:2)=='on';
 print('-dpsc',['mod_',para(1:3),'.ps']);
end
%--------------------------------------------------------------
%
else
 disp(['No Input File = ',infile]);
end
%
