llc/matlab/llc_pcol.m

%function h=llc_pcol(t,fldname,iter)
%function h=llc_pcol(t,fldname,iter,'sphere')
%function h=llc_pcol(t,fldname,iter,'m_map arguments')
% with default iter = 1
%function h=llc_pcol(t,fldname)
%function h=llc_pcol(t,fldname,[],'m_map arguments')
% 
% plots 2D scalar fields v on the MITgcm llc (lat-lon-cap) grid with pcolor.
% The structure of tiles is assumed to be created by rdnctiles, e.g.,
% t = rdnctiles({'grid.*','state.*'},{'XG','YG','Eta'},[],'bytile',[]);
% see help rdnctiles for details 
% (addpath ${mitgcm_rootdir}/utils/matlab/gmt) 
% t must contain XG,YG and the field to be plotted as specified in fldname
%
% The optional flag sphere results in a 3D visualization on the sphere
% without any specific projection. Good for debugging.
%
% llc_pcol return a vector h of patch handles.
% 
% If present 'm_map argurments' are passed directly to m_proj to select
% any projection that m_map allows, then m_grid is called without any
% options (requires m_map, obviously). In order to control the input of
% m_grid, use m_ungrid and the m_grid with proper arguments, e.g.
% h=llc_pcol(x,y,v,'stereo','lon',0,'lat',80,'rad',60)
% plots a stereographic map with center at (0E,80N) and a radius of 60
% degrees as in: m_proj('stereo','lon',0,'lat',80,'rad',60); use
% m_ungrid; m_grid('box','off')
% to turn off the box around the plot.
% If you want coast lines, you can add them with m_coast after calling
% llc_pcol.
% Unfortunatly, cylindrical and conic maps are limited to the [-180 180]
% range. 
function h=llc_pcol(varargin)

