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). PyGMT is an option for those using Python.
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.
If making .csv or .txt files: 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) and make the format for numbers not contain commas. 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.
Alternatively, Matlab is now pretty good about importing Excel native files provided they are straightforward (so not several sheets or plots, etc.). The main wrinkle is making sure you don't include the column titles to be imported (they will map to NaNs).
This should be straightforward if your file is clean. Note that some particulars vary between versions of Matlab; this is referring to 2019a on a Mac. After launching Matlab, choose "Import Data..." (an icon in the menu bar at the top of the window). 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 Import window at the upper left. You should see your data file in the Import window. Choose "Column vectors" from the Output Type: pulldown menu. Now each column should have an editable vector name above the data; you can edit those names if necessary. If the data is selcted except for the first line with the headers, you can click on Import Selection or the big green checkmark. You should see the vector names show up in the Workspace area of the main window. You can now close the Import Window if desired. (If you do not have the column 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). (If you neglect to exclude the column names from the import, there will be NaNs at the top of your file that will cause grief).
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. Some care must be taken if the spread of coordinates is quite large or quite small. So work from what you want: how many points on the Easting coordinate? Call this M. Then, how large is the range of Eastings? Call this R. Then your spacing is R/M = S.
rangeX=min(Easting):S:max(Easting);
rangeY=min(Northing):S:max(Northing);This is often not quite what we want--we often want some extra on the edge (and a lot of times want the coordinates at the edge to not be super funky numbers). When our spacing is close to an integer, this may be accomplished by a pair of commands like
rangeX=floor(min(Easting)):0.1:ceil(max(Easting));
rangeY=floor(min(Northing)):0.1:ceil(max(Northing));You should set the 0.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.
Now if your spacing is way different than 1, you need greater care if you want that buffer (if you don't, just go back up to that earlier equation). Say you are going from 40.000121 to 40.000322 with 100 points. So your spacing would be 0.000201/100 = 0.00000201, so we round to the more even 2x106. Floor and ceil return integers, so you'd end up with points running from 40 to 41 at an 0.000002 interval, or about 500,000 points. Probably not ideal. In this case, we want to use floor and ceil on the digit above our interval, so here 0.00001, or 105. So if we multiply by 105 and then use floor and ceil and divide back by 105, we should be OK:
rangeX=floor(min(Easting)*10^5)/10^5:2*10^-6:ceil(max(Easting)*10^5)/10^5;
rangeY=floor(min(Northing)*10^5)/10^5:2*10^-6:ceil(max(Northing)*10^5)/10^5;This will change 40.000121(the minimum) to 4000012.1, the floor is then 4000012 which when divided by 105, gives 40.00012 for the start. Similarly the maximum goes from 40.000322 to 40.00033. Our spacing 0.000002 from 40.000120 to 40.000330 requires 211 points, which is a fair number...
[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). The polyfit command allows for creation of a polynomial fit to data, which can be useful in noisy datasets.
GEOL4714/5714 public home | C. H. Jones | CIRES | Dept. of Geological Sciences | Univ. of Colorado at Boulder
Last modified at Tuesday, October 14, 2025 5:36 PM
0 visits from