| 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 |
Wednesday, March 12, 2014
[Matlab] make the t-s diagram with a third field shown in color
Labels:
diagram,
figure,
Matlab,
oceanography,
t-s
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 |
Subscribe to:
Posts (Atom)