%$Header: /u/gcmpack/MITgcm_contrib/eh3/llc/matlab/llc_pcol.m,v 1.2 2006/11/10 19:14:14 heimbach Exp $
%$Name:  $

  zlevel = 1;
  ph = [];
  holdstatus=ishold;
  if ~holdstatus; cla; end
  t      =varargin{1};
  if nargin < 2
    error('need at least tile-structure and fldname as input')
  end
  fldname=varargin{2}; 
  if nargin < 3
    iter = 1;
  else
    iter=varargin{3};
    if isempty(iter); iter = 1; end
  end
  % handle map projections
  mapit = 0;
  noproj= 0;
  if nargin > 3
    proj=varargin{4};
    if strcmp('sphere',proj)
      noproj = 1;
    else
      mapit = 1;
    end
  end
  if mapit
    m_proj(varargin{4:end});
  end
  % number of tiles
  ntiles = length(t);
  % loop over tiles
  if noproj
    % no projection at all
    for itile = 1:length(t)
      xc = getfield(t(itile).var,'XG')*pi/180;
      yc = getfield(t(itile).var,'YG')*pi/180;
      fld = getfield(t(itile).var,fldname);
      if ndims(fld) == 3;
        fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
      else
        fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
      end
      [x,y,z] = sph2cart(xc,yc,1);
      ph = [ph;surf(x,y,z,sq(fld))];        
      view(3)
      if (itile == 1); hold on; end
    end
    set(ph,'edgecolor','none')
    shading flat
    axis(gca,'image')
    view(gca,3)
  else
    overlap = 30;
    for itile = 1:length(t)
      fld = getfield(t(itile).var,fldname);
      xc = getfield(t(itile).var,'XG');
      yc = getfield(t(itile).var,'YG');
      if ndims(fld) == 3;
        fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
      else
        fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
      end
      % divide tile into two to avoid 360 degree jump
      ic{1}=find(xc(:)>+overlap);
      ic{2}=find(xc(:)<-overlap);
      if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
      %    if itile==9; keyboard; end
      for k=1:2
        tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
        x = xc; x(ic{k}) = NaN; y = yc; y(ic{k}) = NaN;
        if sum(~isnan(x(:))) > 0
          ph = [ph;pcolor(x',y',sq(tfld)')];
        end
        if (k==1 & itile == 1); hold on; end
      end
      % now do this again with extensions above 180 and below -180 degrees
      xc = getfield(t(itile).var,'XG');
      yc = getfield(t(itile).var,'YG');
      im=find(xc<0); xc(im) = xc(im) + 360;
      ic{1}=find(xc(:)    >+overlap+180);
      ic{2}=find(xc(:)-360<-overlap-180);
      if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
      for k=1:2
        x = xc-(k-1)*360;
        if mapit; x = xc-(k-1)*2*pi; end
        y = yc;
        tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
        x(ic{k}) = NaN;
        y(ic{k}) = NaN;
        if sum(~isnan(x(:))) > 0
          ph = [ph;pcolor(x',y',sq(tfld)')];
        end
%    axis image; shading flat; pause
      end
    end 
    axis(gca,'image')
    set(gca,'box','on','layer','top')
    set(ph,'edgecolor','none')
    if mapit
      m_grid
    else
      axis([-2 2 -1 1]*90)
    end
  end %if
  hold off
  if holdstatus; hold on; end

  if nargout > 0
    h = ph;
  end
  
  return
1	mlosch	1.1	%function h=llc_pcol(t,fldname,iter)
2			%function h=llc_pcol(t,fldname,iter,'sphere')
3			%function h=llc_pcol(t,fldname,iter,'m_map arguments')
4			% with default iter = 1
5			%function h=llc_pcol(t,fldname)
6			%function h=llc_pcol(t,fldname,[],'m_map arguments')
7			%
8			% plots 2D scalar fields v on the MITgcm llc (lat-lon-cap) grid with pcolor.
9			% The structure of tiles is assumed to be created by rdnctiles, e.g.,
10	heimbach	1.2	% t = rdnctiles({'grid.','state.'},{'XG','YG','Eta'},[],'bytile',[]);
11	mlosch	1.1	% see help rdnctiles for details
12			% (addpath ${mitgcm_rootdir}/utils/matlab/gmt)
13			% t must contain XG,YG and the field to be plotted as specified in fldname
14			%
15			% The optional flag sphere results in a 3D visualization on the sphere
16			% without any specific projection. Good for debugging.
17			%
18			% llc_pcol return a vector h of patch handles.
19			%
20			% If present 'm_map argurments' are passed directly to m_proj to select
21			% any projection that m_map allows, then m_grid is called without any
22			% options (requires m_map, obviously). In order to control the input of
23			% m_grid, use m_ungrid and the m_grid with proper arguments, e.g.
24			% h=llc_pcol(x,y,v,'stereo','lon',0,'lat',80,'rad',60)
25			% plots a stereographic map with center at (0E,80N) and a radius of 60
26			% degrees as in: m_proj('stereo','lon',0,'lat',80,'rad',60); use
27			% m_ungrid; m_grid('box','off')
28			% to turn off the box around the plot.
29			% If you want coast lines, you can add them with m_coast after calling
30			% llc_pcol.
31			% Unfortunatly, cylindrical and conic maps are limited to the [-180 180]
32			% range.
33			function h=llc_pcol(varargin)
34
35	mlosch	1.3	%$Header: /u/gcmpack/MITgcm_contrib/eh3/llc/matlab/llc_pcol.m,v 1.2 2006/11/10 19:14:14 heimbach Exp $
36	heimbach	1.2	%$Name: $
37	mlosch	1.1
38			zlevel = 1;
39			ph = [];
40			holdstatus=ishold;
41			if ~holdstatus; cla; end
42			t =varargin{1};
43			if nargin < 2
44			error('need at least tile-structure and fldname as input')
45			end
46			fldname=varargin{2};
47			if nargin < 3
48			iter = 1;
49			else
50			iter=varargin{3};
51			if isempty(iter); iter = 1; end
52			end
53			% handle map projections
54			mapit = 0;
55			noproj= 0;
56			if nargin > 3
57			proj=varargin{4};
58			if strcmp('sphere',proj)
59			noproj = 1;
60			else
61			mapit = 1;
62			end
63			end
64			if mapit
65			m_proj(varargin{4:end});
66			end
67			% number of tiles
68			ntiles = length(t);
69			% loop over tiles
70			if noproj
71			% no projection at all
72			for itile = 1:length(t)
73			xc = getfield(t(itile).var,'XG')*pi/180;
74			yc = getfield(t(itile).var,'YG')*pi/180;
75			fld = getfield(t(itile).var,fldname);
76			if ndims(fld) == 3;
77			fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
78			else
79			fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
80			end
81			[x,y,z] = sph2cart(xc,yc,1);
82			ph = [ph;surf(x,y,z,sq(fld))];
83			view(3)
84			if (itile == 1); hold on; end
85			end
86			set(ph,'edgecolor','none')
87			shading flat
88			axis(gca,'image')
89			view(gca,3)
90			else
91			overlap = 30;
92			for itile = 1:length(t)
93			fld = getfield(t(itile).var,fldname);
94			xc = getfield(t(itile).var,'XG');
95			yc = getfield(t(itile).var,'YG');
96			if ndims(fld) == 3;
97			fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
98			else
99			fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
100			end
101			% divide tile into two to avoid 360 degree jump
102			ic{1}=find(xc(:)>+overlap);
103			ic{2}=find(xc(:)<-overlap);
104			if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
105			% if itile==9; keyboard; end
106			for k=1:2
107			tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
108			x = xc; x(ic{k}) = NaN; y = yc; y(ic{k}) = NaN;
109			if sum(~isnan(x(:))) > 0
110			ph = [ph;pcolor(x',y',sq(tfld)')];
111			end
112			if (k==1 & itile == 1); hold on; end
113			end
114			% now do this again with extensions above 180 and below -180 degrees
115			xc = getfield(t(itile).var,'XG');
116			yc = getfield(t(itile).var,'YG');
117			im=find(xc<0); xc(im) = xc(im) + 360;
118			ic{1}=find(xc(:) >+overlap+180);
119			ic{2}=find(xc(:)-360<-overlap-180);
120			if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
121			for k=1:2
122			x = xc-(k-1)*360;
123			if mapit; x = xc-(k-1)2pi; end
124			y = yc;
125			tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
126			x(ic{k}) = NaN;
127			y(ic{k}) = NaN;
128			if sum(~isnan(x(:))) > 0
129			ph = [ph;pcolor(x',y',sq(tfld)')];
130			end
131			% axis image; shading flat; pause
132			end
133			end
134			axis(gca,'image')
135			set(gca,'box','on','layer','top')
136			set(ph,'edgecolor','none')
137			if mapit
138			m_grid
139			else
140			axis([-2 2 -1 1]*90)
141			end
142			end %if
143			hold off
144			if holdstatus; hold on; end
145
146			if nargout > 0
147			h = ph;
148			end
149
150			return