Occam2DMT_Matlab_occam_matlab图像处理_

%-------------------------------------------------------------------------- % plotOccam2DMT.m % % About: A Program to plot up OCCAM2DMT inversion models and data. % % Usage: Call this program from the matlab prompt and GUI window pops % up. Their are model plot options in the figure menu, as well % as menu items for plotting conductivity profiles and MT % responses and data fits. % % Or call as: % % plotOccam2DMT('iterationfilename') % % where 'iterationfilename' is the name of the iteration to % plot. % % You can also use: % % plotOccam2DMT('lastiter') % % and the last (most recent) iteration will be plotted. % % OPTIONAL PARAMS: % 'plotcenter', [xmin xmax ymin ymax] plots only the range given % 'caxis', [cmin cmax] color axis values % 'plotlog', 'off' uses a linear color scale % 'x0', n offsets the x axis by n (uses: x-n % 'sitectr', 'on' center the plot on the sites % 'axes', h axes to plot in (subplotting) - lots of stuff not done % in this case. % % % Written by: Kerry Key % Scripps Institution of Oceanography % kkey@ucsd.edu % % Version: 1.0 August 2001 % 1.1 July 2003 - Made a function with callable iteration file % 1.2 2006 DGM - various cleanups to keep from crashing. % 1.3 April, 2007. Added 'lastiter' option for file name % input % 1.4, October, 2007. All functions are in this single file. % %-------------------------------------------------------------------------- function plotOccam2DMT(varargin) %-------------------------------------------------------------------------- % Plotting OPTIONS: %-------------------------------------------------------------------------- plt_log = 'on'; % Flag for linear or log resistivity scale: plt_cntr = 'off'; % Flag for plotting just the center: plt_caxis = []; % color axis to apply plt_sitescenter = 'off'; % flag for plotting model as wide as sites plot_prej = 'off'; % do not prejudice model plot_zero = 0; % zero point on horizontal axis % kludge for re-zeroing UTM models exceptioncolor = 'w-'; % color of line to plot for model penalty exceptions x0 = 0; plt_axes = []; hfigure = []; sIterFile = []; if nargin==1 sIterFile = varargin{1}; elseif nargin>1 sIterFile = varargin{1}; for i=2:2:nargin str = varargin{i}; switch( lower(str) ) case {'plotcenter'} plt_cntr = 'on'; plot_center = varargin{i+1}; case {'plotlog'} plt_log = varargin{i+1}; case {'x0'} x0 = varargin{i+1}; case {'caxis'} plt_caxis = varargin{i+1}; case {'sitectr'} plt_sitescenter = varargin{i+1}; case {'axes','axis'} plt_axes = varargin{i+1}; case {'useprejudicefile'} plot_prej = 'on'; case {'usefigure'}; hfigure = varargin{i+1}; end end end %-------------------------------------------------------------------------- % GET ITERATION FILENAME AND PATH: %-------------------------------------------------------------------------- if isempty(sIterFile) || ~exist( sIterFile, 'file' ) % If input file string is 'lastiter', then find last iteration in current % directory and use that. if strcmpi(sIterFile,'lastiter') files = dir('*.iter'); % get all iteration files if isempty(files) files = dir('ITER*'); end [temp isort] = sort({files.date}); isort = fliplr(isort); % descending order sIterFile = files(isort(1)).name; [sIterPath,sF,sX] = fileparts(sIterFile); sIterFile = [sF sX]; else % Ask user for file name cFileSpec = {'*.iter','Iteration files (*.iter)';'*.*','All Files'}; [sIterFile,sIterPath] ... = uigetfile( cFileSpec, 'Select OCCAM iteration file to plot:' ); if ~ischar(sIterFile) return end end else [sIterPath,sF,sX] = fileparts(sIterFile); sIterFile = [sF sX]; end %-------------------------------------------------------------------------- % READ ITERATION FILE: %-------------------------------------------------------------------------- [stIter, nParamValues] = readITER( fullfile( sIterPath, sIterFile ) ); %-------------------------------------------------------------------------- % READ MODEL FILE: %-------------------------------------------------------------------------- [fidmod,modfrm,modnam,moddes,mshfile,mshtyp,stcfile,... prjfile,bndofs,nlayer,irz,nrcol,columns,nexep,l1,... c1,l2,c2,pen]= readMODEL( fullfile( sIterPath, stIter.sModelFile ) ); %-------------------------------------------------------------------------- % READ PREJUDICE FILE %-------------------------------------------------------------------------- if strcmp(plot_prej,'on'); if ~strcmpi(prjfile,'none') fid = fopen(prjfile,'r'); sLine = fgets(fid); [sCode, sValue] = strtok( sLine, ':' ); sCode = lower(strtrim(sCode)); sValue(1) = []; % remove the leading token sValue = strtrim( strtok(sValue, '!%') ); if strcmpi(sValue,'OCCAM2MTPREJ_2.0') sLine = fgets(fid); [sCode, sValue] = strtok( sLine, ':' ); sCode = lower(strtrim(sCode)); sValue(1) = []; % remove the leading token sValue = strtrim( strtok(sValue, '!%') ); nprej = sscanf(sValue,'%i'); prej = fscanf(fid,'%i %g %g\n',[3 nprej]); prej = prej'; else plot_prej='off'; end fclose(fid); end end %-------------------------------------------------------------------------- % READ MESH FILE: %-------------------------------------------------------------------------- [fidmsh,mshdesc,idx,nodey,nodez,nres,nfre,nexec,... fixres,freqs,meshy,meshz,nrc,nrl,mshres]... = readMESH( fullfile( sIterPath, mshfile ) ); %-------------------------------------------------------------------------- % READ DATA FILE: %-------------------------------------------------------------------------- [datfrm,datdes,nsites,sites,offsts,nfre,freqs, ndblks,DATA] ... = readDATA( fullfile( sIterPath, stIter.sDataFile ) ); %-------------------------------------------------------------------------- % MAKE FE MESH ELEMENTS: %-------------------------------------------------------------------------- sumy = cumsum([0; meshy]); sumz = cumsum([0; meshz]); % ADJUST SUMY WITH BINDING OFFSET % BNDOFS is the right terminating edge of the left-most reg. brick temp = bndofs - sumy(columns(1,1)+1); sumy = sumy+temp; % remove x0 offset for plotting: sumy = sumy-x0 - plot_zero; sumy = sumy/1000; sumz = sumz/1000; % convert to km offsts = (offsts - plot_zero) / 1000; x0 = x0/1000; [YY, ZZ, restri] = makeELEMENTS(sumy,sort(-1*sumz),mshres); ZZ = -1*ZZ; %-------------------------------------------------------------------------- % MAKE FE MESH LINE