function[varargout]=polysmooth(varargin) %POLYSMOOTH Smoothing scattered 2D data with local polynomial fitting. % % POLYSMOOTH generates a map from scattered data in two dimensions--- % either on the plane or on the sphere---using a local least squares fit % to a polynomial. % % A numerically efficient algorithm is used that avoids explicit loops. % Also, the data are pre-sorted so that different mapping parameters can % be tried out at little computational expense. % % A reference paper for this function is currently in preparation. In % the meantime, for more details see the following book: % % Fan and Gijbels, 1996. Local polynomial modelling and its % applications. Chapman and Hall. % __________________________________________________________________ % % Smoothing on the plane % % Let's say we have an array Z of data is at locations X,Y, where X,Y, % and Z are all arrays of the same size. The problem is to obtain a % mapped field ZHAT on some regular grid specified by arrays XO and YO. % % Calling POLYSMOOTH is a two-step process: % % [DS,XS,YS,ZS]=TWODSORT(X,Y,Z,XO,YO,CUTOFF); % ZHAT=POLYSMOOTH(DS,XS,YS,ZS,B,P); % % Firstly, TWODSORT which returns ZS, a 3D array of data values at each % grid point, sorted by increasing distance DS, and the corresponding % positions XS and YS. Here XO and YO are arrays specifying the bin % center locations of the grid. See TWODSORT for more details. % % Here CUTOFF determines the maximum distance included in the sorting % and should be chosen to be greater than B. % % Secondly, POLYSMOOTH fits a Pth order polynomial at each gridpoint % within a neighborhood specified by the "bandwidth" B. % % Data and grid point locations are presumed to have the same units as % the bandwidth B (e.g., kilometers). % % P may be chosen as P=0 (fit to a constant), P=1 (fit to a plane), or % else P=2 (fit to a parabolic surface). % % POLYSMOOTH can also be used for data on the sphere, as described next. % __________________________________________________________________ % % Smoothing on the sphere % % POLYSMOOTH supports a local polynomial fit on the sphere. As before % this is a two-step process: % % [DS,XS,YS,ZS]=SPHERESORT(LAT,LON,Z,LATO,LONO,CUTOFF); % ZHAT=POLYSMOOTH(DS,XS,YS,ZS,B,P,'sphere') % % Firstly one calls SPHERESORT, which is the analogue of TWODSORT for the % sphere. LATO and LONO are arrays specifying the bin center locations % of the grid. See SPHERESORT for more details. % % Secondly, POLYSMOOTH fits a Pth order polynomial at each gridpoint % within a neighborhood specified by the "bandwidth" B. The bandwidth % here should have units of kilometers. % % Note that SPHERESORT and POLYSMOOTH both assume the sphere to be the % radius of the earth, as specified by the function RADEARTH. % __________________________________________________________________ % % Additional options % % A number of different options are possible, as descibed below. These % are specified with trailing string arguments, which can be given in any % order provided they are after the numeric arguments. % __________________________________________________________________ % % Weighting function % % POLYSMOOTH weights the data points in the vicinity of each grid point % according to some decaying function of distance. % % 1) Parabolic weighting % % POLYSMOOTH(...,'epa'), the default behavior, uses the Epanechnikov % kernel --- a fancy name for a parabolic weighting function. This % function vanishes outside a radius of the bandwidth B. % % The Epanechnikov kernel has the numerical advantage that distances from % far away points do not have to be computed, and is therefore preferred. % % 2) Gaussian weighting: % % POLYSMOOTH(...,'gau') uses a Gaussian weighting function instead, with % a standard deviation of B and truncated, for numerical reasons, at a % radius of 3B. % % When using the Gaussian weighting function you'll need to call TWODSORT % or SPHERESORT with a cutoff distance of at least 3B. % __________________________________________________________________ % % Algorithms % % POLYSMOOTH has two algorithm choices for optimal performance, both of % which give identical answers. % % 1) Speed optimization % % The default algorithm is optimized for speed but uses a great deal of % memory. It solves the least squares problem for all grid points % simultaneously (by directly solving matrix inversions---see MATINV.) % % This can be greater than an order of magnitude faster than the obvious % approach of looping over the grid points. % % 2) Memory optimization % % If the data is so large that memory becomes a limiting factor, use % % ZHAT=POLYSMOOTH(DS,XS,YS,ZS,B,P,'memory'); % % which performs an explicit loop. % __________________________________________________________________ % % Weighted data points % % POLYSMOOTH can also incorporate a weighting factor on the data points. % Let W be an array of positive values the same size as the data array Z. % Each data point is then treated as if it were W data points. % % To use weighted data points, call TWODSORT or SPHERESORT with W added: % % [DS,XS,YS,ZS,WS]=TWODSORT(X,Y,Z,W,XO,YO,CUTOFF); % or [DS,XS,YS,ZS,WS]=SPHERESORT(LAT,LON,Z,W,LATO,LONO,CUTOFF); % % followed by ZHAT=POLYSMOOTH(DS,XS,YS,ZS,WS,B,P); % or ZHAT=POLYSMOOTH(DS,XS,YS,ZS,WS,B,P,'sphere'); % % for data distributed on the plane or on the sphere, respectively. % % For very large datasets, this is useful in first condensing the data % into a manageable size. % __________________________________________________________________ % % Additional output arguments % % [ZHAT,WEIGHT]=POLYSMOOTH(...) also returns the total weight used in % the computation at each grid point. WEIGHT is the same size as ZHAT. % % Small weights mean that only few or far away data points contributed % to the calculation; therefore, the values of ZHAT at these points % are less reliable than at points where the total weight is large. % % [ZHAT,WEIGHT,BETA]=POLYSMOOTH(...) additionally returns the matrix % of the estimated first P derivatives. BETA is 3D and is the same % size as ZHAT and WEIGHT along its first two dimensions. % % BETA(:,:,1) is the estimated zeroth derivative, i.e. the surface, so % this is the same as ZHAT. % __________________________________________________________________ % % Variable bandwidth % % ZHAT=POLYSMOOTH(...,N,P,'variable') uses the parameter N instead of a % fixed bandwidth. The bandwidth now varies over the grid such that N % points fall within one bandwidth distance from every grid point. % % [ZHAT,WEIGHT,BETA,B]=POLYSMOOTH(...,N,P,'variable') returns the % bandwidth B at each grid point, which is the same size as ZHAT. % % The variable bandwidth algorithm can give good results when the data % spacing is highly uneven, particularly when used with a higher-order % (P=1 or P=2) fit. % % When using the variable bandwidth method, be aware that the bandwidth % necessary to include N points may turn out to be larger than CUTOFF, % the maximum distance value given to TWODSORT or SPHERESORT. % __________________________________________________________________ % % 'polysmooth --t' runs some tests. % 'polysmooth --f' generates some sample figures. % % Usage: [ds,xs,ys,zs]=twodsort(x,y,z,xo,yo,cutoff); % [zhat,weight,beta]=polysmooth(ds,xs,ys,zs,b,p); % --or-- % [ds,xs,ys,zs]=spheresort(lat,lon,z,w,lato,lono,cutoff); % [zhat,weight,beta]=polysmooth(ds,xs,ys,zs,b,p,'sphere'); % __________________________________________________________________ % This is part of JLAB --- type 'help jlab' for more information % (C) 2008--2009 J.M. Lilly --- type 'help jlab_license' for details if strcmp(varargin{1}, '--t') polysmooth_test,return elseif strcmp(varargin{1}, '--f') polysmooth_figures,return end str='speed'; kernstr='epanechnicov'; geostr='cartesian'; varstr='constant'; for i=1:4 if ischar(varargin{end}) tempstr=varargin{end}; if strcmp(tempstr(1:3),'epa')||strcmp(tempstr(1:3),'gau') kernstr=tempstr; elseif strcmp(tempstr(1:3),'car')||strcmp(tempstr(1:3),'sph') geostr=tempstr; elseif strcmp(tempstr(1:3),'spe')||strcmp(tempstr(1:3),'loo')||strcmp(tempstr(1:3),'mem') str=tempstr; elseif strcmp(tempstr(1:3),'var')||strcmp(tempstr(1:3),'con') varstr=tempstr; end varargin=varargin(1:end-1); end end R=varargin{end}; B=varargin{end-1}; na=length(varargin)-2; if strcmp(str(1:3),'spe')|strcmp(str(1:3),'mem') d=varargin{1}; x=varargin{2}; y=varargin{3}; z=varargin{4}; if na==5 w=varargin{5}; else w=ones(size(z)); end elseif strcmp(str(1:3),'loo') x=varargin{1}; y=varargin{2}; z=varargin{3}; if na==5 w=ones(size(z)); xo=varargin{4}; yo=varargin{5}; elseif na==6 w=varargin{4}; xo=varargin{5}; yo=varargin{6}; end end if strcmp(str(1:3),'loo') [beta,weight,B]=polysmooth_slow(x,y,z,w,xo,yo,B,R,varstr,kernstr,geostr); elseif strcmp(str(1:3),'spe') [beta,weight,B]=polysmooth_fast(d,x,y,z,w,B,R,varstr,kernstr); elseif strcmp(str(1:3),'mem') [beta,weight,B]=polysmooth_memory(d,x,y,z,w,B,R,varstr,kernstr); else error(['Algorithm type ' str ' is not supported.']) end %Adjustment for complex-valued data if ~allall(isreal(z)) beta(isnan(real(beta(:))))=nan+sqrt(-1)*nan; end varargout{1}=beta(:,:,1); varargout{2}=weight; varargout{3}=beta; varargout{4}=B; %for k=1:max(size(beta,3),nargout) % varargout{k}=squeeze(beta(:,:,k)); %end function[B]=polysmooth_bandwidth(varstr,B,M,N,d,w) if strcmp(varstr(1:3),'con') if size(B,2)==1 B=vrep(B,N,2); end if size(B,1)==1 B=vrep(B,M,1); end else %%Earlier version before weighting % if ndims(d)<3 % B=d(B); % else % B=squeeze(d(:,:,B)); % end N=B; %Compute number of data points with potential for weighted points if ndims(d)<3 B=nan; cumw=cumsum(double(w),1); index=find(cumw>=N,1,'first'); if ~isempty(index) B=squeeze(d(index)); end else B=nan*zeros(size(d,1),size(d,2)); cumw=cumsum(double(w),3); for i=1:size(d,1); for j=1:size(d,2) index=find(cumw(i,j,:)>=N,1,'first'); if ~isempty(index) B(i,j)=squeeze(d(i,j,index)); end end end end end function[beta,weight,B]=polysmooth_fast(ds,xs,ys,zs,ws,B,R,varstr,kernstr) tol=1e15; M=size(xs,1); N=size(xs,2); B=polysmooth_bandwidth(varstr,B,M,N,ds,ws); B=vrep(B,size(xs,3),3); P=sum(0:R+1); matcond=zeros(M,N); beta=nan*zeros(size(xs,1),size(xs,2),P); weight=nan*zeros(size(xs,1),size(xs,2)); %Initial search to exclude right away points too far away ds(isnan(zs)|ds>B)=nan; numgood=squeeze(sum(sum(~isnan(ds),2),1)); index=max(find(numgood~=0)); if ~isempty(index) vindex(ds,xs,ys,zs,ws,B,1:index,3); W=polysmooth_kern_dist(ds,B,kernstr).*ws; %ws for weighted data ponts vswap(ds,xs,ys,zs,W,nan,0); weight=sum(W,3); X=polysmooth_xmat(xs,ys,R); XtimesW=X.*vrep(W,P,4); mat=zeros(M,N,P,P); vect=zeros(M,N,P); for i=1:size(X,4) mat(:,:,:,i)=squeeze(sum(XtimesW.*vrep(X(:,:,:,i),P,4),3)); end vect=squeeze(sum(XtimesW.*vrep(zs,P,4),3)); for i=1:N for j=1:M matcond(j,i)=cond(squeeze(mat(j,i,:,:))); if R==1 if matcond(j,i).^2>tol mat(j,i,:,:)=nan; end elseif R==2 if matcond(j,i)>tol mat(j,i,:,:)=nan; end end end end if P==1 vswap(mat,0,nan); beta=vect./mat; else beta=matmult(matinv(mat),vect,3); % sizemat=size(mat); % mat=reshape(mat,size(mat,1).*size(mat,2),size(mat,3),size(mat,4)); % vect=reshape(vect,size(vect,1).*size(vect,2),size(vect,3)); % nonnani=~isnan(sum(sum(mat,2),3)); % beta=nan*zeros(size(mat,1),size(mat,2)); % % vindex(mat,vect,nonnani,1); % size(mat) % beta(nonnani,:)=matmult(matinv(mat),vect,2); % beta=reshape(beta,sizemat(1),sizemat(2),sizemat(3)); % vsize(beta,beta1) % % aresame(beta1,beta,1e-6) % toc end end weight=matcond; B=B(:,:,1); function[beta,weight,B]=polysmooth_memory(d,x,y,z,w,B,R,varstr,kernstr) tol=1e-10; P=sum(0:R+1); M=size(x,1); N=size(x,2); B=polysmooth_bandwidth(varstr,B,M,N,d,w); beta=nan*zeros(M,N,P); weight=nan*zeros(M,N,2); for j=1:M disp(['Polysmooth computing map for row ' int2str(j) '.']) for i=1:N [dist,xd,yd,zd,wd]=vsqueeze(d(j,i,:),x(j,i,:),y(j,i,:),z(j,i,:),w(j,i,:)); W=polysmooth_kern_dist(dist,B(j,i),kernstr); index=find(W~=0&~isnan(zd)); if ~isempty(index) vindex(xd,yd,zd,wd,index,1); weight(j,i)=sum(W(index)); W=diag(W(index).*wd); %Times wd for weighted dat points X=squeeze(polysmooth_xmat(xd,yd,R)); XW=conj(X')*W; mat=XW*X; matcond=rcond(mat); if matcond>tol; beta(j,i,:)=inv(mat)*XW*zd; end end end end function[beta,weight,B]=polysmooth_slow(x,y,z,w,xo,yo,B0,R,varstr,kernstr,geostr) %Explicit loop without pre-sorting, for testing purposes only M=length(yo); N=length(xo); vsize(x,y,z,xo,yo,B0,R); tol=1e-10; P=sum(0:R+1); beta=nan*zeros(M,N,P); weight=nan*zeros(M,N); B=nan*zeros(M,N); for j=1:M %Note that the slow method loops over the other direction for sphere %than the memory method disp(['Polysmooth computing map for row ' int2str(j) '.']) for i=1:N if strcmp(geostr(1:3),'sph') [xd,yd,dist]=latlon2xy(x,y,xo(i),yo(j)); else xd=x-xo(i); yd=y-yo(j); dist=sqrt(xd.^2+yd.^2); end zd=z; wd=w; [sortdist,sorter]=sort(dist); B(j,i)=polysmooth_bandwidth(varstr,B0,1,1,sortdist,wd(sorter)); W=polysmooth_kern_dist(dist,B(j,i),kernstr); %size(W),size(zd) index=find(W~=0&~isnan(zd)); if ~isempty(index) vindex(xd,yd,zd,wd,index,1); X=squeeze(polysmooth_xmat(xd,yd,R)); weight(j,i)=sum(W(index)); W=diag(W(index).*wd); %Times wd for weighted data points XW=conj(X')*W; mat=XW*X; matcond=rcond(mat); if matcond>tol; beta(j,i,:)=inv(mat)*XW*zd; end end end end if strcmp(geostr(1:3),'sph') beta=permute(beta,[2 1 3]); weight=permute(weight,[2 1]); B=permute(B,[2 1]); end function[X]=polysmooth_xmat(x,y,R) P=sum(0:R+1); X=ones(size(x,1),size(x,2),size(x,3),P); if R>0 X(:,:,:,2)=x; X(:,:,:,3)=y; end if R>1 X(:,:,:,4)=frac(1,2)*x.^2; %See writeup X(:,:,:,5)=x.*y; X(:,:,:,6)=frac(1,2)*y.^2; %See writeup end % function[W]=polysmooth_kern_xy(x,y,Binv,kernstr); % % % if size(Binv,4)==1 % xp=Binv.*x; % yp=Binv.*y; % elseif size(Binv,4)==2 % xp=Binv(:,:,:,1).*x; % yp=Binv(:,:,:,2).*y; % %elseif size(Binv,3))==4 % % xp=Binv(1,1).*x+Binv(1,2).*y; % % yp=Binv(2,1).*x+Binv(2,2).*y; % end % % distsquared=abs(xp).^2+abs(yp).^2; % % if strcmp(kernstr(1:3),'epa') % W=frac(2,pi).*(1-distsquared).*(1+sign(1-sqrt(distsquared))); % elseif strcmp(lower(kernstr(1:3)),'gau') % W=frac(2,pi).*exp(-frac(distsquared,2));%Not right, should have frac(1,sigx*sigy) % end function[W]=polysmooth_kern_dist(dist,B,kernstr) % if iscell(B) % N=B{2}; % B=B{1}; % else % N=3; % end N=3; dist=frac(dist,B); if strcmp(kernstr(1:3),'epa') W=frac(2,pi).*(1-dist.^2).*(1+sign(1-dist)); elseif strcmp(lower(kernstr(1:3)),'gau') W=frac(2,pi.*B.^2).*exp(-frac(dist.^2,2)).*(1+sign(1-dist./N)); end W=vswap(W,nan,0); function[]=polysmooth_test polysmooth_test_cartesian; polysmooth_test_sphere; function[]=polysmooth_test_cartesian %Use peaks for testing... random assortment [x,y,z]=peaks; index=randperm(length(z(:))); index=index(1:200); [xdata,ydata,zdata]=vindex(x(:),y(:),z(:),index,1); xo=(-3:.5:3); yo=(-3:.6:3); B=2; [ds,xs,ys,zs]=twodsort(xdata,ydata,zdata,xo,yo,B); tic;z1=polysmooth(xdata,ydata,zdata,xo,yo,B,0,'epan','loop');etime1=toc; tic;z2=polysmooth(ds,xs,ys,zs,B,0,'epan','memory');etime2=toc; tic;z3=polysmooth(ds,xs,ys,zs,B,0,'epan','speed');etime3=toc; disp(['POLYSMOOTH was ' num2str(etime1./etime3) ' times faster than loop, constant fit.']) disp(['POLYSMOOTH was ' num2str(etime2./etime3) ' times faster than memory method, constant fit.']) tol=1e-8; b1=aresame(z1,z3,tol); b2=aresame(z2,z3,tol); reporttest('POLYSMOOTH output size matches input size',aresame(size(z1),[size(xs,1) size(xs,2)])) reporttest('POLYSMOOTH speed and loop methods are identical for constant fit',b1) reporttest('POLYSMOOTH speed and memory methods are identical for constant fit',b2) tic;[z1,w1,beta1]=polysmooth(xdata,ydata,zdata,xo,yo,B,1,'epan','loop');etime1=toc; tic;[z2,w2,beta2]=polysmooth(ds,xs,ys,zs,B,1,'epan','memory');etime2=toc; tic;[z3,w3,beta3]=polysmooth(ds,xs,ys,zs,B,1,'epan','speed');etime3=toc; disp(['POLYSMOOTH was ' num2str(etime1./etime3) ' times faster than loop, linear fit.']) disp(['POLYSMOOTH was ' num2str(etime2./etime3) ' times faster than memory method, linear fit.']) tol=1e-8; b1=aresame(beta1,beta3,tol); b2=aresame(beta2,beta3,tol); reporttest('POLYSMOOTH speed and loop methods are identical for linear fit',b1) reporttest('POLYSMOOTH speed and memory methods are identical for linear fit',b2) function[]=polysmooth_test_sphere %Use peaks for testing... random assortment [x,y,z]=peaks; index=randperm(length(z(:))); index=index(1:200); %Convert to Lat and Longitude lon=x.*60; lat=y.*30; [latdata,londata,zdata]=vindex(lat(:),lon(:),z(:),index,1); wdata=10*abs(rand(size(zdata))); %For heavy data point test lono=(-180:20:180); lato=(-90:20:90); B=2000; [ds,xs,ys,zs,ws]=spheresort(latdata,londata,zdata,wdata,lato,lono,B,'loop'); tic;z1=polysmooth(latdata,londata,zdata,lato,lono,B,0,'epan','sphere','loop');etime1=toc; tic;z2=polysmooth(ds,xs,ys,zs,B,0,'epan','memory');etime2=toc; tic;z3=polysmooth(ds,xs,ys,zs,B,0,'epan','speed');etime3=toc; disp(['POLYSMOOTH was ' num2str(etime1./etime3) ' times faster than loop, constant fit on a sphere.']) disp(['POLYSMOOTH was ' num2str(etime2./etime3) ' times faster than memory method, constant fit on a sphere.']) tol=1e-8; b1=aresame(z1,z3,tol); b2=aresame(z2,z3,tol); reporttest('POLYSMOOTH output size matches input size on a sphere',aresame(size(z1),[size(xs,1) size(xs,2)])) reporttest('POLYSMOOTH speed and loop methods are identical for constant fit on a sphere',b1) reporttest('POLYSMOOTH speed and memory methods are identical for constant fit on a sphere',b2) tic;z1=polysmooth(latdata,londata,zdata,wdata,lato,lono,B,0,'epan','sphere','loop');etime1=toc; tic;z2=polysmooth(ds,xs,ys,zs,ws,B,0,'epan','memory');etime2=toc; tic;z3=polysmooth(ds,xs,ys,zs,ws,B,0,'epan','speed');etime3=toc; tol=1e-8; b1=aresame(z1,z3,tol); b2=aresame(z2,z3,tol); reporttest('POLYSMOOTH speed and loop methods are identical with heavy data points',b1) reporttest('POLYSMOOTH speed and memory methods are identical with heavy data points',b2) B0=10; [ds,xs,ys,zs]=spheresort(latdata,londata,zdata,lato,lono,inf); tic;[z1,weight1,beta1,B1]=polysmooth(latdata,londata,zdata,lato,lono,B0,0,'epan','sphere','variable','loop');etime1=toc; tic;[z2,weight2,beta2,B2]=polysmooth(ds,xs,ys,zs,B0,0,'epan','memory','variable');etime2=toc; tic;[z3,weight3,beta3,B3]=polysmooth(ds,xs,ys,zs,B0,0,'epan','speed','variable');etime3=toc; disp(['POLYSMOOTH was ' num2str(etime1./etime3) ' times faster than loop, constant 10-point fit on a sphere.']) disp(['POLYSMOOTH was ' num2str(etime2./etime3) ' times faster than memory method, constant 10-point fit on a sphere.']) tol=1e-8; b1=aresame(z1,z3,tol); b2=aresame(z2,z3,tol); reporttest('POLYSMOOTH output size matches input size on a sphere',aresame(size(z1),[size(xs,1) size(xs,2)])) reporttest('POLYSMOOTH speed and loop methods are identical for constant 10-point fit on a sphere',b1) reporttest('POLYSMOOTH speed and memory methods are identical for constant 10-point fit on a sphere',b2) function[]=polysmooth_figures disp('Example of polysmooth using Matlab''s ''Peaks'' function comparing constant, linear, and quadratic fits.') %Use peaks for testing... random assortment [x,y,z]=peaks; index=randperm(length(z(:))); index=index(1:600); [xdata,ydata,zdata]=vindex(x(:),y(:),z(:),index,1); xo=(-3:.125:3); yo=(-3:.125:3); B=1; [ds,xs,ys,zs]=twodsort(xdata,ydata,zdata,xo,yo,B); z0=polysmooth(ds,xs,ys,zs,B,0,'epan'); z1=polysmooth(ds,xs,ys,zs,B,1,'epan'); z2=polysmooth(ds,xs,ys,zs,B,2,'epan'); figure subplot(2,2,1),contourf(xo,yo,z,(-10:10)),title('Gaussian Topography') subplot(2,2,2),contourf(xo,yo,z0,(-10:10)),title('Constant Fit') subplot(2,2,3),contourf(xo,yo,z1,(-10:10)),title('Linear Fit') subplot(2,2,4),contourf(xo,yo,z2,(-10:10)),title('Quadratic Fit') for i=1:4 subplot(2,2,i),hold on,caxis([-7 7]),nocontours,plot(xdata,ydata,'k+') plot(B*rot(0:.05:2*pi+.1),'k'),axis square,axis equal,axis tight xtick(-3:1:3),ytick(-3:1:3) colorbar('vertical') end figure subplot(2,2,1),contourf(xo,yo,z,(-10:10)),caxis([-7 7]),title('Gaussian Topography') subplot(2,2,2),contourf(xo,yo,z0-z,(-10:10)/3),caxis([-3 3]),title('Constant Deviation') subplot(2,2,3),contourf(xo,yo,z1-z,(-10:10)/3),caxis([-3 3]),title('Linear Deviation') subplot(2,2,4),contourf(xo,yo,z2-z,(-10:10)/3),caxis([-3 3]),title('Quadratic Deviation') for i=1:4 subplot(2,2,i),hold on,nocontours,plot(xdata,ydata,'k+') plot(B*rot(0:.05:2*pi+.1),'k'),axis square,axis equal,axis tight xtick(-3:1:3),ytick(-3:1:3) colorbar('vertical') end %orient portrait %print -dpng polysmooth_peaks_deviation.png figure, subplot(1,3,1),plot(z,z0,'.'),title('Constant Fit') subplot(1,3,2),plot(z,z1,'.'),title('Linear Fit') subplot(1,3,3),plot(z,z2,'.'),title('Quadratic Fit') for i=1:3 subplot(1,3,i), axis square,axis equal,axis([-7 7 -7 7]) xtick(-6:2:6),ytick(-6:2:6),hlines(0,'k:'),vlines(0,'k:') plot([-1-sqrt(-1);1+sqrt(-1)]*7,'k') xlabel('Original Surface'),ylabel('Estimated Surface') end %orient portrait %print -dpng polysmooth_peaks_scatter.png function[]=polysmooth_test2 %[indexfull]=twodsort(xdata,ydata,xbin,ybin) load labsea_topex use labsea_topex mssh=vmean(ssh,3); %mssh=squeeze(ssh(:,:,137)); std=vstd(ssh,3); lato=(52:1/4:65); lono=(-65:1/4:-40); B=200; [ds,xs,ys,zs,index]=spheresort(lat,lon,mssh,lato,lono,B); [ds,xs,ys,stds,index]=spheresort(lat,lon,std,lato,lono,B); eke=ssh2eke(lat,atd,mssh,10); [ds,xs,ys,ekes,index]=spheresort(lat,lon,eke,lato,lono,B); stds=nan*zeros(size(ds)); stds(~isnan(index))=std(nonnan(index)); [z0,w0,beta0,r0]=polysmooth(ds,xs,ys,stds,24,0,'epan','speed','variable'); stds=nan*zeros(size(ds)); size(index) size(ekes) [z0,w0,beta0,r0]=polysmooth(ds,xs,ys,ekes,24,0,'epan','speed','variable'); [z0,w0,beta0,r0]=polysmooth(ds,xs,ys,zs,24,0,'epan','speed','variable'); eke=ssh2eke(lat,atd,ssh,N); z2=polysmooth(ds,xs,ys,zs,24,2,'epan','speed','variable'); plot(lato,vmean(r0,2)) pcolor(lono,lato,r0),hold on,plot(lon,lat),shading flat z0=polysmooth(ds,xs,ys,zs,100,0,'epan','speed'); %z1=polysmooth(ds,xs,ys,zs,100,1,'epan','speed'); z2=polysmooth(ds,xs,ys,zs,150,2,'epan','speed'); z3=polysmooth(ds,xs,ys,zs,100,0,'epan'); [z2,w2,beta2,r2]=polysmooth(ds,xs,ys,zs,96,2,'epa','speed','variable'); fc=corfreq(56.5)./3600*2*pi; g=9.81; v=frac(100,1000*100)*frac(g,fc).*beta0(:,:,2); u=-frac(100,1000*100)*frac(g,fc).*beta0(:,:,3); figure,contourf(xbin,ybin,squeeze(b1(:,:,1)),(0:.1:2.5)),hold on,plot(xdata,ydata) beta0=polysmooth(xdata,ydata,mssh,xbin,ybin,B,0); beta1=polysmooth(xdata,ydata,mssh,xbin,ybin,B,1); B=diag([150 150]); [beta2,weight]=polysmooth(xdata,ydata,mssh,xbin,ybin,B,2); z0=squeeze(beta0(:,:,1)); z1=squeeze(beta1(:,:,1)); z2=squeeze(beta2(:,:,1)); figure,contourf(xbin,ybin,z0,(0:.1:2.5)),hold on,plot(xdata,ydata) figure,contourf(xbin,ybin,z1,(0:.1:2.5)),hold on,plot(xdata,ydata) figure,contourf(xbin,ybin,z2,(0:.1:2.5)),hold on,plot(xdata,ydata) figure, contourf(xbin,ybin,sqrt(squeeze(beta2(:,:,2)).^2+squeeze(beta2(:,:,3)).^2),(0:.05:1)/50) hold on,plot(xdata,ydata) pcolor(xbin,ybin,z1-z0) caxis([-1 1]/10) hold on,plot(xdata,ydata) %reporttest('POLYSMOOTH',aresame())