| 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 |
% t = rdnctiles({'grid.*','state.*'},{'XG','YG','Eta'},[],'bytile',[]); |
| 11 |
% 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 |
%$Header: /u/gcmpack/MITgcm_contrib/eh3/llc/matlab/llc_pcol.m,v 1.2 2006/11/10 19:14:14 heimbach Exp $ |
| 36 |
%$Name: $ |
| 37 |
|
| 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)*2*pi; 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph