Wednesday, March 12, 2014

[Matlab] make the t-s diagram with a third field shown in color

function [hc]=theta_sdiag(theta,s,varargin)
% make the t-s plot given temperature and salinity data with ability to show a third field (e.g., density, depth) in color
% need sw_dens.m from the seawater package
% usage:
%         hc=theta_sdiag(theta,s,varargin)
%             varargin:
%               ['color',value] : if value=1, the density will be shown in color;
%                                        otherwise value should be some filed (e.g., depth) with same size as theta/s
%         ['trange',Trange] : the range of temperature shown on the plot (ylim)
%         ['srange',Srange] : the range of salinity shown on the plot (xlim)
%  ['markersize',size-of-marker] : as it tells, the marker size, check scatter for its definition
%     ['caxis',color-range] : the range of the third field shown in color
% e.g.,
%       hc=theta_sdiag(t,s,'color',dep,'caxis',[10 500],'trange',[-2 10],'srange',[31.5 34.8]);

% history:
%
% April 2009: original code by Vihang bhatt
%                    http://www.mathworks.com/matlabcentral/fileexchange/23796-t-s-diagram/content/theta_sdiag.m
% March 2014: enable to show a thrid field in color and more other options by xianmin (xianmin@ualberta.ca)
%                    http://scriptdemo.blogspot.com

if nargin<2
   help theta_sdiag
   return
end

isColor=0;
markerS=16;
if nargin>3
    while(size(varargin,2)>0)
        switch lower(varargin{1})
            case {'iscolor','color'}
                if numel(varargin{2})~=1
                    colorVar=varargin{2};
                    isColor=2;
                else
                    isColor=varargin{2};
                    if isColor==2
                        error('isColor should not equal 2 if no variable specified for the scatter color')
                    end
                end
                varargin(1:2)=[];
            case {'dt','deltat'}
                deltaT=varargin{2};varargin(1:2)=[];
            case {'ds','deltas'}
                deltaS=varargin{2};varargin(1:2)=[];
            case {'tmin','mint'}
                thetamin=varargin{2};varargin(1:2)=[];
            case {'tmax','maxt'}
                thetamax=varargin{2};varargin(1:2)=[];
            case {'trange','tlimit','ylim','tlim'}
                tRange=varargin{2}; thetamax=tRange(2); thetamin=tRange(1); varargin(1:2)=[]; clear tRange
            case {'smin','mins'}
                smin=varargin{2};varargin(1:2)=[];
            case {'smax','maxs'}
                smax=varargin{2};varargin(1:2)=[];
            case {'srange','slimit','xlim','slim'}
                sRange=varargin{2}; smax=sRange(2); smin=sRange(1); varargin(1:2)=[]; clear sRange
            case {'markersize'}
                markerS=varargin{2};varargin(1:2)=[];
            case {'caxis','mycaxis'}
                myCAXIS=varargin{2};varargin(1:2)=[];
            otherwise
        end
    end

end

theta=theta(:);
s=s(:);
if ~exist('smin','var'), smin=min(s)-0.01.*min(s); end
if ~exist('smax','var'), smax=max(s)+0.01.*max(s); end
if ~exist('thetamin','var'), thetamin=min(theta)-0.1*max(theta); end
if ~exist('thetamax','var'), thetamax=max(theta)+0.1*max(theta); end
if ~exist('deltaS','var'), deltaS=0.05; end
if ~exist('deltaT','var'), deltaT=0.25; end

xdim=round((smax-smin)/deltaS+1);
ydim=round((thetamax-thetamin)/deltaT+1);
dens=zeros(ydim,xdim);
thetai=((1:ydim)-1)*deltaT+thetamin;
si=((1:xdim)-1)*deltaS+smin;

for j=1:ydim
    for i=1:xdim
        dens(j,i)=sw_dens(si(i),thetai(j),0);
    end
end

dens=dens-1000;
[c,h]=contour(si,thetai,dens,'k');
clabel(c,h,'LabelSpacing',1000);
xlabel('Salinity','FontWeight','bold','FontSize',12,'fontname','Nimbus Sans L')
ylabel('Theta (^oC)','FontWeight','bold','FontSize',12,'fontname','Nimbus Sans L')
set(gca,'fontname','Nimbus Sans L')

%% plotting scatter plot of theta and s;
hold on;

if isColor==1
   mydens=sw_dens(s,theta,0)-1000;
   hc=scatter(s,theta,markerS,mydens,'o','fill');
   set(hc,'markeredgecolor','none');
elseif isColor==2 % use the depth as color
   if ~exist('myCAXIS','var')
       myCAXIS=[nanmin(colorVar(:)) nanmax(colorVar(:))];
   end
   caxis(myCAXIS);
   hc=scatter(s,theta,markerS,colorVar(:),'o','fill');
   set(hc,'markeredgecolor','none');
else
   hc=scatter(s,theta,'.');
end

if nargout==0
   clear hc
