One of the all-purpose software packages for numerical analysis of data is Matlab. It is one of the packages we have on the computers in Benson (and it is also on many of the Engineering labs). Sadly, it is awkward in the way that it will contour irregularly spaced data. (Some other packages such as Aabel and Surfer will contour irregular data much more faithfully, but we lack these packages on lab computers. GMT (Generic Mapping Tools) will do this well and is free, but a command-line Unix package that can be challenging to understand and install; iGMT helps but again is not the easiest thing to install. If you plan to do a lot of mapping, learning GMT is probably worthwhile).
An overview of using Matlab for geological applications in general is at ASU.
To make a regular contour plot, you will first need to make a file that will import smoothly into Matlab, then resample the data to a uniform grid, and then make your contour map.
Assuming you are in Excel with columns for latitude (or northing) and longitude (or easting) and the data you wish to contour (e.g., complete Bouguer anomaly, elevation, etc.), you first want to strip out any text columns (e.g., station names). Matlab has the nasty tendency to decide that you cannot import columns as vectors if there is anything odd in text, so it is best to avoid. Ideally rows or columns with blank values should be dropped or something filled in. You need to leave a set of labels at the top of the file (though remove spaces from these also to avoid trouble later on; these will become your variable names). Export the data as either a comma delimited (csv) or tab-delimted text file.
This should be straightforward if your file is clean. Note that some particulars vary between versions of Matlab; this is referring to 2009b on a Mac. After launching Matlab, choose "Import Data..." (typically in the File menu). Matlab will probably correctly figure out the number of text lines at the start and the delimiter (comma or tab) used in the file, but if not, correct that in the first screen. You should see that there is data, textdata, and colheaders in the preview if all is right. In the second screen (after hitting "Next >") you should have the option to "Create vectors from each column using column names."--select this option. Hit "Finish". (If you do not have the create vectors option, there is probably a mismatch somewhere because of some text down in the number, you are missing column headers or there is a blank cell or something like that).
Ideally you would not do this, but while there is a delaunay command in Matlab, it assigns colors to triangles in unsatisfactory fashion. In gridding, you want to have a grid that is uniformly spaced in X and Y. Thus it is preferable that you have imported data with a Northing and Easting of some kind rather than Latitude and Longitude, though you can work with either. For this tutorial, we will assume you have a vector with distance east called Easting, a vector with distance north called Northing and a vector of data to contour called myData. Note that commands below are typed into the Command Window. Bring this to the front if not visible.
The first task to to make the grid to be used. This may be accomplished by a pair of commands like
rangeX=floor(min(Easting)):.1:ceil(max(Easting));
rangeY=floor(min(Northing)):.1:ceil(max(Northing));You should set the .1 value to the spacing you want between grid points. If you are using latitude and longitude, you should use different spacings for X and Y as a degree of latitude is longer than one of longitude, otherwise your grid will be somewhat distorted.
[If you are unfamiliar with Matlab commands, what we are doing is making two vectors, rangeX and rangeY, which will contain the X and Y positions of each line of nodes. In this case we are making rangeX start at the integer value below the lowest value in the Easting vector (this is floor(min(Easting))) and go to the integer above the highest value of Easting (ceil(max(Easting))) with an interval in between of 0.1. The colons mean to make an evenly spaced vector between the two end values with a spacing of the value in the middle. The semicolon at the end means that we don't want to see the result echoed to the screen.]
We now build a pair of big rectangular matrices that have the X and Y values of each grid point (as opposed to just listing them all).
[Xg,Yg]=meshgrid(rangeX,rangeY);
Xg is the matrix of X (easting) values at all points in the grid and Yg is the matrix of Y (northing) values. If you look at these matrices, you will see they are pretty dull, Xg looking like rangeX repeated a bunch of times up and down and Yg looking like rangeY repeated a bunch side-to-side. We now will grid our original data using the command
dataGrid=griddata(Easting,Northing,myData,Xg,Yg);
The new matrix dataGrid has values interpolated to each grid point from the irregular collection of points defined at Easting(i),Northing(i), myData(i) for all the i points in the original vectors. griddata has a number of options for gridding (you can use linear or cubic interpolation, nearest neighbor, or a special Matlab interpolation 'v4'; linear interpolation is the default. You can get cubic interpolation, for instance, by entering dataGrid=griddata(Easting,Northing,myData,Xg,Yg,'cubic');). You now have the three matrices needed to make contour plots (and other fun plots.
To make a contour plot, you could simply type:
[C,h]=contour(rangeX,rangeY,dataGrid);
The new variables on the left are handy for adding labels and such not. You can add labels interactively by issuing the command clabel(C,h, 'manual'). This lets you pick the spots where contours are labeled. If you prefer, you could specify the contour interval directly; for instance, if you wanted contours at -5, -4, -2, 0, 1, and 2 you could instead type:
[C,h]=contour(rangeX,rangeY,dataGrid,[-5 -4 -2 0 1 2]);
A more elegant way to specify evenly spaced intervals is similar to the way we made the rangeX variable; if we wanted contours ever 0.5 units, we could type:
rangeD=floor(min(myData)):.5:ceil(max(myData));
[C,h]=contour(rangeX,rangeY,dataGrid,rangeD);You might desire to see where the original points were that were used in making the grid. This is accomplished by preventing the plot screen from refreshing (hold on) and then issuing a simple x-y plot command:
hold on;
plot(Easting,Northing,'x')The 'x' indicates to use an x in plotting the points. 'o' would make circles, '+' a plus, etc. You can learn more on how to customize the plot using the help system in Matlab.
To make the plot a properly scaled map, we use the daspectcommand. This forces the axes to be of a specified length. If using distance east and north, then we simply use
daspect([1 1 1]);
However, if we prefer to use latitude and longitude, we would use something like
daspect([111.3 85.3 1]);
as this will scale down the x-axis and y-axis to ratios appropriate at 40° N (note that the values are swapped from the lengths in km of a degree; this is just to make it easier to read. You can view that what you want is the number of degrees/km for each axis, which is ([1/85.3 1/111.3 Z]), where Z is immaterial here, and then we just multiply by both 85.3 and 111.3 to get the command above).
You can add labels to the points, though this usually makes for a real mess. If you had a vector named "Station" in this instance, then adding the text command would add labels to all the points:
text(Easting,Northing,Station)
You can add name-value pairs like 'Color', 'FontSize', 'VerticalAlignment' and 'HorizontalAlignment' names to help. But perhaps the most useful way to use the labels when there are lots of points is to zoom in with the zoom tool in Matlab's figure to see what individual points are. Along those lines, you can also label points by numerical values using the num2str command in the third position in the text command.
With these arrays, there are a bunch of other plots that can be made. Most trivially, contourf will make filled contours plots instead of line plots. surfwill make surface plots that can be rotated and played with (you have to use plot3 to then get the points located on the plot). delaunay and trisurf allow creation of a surface connecting the actual points without interpolation, but you cannot easily turn this into contours. contour3 allows contours to be made in 3-D, and these can be combined with a surface (see, for instance, the bottom of this page).
GEOL4714/5714 home | C. H. Jones | CIRES | Dept. of Geological Sciences | Univ. of Colorado at Boulder
Last modified at Tuesday, February 26, 2019 11:16 AM