<p><b>mpetersen@lanl.gov</b> 2012-05-11 09:55:34 -0600 (Fri, 11 May 2012)</p><p>Adding new tool: branches/tools/plot_mpas_cross_sections<br>
<br>
Adding a new tool to view cross-sections of mpas output data in<br>
matlab, by finding a path of cells that most closely follows a line in<br>
lat/lon coordinates. This plotting tool has been tested for ocean<br>
output files, but should nearly work for other cores as well. See<br>
README file.<br>
</p><hr noshade><pre><font color="gray">Added: branches/tools/plot_mpas_cross_sections/README
===================================================================
--- branches/tools/plot_mpas_cross_sections/README         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/README        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,33 @@
+plot_mpas_cross_sections README
+
+This is a matlab post-processing tool to view cross sections of
+MPAS-Ocean data. To begin, change the parameters at the top of
+plot_mpas_cross_sections.m. You will need to change the text strings
+wd, sim(i).dir, and sim(i).netcdf_file so that the text string
+
+[wd '/' sim(i).dir '/' sim(i).netcdf_file ] is the file path,
+
+where wd is the working directory and dir is the run directory.
+Details of the section coordinates and variables may be specified in
+plot_mpas_cross_sections.m.
+
+The data files only need to contain a small number of variables.
+You may need to reduce the file size before copying to a local
+machine using:
+
+ncks -v acc_uReconstructMeridional,acc_uReconstructZonal, \
+nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+file_in.nc file_out.nc
+
+After plot_mpas_cross_sections is executed in matlab, plots are placed
+in the docs/*.nc_dir/f subdirectory, and a latex file is created and
+compiled to produce a docs/*.nc_dir/*.nc.pdf file that includes all
+the figures.
+
+This plotting tool has been tested for ocean output files, but should
+nearly work for other cores as well. A few lines will need to be
+changed. We compute depth from the hZLevel variable, and plot cross
+sections with depth increasing downward.
+
+Mark Petersen, MPAS-Ocean Team, LANL, May 2012
Added: branches/tools/plot_mpas_cross_sections/find_cell_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/find_cell_sections.m         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/find_cell_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,208 @@
+function [sectionCellIndex, nCellsInSection, latSection,lonSection, ...
+         depth, latCellDeg,lonCellDeg] = find_cell_sections ...
+ (wd,dir,netcdf_file,sectionText,coord)
+
+% This function reads grid data from an MPAS-Ocean grid or restart
+% netCDF file, and finds a path of cells that connect the endpoints
+% specified in coord. The path is forced to travel through cells
+% that are closest to the line connecting the beginning and end
+% cells.
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% The text string [wd '/' dir '/' netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+% sectionText a cell array with text describing each section
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionCellIndex(maxCells,nSections) cell index of each section
+% nCellsInSection(nSections) number of cells in each section
+% latSection(max(nCellsInSection),nSections) lat coordinates of each section
+% lonSection(max(nCellsInSection),nSections) lon coordinates of each section
+% depth(nVertLevels) depth of center of each layer, for plotting
+% latCellDeg(nCells) lat arrays for all cells
+% lonCellDeg(nCells) lon arrays for all cells
+
+%%%%%%%%%% parameters internal to this function %%%%%%%%%%
+
+% maxCells specifies the maximum number of Cells attempted along
+% the path to the end-cell before stopping with a warning.
+maxCells = 1500;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Read cell and edge data from grid file
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf(['** find_cell_sections, simulation: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+latCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'latCell'));
+lonCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lonCell'));
+cellsOnCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'cellsOnCell'));
+nEdgesOnCell = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'nEdgesOnCell'));
+hZLevel = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'hZLevel'));
+[dimname,nCells]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nCells'));
+[dimname,nVertLevels]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertLevels'));
+netcdf.close(ncid)
+
+nSections = size(coord,1);
+
+% Compute depth of center of each layer, for plotting
+depth(1) = hZLevel(1)/2;
+for i=2:nVertLevels
+ depth(i) = depth(i-1) + 0.5*(hZLevel(i) + hZLevel(i-1));
+end
+
+if (min(lonCell)>-1e-8)
+ for iSection=1:nSections
+ if (coord(iSection,2)<0)
+ coord(iSection,2) = coord(iSection,2) + 360;
+ end
+ if (coord(iSection,4)<0)
+ coord(iSection,4) = coord(iSection,4) + 360;
+ end
+ end
+end
+
+% The native grid variables use:
+% lat varies from -pi/2:pi/2
+% lon varies from -pi:pi
+% convert to degrees for plotting:
+latCellDeg = latCell*180/pi;
+lonCellDeg = lonCell*180/pi;
+
+sectionCellIndex = zeros(maxCells,nSections);
+nCellsInSection = zeros(1,nSections);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Find cells that connect beginning and ending points
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for iSection=1:nSections
+ latCoord = [coord(iSection,1) coord(iSection,3)]/180*pi;
+ lonCoord = [coord(iSection,2) coord(iSection,4)]/180*pi;
+
+ % Find cell closest to start and end coordinates.
+ % The seed cell array simply stores start and end index.
+ minDist = 1e10*ones(1,2);
+ seedCellIndex = zeros(2,nSections);
+ for iCell = 1:nCells
+ for i=1:2
+ dist = sqrt( ...
+         (lonCoord(i) - lonCell(iCell))^2 ...
+        + (latCoord(i) - latCell(iCell))^2);
+ if (dist<minDist(i))
+         minDist(i) = dist;
+         seedCellIndex(i,iSection) = iCell;
+ end
+ end
+ end
+
+ % Rename first cell on route:
+ lonBeginCell = lonCell(seedCellIndex(1,iSection));
+ latBeginCell = latCell(seedCellIndex(1,iSection));
+ beginCellIndex = seedCellIndex(1,iSection);
+
+ % Assign first cell on route to cell index array:
+ sectionCellIndex(1,iSection) = beginCellIndex;
+
+ % Rename last cell on route:
+ lonEndCell = lonCell(seedCellIndex(2,iSection));
+ latEndCell = latCell(seedCellIndex(2,iSection));
+ endCellIndex = seedCellIndex(2,iSection);
+
+ % Assign distance from beginCell to endCell. This is the
+ % distance to 'beat' by candidate neighbor cells.
+ distLastCell = sqrt( ...
+ (lonBeginCell - lonEndCell)^2 ...
+ + (latBeginCell - latEndCell)^2 );
+
+ % loop through cells in search for section
+ for i=1:maxCells
+ additionalCellFound = 0;
+ minDistToLine = 1e10;
+
+ for j = 1:nEdgesOnCell(sectionCellIndex(i,iSection))
+ iCell = cellsOnCell(j,sectionCellIndex(i,iSection));
+
+ if (iCell<nCells+1 && iCell>0)
+         % I am using lat/lon Cartesian distance.
+         % This is distance to the final cell location.
+         dist = sqrt( ...
+         (lonCell(iCell) - lonEndCell)^2 ...
+         + (latCell(iCell) - latEndCell)^2 );
+
+         % check if this cell is closer to the end cell than the
+ % last cell. If so, it is a candidate, and we can continue.
+         if (dist<distLastCell)
+        
+         % distance formula from http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
+         % distToLine = abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) ...
+         % / sqrt((x2-x1)^2 + (y2-y1)^2);
+
+         distToLine = abs(...
+ (lonEndCell -lonBeginCell )*(latBeginCell-latCell(iCell)) ...
+ - (lonBeginCell-lonCell(iCell))*(latEndCell -latBeginCell) ) ...
+         / sqrt((lonEndCell-lonBeginCell)^2 + (latEndCell-latBeginCell)^2);
+
+ % Find the cell that is closest to the line connecting
+ % beginning and ending cells.
+ if (distToLine<minDistToLine)
+ additionalCellFound = 1;
+         distChosenCell = dist;
+         minDistToLine = distToLine;
+         minDistCellIndex = iCell;
+         end
+        
+         end
+ end
+
+ end % nEdgesOnCell
+
+ distLastCell = distChosenCell;
+
+ sectionCellIndex(i+1,iSection) = minDistCellIndex;
+ nCellsInSection(iSection) = i+1;
+
+ if (minDistCellIndex==seedCellIndex(2,iSection))
+ break
+ end
+
+ end % maxCells
+
+ temptext = char(sectionText(iSection));
+ fprintf(['Section %3i, ' temptext ' '],iSection)
+ if (minDistCellIndex==seedCellIndex(2,iSection))
+ fprintf(' complete with %g cells.</font>
<font color="blue">',nCellsInSection(iSection))
+ else
+ fprintf(['WARNING: Did not complete path to end' ...
+         ' cell with %g maxCells. </font>
<font color="blue">'],maxCells)
+ end
+
+end % iSections
+
+maxNumberCells = max(nCellsInSection);
+
+% reduce size of sectionCellIndex array
+sectionCellIndex = sectionCellIndex(1:maxNumberCells,:);
+
+% assign lat/lon on section
+latSection = zeros(maxNumberCells,nSections);
+lonSection = zeros(maxNumberCells,nSections);
+for iSection = 1:nSections
+ for i=1:nCellsInSection(iSection)
+ iCell = sectionCellIndex(i,iSection);
+ latSection(i,iSection) = latCellDeg(iCell);
+ lonSection(i,iSection) = lonCellDeg(iCell);
+ end
+end
+
+fprintf(['</font>
<font color="gray">'])
Added: branches/tools/plot_mpas_cross_sections/load_large_variables.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/load_large_variables.m         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/load_large_variables.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,73 @@
+function [sectionData] = load_large_variables ...
+ (wd,dir,netcdf_file, var_name, var_conv_factor, ...
+ sectionCellIndex, nCellsInSection)
+
+% Load large variables from netcdf file
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% The text string [wd '/' dir '/' netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% var_conv_factor multiply each variable by this unit conversion.
+% sectionCellIndex(maxCells,nSections) cell index of each section
+% nCellsInSection(nSections) number of cells in each section
+
+%%%%%%%%%% output arguments %%%%%%%%%
+% sectionData(nVertLevels,max(nCellsInSection),nSections,nVars)
+% data in each cross-section for each variable
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Load large variables
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+fprintf('</font>
<font color="blue">')
+fprintf(['** load_large_variables simulation: ' dir '</font>
<font color="blue">'])
+
+filename = [wd '/' dir '/' netcdf_file ];
+ncid = netcdf.open(filename,'nc_nowrite');
+
+nAccumulate = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'nAccumulate'));
+
+[dimname,nVertLevels]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nVertLevels'));
+[dimname,nCells]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'nCells'));
+[dimname,nTimeSlices]= netcdf.inqDim(ncid,netcdf.inqDimID(ncid,'Time'));
+nVars = length(var_name);
+nSections = length(nCellsInSection);
+
+maxNumberCells = max(nCellsInSection);
+sectionData = zeros(nVertLevels,maxNumberCells,nSections,nVars);
+
+for iVar=1:nVars
+ temptext = char(var_name(iVar));
+ fprintf(['loading: ' temptext '</font>
<font color="blue">'])
+ if (temptext(1:6)=='ke_acc')
+ % assume zonal and meridional velocity are in index 1 and 2.
+ sectionData(:,:,:,iVar) = sqrt(sectionData(:,:,:,1).^2 + sectionData(:,:,:,2).^2)/2;
+ else
+
+ acc_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,char(var_name(iVar))));
+ mean_var = zeros(nVertLevels, nCells);
+ for i=1:nTimeSlices
+ mean_var = mean_var + nAccumulate(i)*squeeze(acc_var(:,:,i));
+ end
+ mean_var = mean_var/sum(nAccumulate)*var_conv_factor(iVar);
+                        
+ for iSection = 1:nSections
+ for i=1:nCellsInSection(iSection)
+ iCell = sectionCellIndex(i,iSection);
+ for k=1:nVertLevels
+ sectionData(k,i,iSection,iVar) = mean_var(k,iCell);
+ end
+ end
+ end
+ end
+
+end
+netcdf.close(ncid)
+
+fprintf('</font>
<font color="gray">')
+
Added: branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex
===================================================================
--- branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/mpas_sections.head.tex        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,16 @@
+\documentclass[11pt]{article}
+\usepackage{color,amsmath,amsfonts,amssymb,amsthm}
+
+\usepackage[pdftex]{graphicx}
+\DeclareGraphicsExtensions{.pdf, .jpg, .tif, .png}
+
+\textwidth=7.5in
+\textheight=9.8in
+\topmargin=-0.5in
+\headsep=0.0in
+\headheight=0.0in
+\oddsidemargin=-0.25in
+\evensidemargin=-0.25in
+
+\begin{document}
+
Added: branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/plot_mpas_cross_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,205 @@
+% script plot_mpas_cross_sections
+
+% Specify data files, coordinates and text, then call functions
+% to find sections, load data, and plot cross-sections.
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify data files
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% all plots are placed in the f directory. Comment out if not needed.
+unix('mkdir -p f docs');
+
+% The text string [wd '/' sim(i).dir '/' sim(i).netcdf_file ] is the file path,
+% where wd is the working directory and dir is the run directory.
+
+wd = '/var/tmp/mpeterse/runs/todds_runs/P1';
+
+% These files only need to contain a small number of variables.
+% You may need to reduce the file size before copying to a local
+% machine using:
+% ncks -v acc_uReconstructMeridional,acc_uReconstructZonal, \
+% nAccumulate,latVertex,lonVertex,verticesOnEdge,edgesOnVertex,hZLevel,\
+% dvEdge,latCell,lonCell,cellsOnCell,nEdgesOnCell \
+% file_in.nc file_out.nc
+
+sim(1).dir = 'x5.NA.75km_15km';
+sim(1).netcdf_file = ['o.x5.NA.75km_15km.0027-01-01_00.00.00_acc_uReconstruct.nc'];
+
+sim(2).dir = 'x5.NA.75km_15km.leith';
+sim(2).netcdf_file = ['o.x5.NA.75km_15km.leith.0027-01-01_00.00.00_acc_uReconstruct.nc'];
+
+sim(3).dir = 'x1.15km';
+sim(3).netcdf_file = ['o.x1.15km.0018-11-03_00.00.00_acc_uReconstruct.nc'];
+
+sim(4).dir = 'x5.NA.37.5km_7.5km';
+sim(4).netcdf_file = ['o.x5.NA.37.5km_7.5km.0007-01-01_00.00.00_acc_uReconstruct.nc'];
+
+%clear sim
+%sim(1).dir='x1.120km';
+%sim(1).netcdf_file = 'output120km.0001-10-22_12:30:00.nc';
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify section coordinates and text
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% sectionText a cell array with text describing each section
+sectionText = {
+'N Atlantic 26N lat',...
+'N Atlantic 36N lat',...
+'N Atlantic 41N lat',...
+'N Atlantic 46N lat',...
+'N Atlantic 56N lat',...
+'N Atlantic 70W lon',...
+'N Atlantic 65W lon',...
+'N Atlantic 60W lon',...
+'N Atlantic 50W lon',...
+'N Atlantic 40W lon',...
+'N Atlantic 30W lon',...
+'Eq Pacific 140W lon',...
+'Eq Pacific 0 lat ',...
+'Drake Pass 65W lon',...
+'ACC Tasman 147E lon',...
+         };
+
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% Traverse from south to north, and from east to west.
+% Then positive velocities are eastward and northward.
+coord = [...
+ 26 -80 26 -15;... % N Atl Zonal
+ 36 -76 36 -10;... % N Atl Zonal
+ 41 -72 41 -10;... % N Atl Zonal
+ 46 -60 46 -10;... % N Atl Zonal
+ 56 -60 56 -10;... % N Atl Zonal
+ 20 -70 44 -70;... % N Atl Meridional
+ 19 -65 44 -65;... % N Atl Meridional
+ 8.5 -60 46 -60;... % N Atl Meridional
+ 1.8 -50 62 -50;... % N Atl Meridional
+ -3 -40 65 -40;... % N Atl Meridional
+ -5 -30 68.2 -30;... % N Atl Meridional
+ -8 -140 8 -140;... % Eq Pac Meridional
+ 0 140 0 -95;... % Eq Pac Zonal
+ -65 -65 -55 -65;... % Drake
+ -67 147 -43.5 147;... % S Oc Tas
+ ];
+
+% plotDepth(nSections) depth to which to plot each section
+plotDepth = 4000*ones(1,size(coord,1));
+plotDepth(12:13) = 400;
+
+% For plotting, only four plots are allowed per row.
+% Choose sections above for each page.
+% page.name name of this group of sections
+% sectionID section numbers for each row of this page
+page(1).name = 'nAtlantic_Zon_p1';
+page(1).sectionID = [1:4];
+page(2).name = 'nAtlantic_Zon_p2';
+page(2).sectionID = [5];
+
+page(3).name = 'nAtlantic_Mer_p1';
+page(3).sectionID = [6:9];
+page(4).name = 'nAtlantic_Mer_p2';
+page(4).sectionID = [10:11];
+                
+page(5).name = 'eq_Pac_ACC';
+page(5).sectionID = [12:15];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Specify variables to view
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% var_conv_factor multiply each variable by this unit conversion.
+% var_lims(nVars,3) contour line definition: min, max, interval
+
+var_name = {...
+'acc_uReconstructZonal',...
+'acc_uReconstructMeridional',...
+'ke_acc_uReconstruct'};
+
+var_conv_factor = [100 100 1]; % convert m/s to cm/s for velocities
+
+var_lims = [-10 10 2.5; -10 10 2.5; 0 20 2.5];
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Begin main code. Normally this does not need to change.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%close all
+
+for iSim = 1:length(sim)
+
+ fprintf(['**** simulation: ' sim(iSim).dir '</font>
<font color="blue">'])
+ unix(['mkdir -p docs/' sim(iSim).netcdf_file '_dir/f']);
+ fid_latex = fopen('temp.tex','w');
+ fprintf(fid_latex,['%% file created by plot_mpas_cross_sections, ' date '</font>
<font color="black"></font>
<font color="blue">']);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Find cells that connect beginning and end points of section
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ [sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection, sim(iSim).latSection,sim(iSim).lonSection, ...
+ depth, sim(iSim).latCellDeg,sim(iSim).lonCellDeg] ...
+ = find_cell_sections(wd,sim(iSim).dir,sim(iSim).netcdf_file, ...
+ sectionText,coord);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Plot cell section locations on world map
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ sub_plot_cell_sections(sim(iSim).dir,coord, ...
+ sim(iSim).latCellDeg,sim(iSim).lonCellDeg, ...
+ sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection,fid_latex)
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Load large variables from netcdf file
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ [sim(iSim).sectionData] = load_large_variables ...
+ (wd,sim(iSim).dir,sim(iSim).netcdf_file, var_name, var_conv_factor, ...
+ sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %
+ % Plot data on cross-sections
+ %
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ for iPage = 1:length(page)
+
+ sub_plot_cross_sections(sim(iSim).dir,sectionText, ...
+ page(iPage).name, page(iPage).sectionID,sim(iSim).sectionData,depth,...
+ sim(iSim).sectionCellIndex, sim(iSim).nCellsInSection, ...
+ sim(iSim).latSection,sim(iSim).lonSection,coord,plotDepth,...
+ var_name,var_lims,fid_latex)
+
+ end % iPage
+ fprintf(fid_latex,['</font>
<font color="black">\\end{document}</font>
<font color="black"></font>
<font color="gray">']);
+ fclose(fid_latex);
+ dir = ['docs/' sim(iSim).netcdf_file '_dir/'];
+ filename = [sim(iSim).netcdf_file '.tex'];
+ unix(['cat mpas_sections.head.tex temp.tex > ' dir filename ]);
+ unix(['mv f/*jpg ' dir 'f']);
+ cd(dir);
+ unix(['latex ' filename]);
+ cd('../..');
+end % iSim
+
+
Added: branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/sub_plot_cell_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,110 @@
+function sub_plot_cell_sections(dir,coord, ...
+ latCellDeg,lonCellDeg, ...
+ sectionCellIndex, nCellsInSection,fid_latex)
+
+% Plot cell section locations on world map
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir text string, name of simulation
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% latCellDeg(nCells) lat arrays for all cells
+% lonCellDeg(nCells) lon arrays for all cells
+% sectionCellIndex(maxCells,nSections) cell index of each section
+% nCellsInSection(nSections) number of cells in each section
+% fid_latex file ID of latex file
+
+fprintf(['** sub_plot_cell_sections, on figure 1.</font>
<font color="blue">'])
+
+nSections = size(coord,1);
+
+figure(1); clf
+
+if (min(lonCellDeg)>-1e-8)
+ minLon = 0.0;
+ latTrans = 360;
+else
+ minLon = -180.0;
+ latTrans = 0.0;
+end
+
+ % plot topo data of the earth. This is just low-rez one deg
+ % data for visual reference.
+ load('topo.mat','topo','topomap1');
+ if minLon==-180
+ topoNew(:,1:180) = topo(:,181:360);
+ topoNew(:,181:360) = topo(:,1:180);
+ image([-180 180],[-90 90],topoNew,'CDataMapping', 'scaled');
+ else
+ image([0 360],[-90 90],topo,'CDataMapping', 'scaled');
+ end
+
+ colormap(topomap1);
+ set(gca,'YDir','normal')
+
+ hold on
+
+ % world
+ axis tight
+ set(gca,'XTick',30*[-10:12])
+ set(gca,'YTick',15*[-20:20])
+
+ % half world
+ axis([-240+latTrans 0+latTrans -80 70])
+ set(gca,'XTick',20*[-10:20])
+ set(gca,'YTick',10*[-20:20])
+
+ % N Atlantic
+% axis([-90+latTrans -5+latTrans -5 70])
+% set(gca,'XTick',[-100:5:360])
+% set(gca,'YTick',[-90:5:90])
+
+ % Drake passage
+% axis([-90+latTrans,-50+latTrans,-75,-50])
+% set(gca,'XTick',[-100:2:360])
+ % set(gca,'YTick',[-200:2:200])
+
+ % Pacific
+% axis([130 260 -10 10])
+% set(gca,'XTick',[0:1:300])
+% set(gca,'YTick',[-20:.1:20])
+
+
+ % plot cells. This is just done for debugging.
+ %plot(lonCellDeg,latCellDeg,'.y')
+
+ grid on
+
+ for iSection=1:nSections
+ latCoordDeg = [coord(iSection,1) coord(iSection,3)];
+ lonCoordDeg = [coord(iSection,2) coord(iSection,4)];
+
+ h=plot([mod(lonCoordDeg,360)],[latCoordDeg],'*-');
+ set(h,'Color','y','LineWidth',1)
+
+ for i=1:nCellsInSection(iSection)-1
+        h = line([lonCellDeg(sectionCellIndex(i,iSection)) lonCellDeg(sectionCellIndex(i+1,iSection))],...
+                 [latCellDeg(sectionCellIndex(i,iSection)) latCellDeg(sectionCellIndex(i+1,iSection))]);
+        set(h,'Color',[1 0 0],'LineWidth',2)
+ end
+ end
+
+ ylabel('latitude')
+ xlabel('longitude')
+ title(['Domain: ' regexprep(dir,'_','\\_') ' Cells in cross sections. '])
+
+ set(gcf,'PaperPositionMode','auto','color',[.8 1 .8], ...
+         'PaperPosition',[0.25 0.25 8 8])
+
+ subplot('position',[0 .95 1 .05]); axis off
+ text(.005,.7,[ date ]);
+
+ filename=['f/' regexprep(dir,'\.','_') '_cell_map' ];
+ print('-djpeg',[filename '.jpg'])
+
+ % put printing text in a latex file
+ fprintf(fid_latex,...
+ ['\\begin{figure}[btp] \\center </font>
<font color="blue"> \\includegraphics[width=7.5in]{'...
+ filename '.jpg} </font>
<font color="black">\\end{figure} </font>
<font color="gray">']);
Added: branches/tools/plot_mpas_cross_sections/sub_plot_cross_sections.m
===================================================================
--- branches/tools/plot_mpas_cross_sections/sub_plot_cross_sections.m         (rev 0)
+++ branches/tools/plot_mpas_cross_sections/sub_plot_cross_sections.m        2012-05-11 15:55:34 UTC (rev 1898)
@@ -0,0 +1,109 @@
+function sub_plot_cross_sections(dir,sectionText,pageName,sectionID,sectionData,depth,...
+ sectionCellIndex, nCellsInSection, latSection,lonSection, coord, plotDepth,...
+ var_name,var_lims,fid_latex)
+
+% Plot data on cross-sections
+
+% Mark Petersen, MPAS-Ocean Team, LANL, May 2012
+
+%%%%%%%%%% input arguments %%%%%%%%%
+% dir text string, name of simulation
+% sectionText a cell array with text describing each section
+% sectionID section numbers for each row of this page
+% pageName name of this group of sections
+% sectionData(nVertLevels,max(nCellsInSection),nSections,nVars)
+% data in each cross-section for each variable
+% depth(nVertLevels) depth of center of each layer, for plotting
+% sectionCellIndex(maxCells,nSections) cell index of each section
+% nCellsInSection(nSections) number of cells in each section
+% latSection(max(nCellsInSection),nSections) lat coordinates of each section
+% lonSection(max(nCellsInSection),nSections) lon coordinates of each section
+% coord(nSections,4) endpoints of sections, with one section per row as
+% [startlat startlon endlat endlon]
+% plotDepth(nSections) depth to which to plot each section
+% var_lims(nVars,3) contour line definition: min, max, interval
+% var_name(nVars) a cell array with text for each variable to
+% load or compute.
+% fid_latex file ID of latex file
+
+fprintf(['** sub_plot_cross_sections simulation: ' dir ...
+         ' plotting page: ' pageName '</font>
<font color="blue">'])
+
+px = [.28 .78];
+py=linspace(.84,.13,4); % Midpoint position of plots
+pw = [.4]; ph=[.18]; % width and height of plots
+nSections = length(nCellsInSection);
+nVars = length(var_name);
+
+for iVar=1:nVars
+ figure(iVar+1); clf
+ set(gcf,'Position',[100+(iVar*100) 600-iVar*100 715 975])
+ temptext2=char(var_name(iVar));
+
+ for iRow = 1:length(sectionID)
+ iSection = sectionID(iRow);
+ if coord(iSection,1)==coord(iSection,3) % meridional section
+        xtext = 'longitude';
+        xaxis = lonSection(1:nCellsInSection(iSection),iSection);
+ else % zonal section
+        xtext = 'latitude';
+        xaxis = latSection(1:nCellsInSection(iSection),iSection);
+ end
+
+ % left column
+ ha=subplot('position',[px(1)-pw/2 py(iRow)-ph/2 pw ph]);
+ temptext = char(sectionText(iSection));
+ if temptext2(1:6)=='ke_acc'
+ h=surf(xaxis, depth,log10(sectionData(:,1:nCellsInSection(iSection),iSection,iVar)));
+ set(gca,'CLim',[-1 1.2])
+ else
+ h=surf(xaxis, depth,sectionData(:,1:nCellsInSection(iSection),iSection,iVar));
+ end
+
+ set(h,'EdgeColor','none')
+ view(0,-90)
+ title([temptext ', cm/s'])
+ ylabel('depth, m')
+ xlabel(xtext)
+ axis tight
+ set(gca,'YLim',[0 plotDepth(iSection)])
+ h=colorbar ;
+ if temptext2(1:6)=='ke_acc'
+        set(h,'YTick',[-1:1:1.2],'YTickLabel',[0.1 1 10])
+ end
+
+ % right column
+        
+ ha=subplot('position',[px(2)-pw/2 py(iRow)-ph/2 pw ph]);
+ temptext = char(sectionText(iSection));
+ hold on
+ contour(xaxis, depth,sectionData(:,1:nCellsInSection(iSection),iSection,iVar), ...
+         [-30:var_lims(iVar,3):30]);
+ set(gca,'CLim',var_lims(iVar,1:2))
+ set(gca,'YDir','reverse')
+ title([temptext ', ' num2str(var_lims(iVar,3)) 'cm/s int'])
+ ylabel('depth, m')
+ xlabel(xtext)
+ axis tight
+ set(gca,'YLim',[0 plotDepth(iSection)])
+ grid on
+ h=colorbar;
+
+ end
+
+ set(gcf,'PaperPositionMode','auto','color',[.8 1 .8], ...
+ 'PaperPosition',[0.25 0.25 8 10])
+ subplot('position',[0 .95 1 .05]); axis off
+ title_txt = [regexprep(char(var_name(iVar)),'_','\\_') ', ' regexprep(dir,'_','\\_')];
+ h=text(.55,.4,title_txt);
+ set(h,'HorizontalAlignment','center','FontWeight','bold','FontSize',14)
+ text(.005,.7,[ date ]);
+
+ filename=['f/' regexprep(dir,'\.','_') '_' pageName '_var' num2str(iVar)];
+ print('-djpeg',[filename '.jpg'])
+ fprintf(fid_latex,['\\begin{figure}[btp] \\center </font>
<font color="blue"> \\includegraphics[width=7.5in]{'...
+         filename '.jpg} </font>
<font color="black">\\end{figure} </font>
<font color="blue">']);
+ % print('-depsc2',[filename '.eps'])
+
+end
+
</font>
</pre>