end

Thursday, September 5, 2013

[MATLAB] read the byte-binary ice concentration data from NSIDC

function myice=readgsfcice(fname,fnamemask)
% read the binary ice concentration data from gsfc
% similar to the fortran version here
% data description:
%     http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
% data :
%     ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/
% mask :
%     ftp://sidads.colorado.edu/pub/DATASETS/brightness-temperatures/polar-stereo/tools/masks/
% usage:
%     myice=readgsfcice(ice_filename, mask_filename)
%     myice: ranging 0 to 1, land value is nan if the mask file is specified
%     http://scriptdemo.blogspot.com
 

isMask=1;
if nargin<1
   help readgsfcice
   return
elseif nargin==2
   if ~exist(fnamemask,'file')
      disp(['mask file not found: ',fnamemask])
      isMask=0;
   end
end
if nargin==1
   isMask=0;
end
if ~exist(fname,'file')
   error(['file not found: ',fname])
end

% read the data
fid=fopen(fname,'r');
% read the header
myheader=(char(fread(fid,300,'uchar')))';
nrow=str2num(myheader(13:18));
ncol=str2num(myheader(7:12));
myscale=str2num(myheader(121:126));
%read data block
myice=fread(fid,'uint8');
fclose(fid);

% read the mask file
if isMask==1
   fid=fopen(fnamemask,'r');
   mymask=fread(fid,'uint8');
   fclose(fid);
   myice(mymask==1)=nan;
end
myice=reshape(myice,ncol,nrow)/myscale;

[Fortran] dump the byte-binary ice concentration data from NSIDC

! @ http://scriptdemo.blogspot.com
! tested on  Linux 2.6.32.26-175.fc12.i686
! convert the byte-binary ice-concentration data file from NSIDC to ascii text format
! data description:
!         http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
! data download:
!         ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/

program dumpGSFCIce
integer, parameter                    :: lenHeader=300
character(len=40)                     :: fname, myfmt, txtfname
character(len=lenHeader)          :: header
integer                                    :: numLen, j,fid,npos, nrow, ncol, getnum, myscale,nlen
character(len=1),allocatable,dimension(:)    :: dblock
integer, allocatable,dimension(:)                 :: myice
logical                                                      :: isDebug=.false.

write(*,*) 'input the bin file name to open:'
read(*,'(a)') fname
write(*,*) 'input the txt file name to save:'
read(*,'(a)') txtfname

fid=7
open(fid,file=trim(fname),form='unformatted',access='direct',recl=4*lenHeader)
read(fid,rec=1) header
close(fid)

ncol=getnum(header(7:12))
nrow=getnum(header(13:18))
myscale=getnum(header(121:126))
if (isDebug) then
   write(*,*) '------------------ begin file header info ------------------------'
   write(*,'(a20,a)') ' file name: ', header(127:150)
   write(*,'(a20,a)') ' image title: ', header(151:230)
   write(*,'(a20,a)') ' year: ', header(103:108)
   write(*,'(a20,a)') ' julian day: ', header(109:114)
   write(*,'(a20,a)') ' scale factor: ', header(121:126)
   write(*,'(a20,a)') ' instrument: ', header(55:60)
   write(*,'(a20,a)') ' other info: ', header(231:300)
   write(*,'(a20,i0)')' num of col: ', ncol
   write(*,'(a20,i0)')' num of row: ', nrow
   write(*,'(a20,i0)')' scale factor: ', myscale
   write(*,*) '------------------ end file header info ------------------------'
endif

! declaration
nlen=ncol*nrow
numLen=lenHeader+nlen
allocate(dblock(numLen))
allocate(myice(nlen))

open(fid,file=trim(fname),form='unformatted',access='direct',recl=4*numLen)
read(fid,rec=1) dblock(:)
close(fid)

! convert the 1-byte char to integer
npos=lenHeader
Do j=1,nlen
     npos=npos+1
     myice(j)=ichar(dblock(npos))
Enddo
deallocate(dblock)

! write the file
write(myfmt,'(a,i0,a)') '(',ncol,'i4)'
open(fid,file=trim(txtfname),form='formatted')
write(fid,trim(myfmt)) myice
close(fid)
deallocate(myice)
end program dumpGSFCIce

function getnum(instr)
character(len=6), intent(in) :: instr
character(len=6) :: newstr
integer :: getnum
integer :: i, n
logical :: isnum
n=0
Do i=1,6
     if (isnum(instr(i:i))) then
        n=n+1
        newstr(n:n)=instr(i:i)
     endif
Enddo

if (n==0) then
   getnum=0
else
   read(newstr(1:n),*) getnum
endif
end function getnum

function isnum(myc)
character(len=1), intent(in) :: myc
logical :: isnum
isnum=.false.
if ( (ichar(myc).ge.ichar('0')) .and. (ichar(myc).le.ichar('9'))) isnum=.true.
end function isnum

ShowCalendar