Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
96 lines (88 sloc) 2.59 KB
function varargout=topofilt(L,wlen,pono,lmcosi)
% [data,data2]=topofilt(L,wlen,pono,lmcosi)
%
% Makes filtered map of the EGM2008 topography expansion
%
% INPUT:
%
% L Lowpass corner frequency [default: 180]
% wlen Filter length [default: 5]
% degres Grid cell size [default: 1/4 degree]
% pono 0 Data will not be plotted
% 1 Data will be plotted [default]
% lmcosi Coefficients to be filtered [loaded if not supplied]
%
% OUTPUT:
%
% data The data on a grid with requested resolution
% data2 The norm of the gradient of the data
%
% SEE ALSO: GRAVIFILT
%
% Last modified by fjsimons-at-alum.mit.edu, 12/06/2009
% Supply defaults
defval('L',180)
defval('wlen',5)
defval('degres',1/4)
defval('pono',1)
% If the variable doesn't exist
if exist('lmcosi','var')~=1
% Read the file
fnpl=fullfile(getenv('IFILES'),...
'EARTHMODELS','EGM2008','EGM2008_topo_720.mat');
% If the file doesn't exist
if exist(fnpl,'file')~=2
% Make the file
lmcosi=fralmanac('EGM2008_Topography','SHM');
% Restrict it to more manageable size
lmcosi=lmcosi(1:addmup(720),:);
save(fnpl,'lmcosi')
else
load(fnpl)
end
end
defval('percs',[0 5 25 50 75 95 100]);
% Perform the filtering
lmcosif=plmfilt(lmcosi,L,wlen);
% Perform the expansion and/or plot
if pono==0
data=plm2xyz(lmcosif,degres);
else
clf
[ah,ha,H]=krijetem(subnum(2,1));
% Plot the topography
axes(ah(1))
data=plotplm(lmcosif,[],[],4,degres);
t(1)=title(sprintf(...
'Earth''s topography and bathymetry filtered to L = %i',L));
cb(1)=colorbar('hor');
vmd=prctile(data(:),percs);
tix=unique([0 vmd([1 find(percs==50) length(vmd)])]);
set(cb(1),'xtick',tix,'xtickl',round(tix),'xlim',vmd([1 length(vmd)]))
% Calculate the norm of the surface gradient, DT (B.167)
lmcosifdelnorm=lmcosif;
[dems,dels]=addmon(lmcosif(end,1));
% Divide by radius of Earth in meters for scale
R=6378136.3;
lmcosifdelnorm(:,3)=lmcosif(:,3).*sqrt(dels.*(dels+1))/R*100;
lmcosifdelnorm(:,4)=lmcosif(:,4).*sqrt(dels.*(dels+1))/R*100;
% Plot norm of surface gradient which is now in percent
axes(ah(2))
data2=plotplm(lmcosifdelnorm,[],[],4,degres);
t(2)=title(sprintf('Signed norm of the surface gradient (%s)','%'));
cb(2)=colorbar('hor');
vmd2=prctile(data2(:),percs);
tix=unique([vmd2([1 find(percs==50) length(vmd2)])]);
set(cb(2),'xtick',tix,'xtickl',round(tix*10)/10,...
'xlim',vmd2([1 length(vmd2)]))
end
% Cosmetics
fig2print(gcf,'tall')
movev(cb,-0.075)
movev(t,0.075)
set(t,'FontS',15)
longticks([ah cb],2)
figdisp([],L)
% Output only if desired
vars={data,data2};
varargout=vars(1:nargout);