| 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 |
|
|
% 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: $ |
| 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 |