Find all points on a gridded surface Z that are farther than some distance from an array of points given by (xpt,|ypt|), then replace the distant points in Z with NaN|s or some declared |replacementVal. If no replacementVal is declared, NaN is assumed. Distances are calculated along the x,y plane only.
Z = replacefartherthan(X,Y,Z,xpt,ypt,maxdist) Z = replacefartherthan(...,replacementVal) [Z,ind] replacefartherthan(...)
Z = replacefartherthan(X,Y,Z,xpt,ypt,maxdist) returns gridded Z where all values in Z corresponding to (X,|Y|) grid points farther than maxdist from points given by (xpt,|ypt|) are replaced by NaN. Z must be a matrix corresponding to X and Y. X and Y can be vectors or grids.
Z = replacefartherthan(...,replacementVal) uses replacementVal to replace points in Z farther than maxdist from (xpt,|ypt|). Default replacementVal is NaN.
[Z,ind] replacefartherthan(...) also returns the indices ind of points in Z farther than maxdist from (xpt,|ypt|).
Create a 501x501 surface Z from 28 GPS elevation measurements taken at the points (xpt,|ypt|). After creating a surface from these observations using griddata, you decide that points in Z farther than 150 meters from any observation point are invalid, so you wish to replace them with |NaN|s.
% Your 28 GPS elevation measurements: xpt = 1000*rand(28,1)+1000; ypt = 1000*rand(28,1); zpt = 1200+100*sin(xpt/100)-50*cos(ypt/150); % Create a surface from the elevation measurements: [X,Y,Z] = griddata(xpt,ypt,zpt,1000:2:2000,(1000:-2:0)'); % Plot your original surface and observational data points: figure('position',[100 100 1000 400]) subplot(1,2,1) pcolor(X,Y,Z) colormap(autumn(256)) shading interp cb = colorbar; ylabel(cb,'surface elevation (m)') axis equal hold on plot(xpt,ypt,'b*','markersize',12) title('complete surface')
Now redefine any points in Z that are farther than 100 m from observational data as NaN:
Z2 = replacefartherthan(X,Y,Z,xpt,ypt,150); subplot(1,2,2) pcolor(X,Y,Z2) colormap(autumn(256)) shading interp cb = colorbar; ylabel(cb,'surface elevation (m)') axis equal hold on plot(xpt,ypt,'b*','markersize',12) title('surface w/ points >150 m from obs. removed.')
This function was written by Chad A. Greene of the University of Texas at Austin's Institute for Geophysics on June 2, 2014. This script works quite efficiently on Matlab R2012b thanks to significant input from Matt J., who wrote the heavy-lifting components of this script during some fruitful discussions with José-Luis, Jan Simon, and Sean de Wolski on the Matlab Answers forum here. Thanks for your help, Matt, José, Jan, and Sean.