Warning, /utils/matlab/cs_grid/latloncap/llc_pcol.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
4bc51130a8 Mart*0001 %function h=llc_pcol(t,fldname,iter)
0002 %function h=llc_pcol(t,fldname,iter,'sphere')
0003 %function h=llc_pcol(t,fldname,iter,'m_map arguments')
0004 % with default iter = 1
0005 %function h=llc_pcol(t,fldname)
0006 %function h=llc_pcol(t,fldname,[],'m_map arguments')
0007 %
0008 % plots 2D scalar fields v on the MITgcm llc (lat-lon-cap) grid with pcolor.
0009 % The structure of tiles is assumed to be created by rdnctiles, e.g.,
0010 % t = rdnctiles({'grid.*','state.*'},{'XG','YG','Eta'},[],'bytile',[]);
0011 % see help rdnctiles for details
0012 % (addpath ${mitgcm_rootdir}/utils/matlab/gmt)
0013 % t must contain XG,YG and the field to be plotted as specified in fldname
0014 %
0015 % The optional flag sphere results in a 3D visualization on the sphere
0016 % without any specific projection. Good for debugging.
0017 %
0018 % llc_pcol return a vector h of patch handles.
0019 %
0020 % If present 'm_map argurments' are passed directly to m_proj to select
0021 % any projection that m_map allows, then m_grid is called without any
0022 % options (requires m_map, obviously). In order to control the input of
0023 % m_grid, use m_ungrid and the m_grid with proper arguments, e.g.
0024 % h=llc_pcol(x,y,v,'stereo','lon',0,'lat',80,'rad',60)
0025 % plots a stereographic map with center at (0E,80N) and a radius of 60
0026 % degrees as in: m_proj('stereo','lon',0,'lat',80,'rad',60); use
0027 % m_ungrid; m_grid('box','off')
0028 % to turn off the box around the plot.
0029 % If you want coast lines, you can add them with m_coast after calling
0030 % llc_pcol.
0031 % Unfortunatly, cylindrical and conic maps are limited to the [-180 180]
0032 % range.
0033 function h=llc_pcol(varargin)
0034
0035 zlevel = 1;
0036 ph = [];
0037 holdstatus=ishold;
0038 if ~holdstatus; cla; end
0039 t =varargin{1};
0040 if nargin < 2
0041 error('need at least tile-structure and fldname as input')
0042 end
0043 fldname=varargin{2};
0044 if nargin < 3
0045 iter = 1;
0046 else
0047 iter=varargin{3};
0048 if isempty(iter); iter = 1; end
0049 end
0050 % handle map projections
0051 mapit = 0;
0052 noproj= 0;
0053 if nargin > 3
0054 proj=varargin{4};
0055 if strcmp('sphere',proj)
0056 noproj = 1;
0057 else
0058 mapit = 1;
0059 end
0060 end
0061 if mapit
0062 m_proj(varargin{4:end});
0063 end
0064 % number of tiles
0065 ntiles = length(t);
0066 % loop over tiles
0067 if noproj
0068 % no projection at all
0069 for itile = 1:length(t)
0070 xc = getfield(t(itile).var,'XG')*pi/180;
0071 yc = getfield(t(itile).var,'YG')*pi/180;
0072 fld = getfield(t(itile).var,fldname);
0073 if ndims(fld) == 3;
0074 fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
0075 else
0076 fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
0077 end
0078 [x,y,z] = sph2cart(xc,yc,1);
0079 ph = [ph;surf(x,y,z,sq(fld))];
0080 view(3)
0081 if (itile == 1); hold on; end
0082 end
0083 set(ph,'edgecolor','none')
0084 shading flat
0085 axis(gca,'image')
0086 view(gca,3)
0087 else
0088 overlap = 30;
0089 for itile = 1:length(t)
0090 fld = getfield(t(itile).var,fldname);
0091 xc = getfield(t(itile).var,'XG');
0092 yc = getfield(t(itile).var,'YG');
0093 if ndims(fld) == 3;
0094 fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,iter);
0095 else
0096 fld = fld(1:size(xc,1)-1,1:size(yc,2)-1,zlevel,iter);
0097 end
0098 % divide tile into two to avoid 360 degree jump
0099 ic{1}=find(xc(:)>+overlap);
0100 ic{2}=find(xc(:)<-overlap);
0101 if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
0102 % if itile==9; keyboard; end
0103 for k=1:2
0104 tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
0105 x = xc; x(ic{k}) = NaN; y = yc; y(ic{k}) = NaN;
0106 if sum(~isnan(x(:))) > 0
0107 ph = [ph;pcolor(x',y',sq(tfld)')];
0108 end
0109 if (k==1 & itile == 1); hold on; end
0110 end
0111 % now do this again with extensions above 180 and below -180 degrees
0112 xc = getfield(t(itile).var,'XG');
0113 yc = getfield(t(itile).var,'YG');
0114 im=find(xc<0); xc(im) = xc(im) + 360;
0115 ic{1}=find(xc(:) >+overlap+180);
0116 ic{2}=find(xc(:)-360<-overlap-180);
0117 if mapit; [xc yc]=m_ll2xy(xc,yc,'clipping','on'); end
0118 for k=1:2
0119 x = xc-(k-1)*360;
0120 if mapit; x = xc-(k-1)*2*pi; end
0121 y = yc;
0122 tfld = fld([1:end 1],[1:end 1]); tfld(ic{k}) = NaN;
0123 x(ic{k}) = NaN;
0124 y(ic{k}) = NaN;
0125 if sum(~isnan(x(:))) > 0
0126 ph = [ph;pcolor(x',y',sq(tfld)')];
0127 end
0128 % axis image; shading flat; pause
0129 end
0130 end
0131 axis(gca,'image')
0132 set(gca,'box','on','layer','top')
0133 set(ph,'edgecolor','none')
0134 if mapit
0135 m_grid
0136 else
0137 axis([-2 2 -1 1]*90)
0138 end
0139 end %if
0140 hold off
0141 if holdstatus; hold on; end
0142
0143 if nargout > 0
0144 h = ph;
0145 end
0146
0147 return
