Golden Software recently hosted a training class on gridding and interpolating data in Surfer. Before the class was held, a user asked if we’d specifically cover the best options for gridding airborne geophysical data. At the time, it was not in the schedule, but as I looked at the data I thought this type of data could be very common and would make a great example. In this type of data, the data is taken in lines, where the data points along the lines are much closer together than the spacing between the lines. Users generally want to interpolate the data to create a smooth color-filled image map while maintaining the data at sufficient resolution to show important anomalies.
A popular source of airborne geophysical survey data is the USGS aeromagnetic surveys, available at http://mrdata.usgs.gov/magnetic/surveys.php. Some of the data sets have over one million data points, which is quite cumbersome to work with unless you have a lot of memory available on your computer. One data set that is more manageable to work with, and which we’ll use for this example, is for Cuprite, Nevada. This is a tough data set to interpolate, because of the high density of data in one direction (along the flight lines) versus the other direction (between the flight lines).
To find the best gridding method for this (or any) data set, these are the steps I perform when evaluating gridding options:
In detail, the steps I followed with this data set are outlined below.
Step 1: Create a Classed Post Map
The first step I always do is to create a classed post map of the data. I find that it is very useful to get a feel for the spatial distribution of the data, so I know what to expect when gridding the data. For example, areas with fewer data points may generate some odd-looking contours or contours that are significantly different depending upon gridding method. Or if the data points are regularly spaced, I know that Nearest Neighbor might be the best gridding method. In this case, I also want to find out how close the data points are along the lines and between the lines.
Step 2: Run the Script
After getting a feel for the data, I run the script GridData_Comparison.bas with the data file using columns 3, 4, 11 (C, D, K), and compare the results to see the best potential gridding methods. The script grids the data with the eight most popular gridding methods and creates a map from the result. It uses the default options, so it’s not perfect, but gives us a good idea for what the data looks like gridded with those methods.
Right away, I can eliminate Modified Shepard’s Method and Nearest Neighbor because I do not think the interpolation matches what I want the map to look like. I delete those maps. I also do not like how jagged the map is using Radial Basis Function. I delete that map also.
That still leaves five potential gridding method options. To try and eliminate others, I next check how well the contours match the original data.
Step 3: Grid the Data Manually with Options Specific to Data
With the 4 methods left (Inverse Distance, Kriging, Natural Neighbor, and Minimum Curvature), I wanted to grid them with the actual grid spacing of the data to help pick up any anomalies that the coarser default grid spacing missed. I do this manually.
Step 4: Calculate Residuals
At this point, all three remaining maps look pretty similar to me visually, and they all look “nice”. The next step is to see if I can tell a difference between the maps (and hence grids) mathematically, using the residuals. This is important because although I want the map to look nice, I also want it to be the most accurate relative to the data. For this, I will calculate the sum of the residuals for those three grids.
Based on the residuals, the Kriging grid best represents the data with a higher density of grid nodes. So my conclusion is that the best gridding parameters to use for this data set is the Kriging gridding method with a grid node spacing of 0.0001 in the X direction and 0.001 in the Y direction.
Although much time was spent to conclude the best gridding parameters, most of the time was spent eliminating the other methods or parameters. I can now be confident that I’m using an accurate set of gridding parameters to best display the data.
The map above is created with proportional XY scaling (and the maps created by the script are scaled to a particular size, unrelated to their map coordinates). As many of you know, 1° of latitude does not equal 1° of longitude, except at the equator. Therefore, the map above is a little compressed in the Y direction. To scale lat/long maps correctly, multiply the X scaling by the cosine of the latitude, and then multiply the Y scaling by that value. Please see our KB article: How can I scale my Lat/Long maps correctly?
Another option, suggested by a user, is to convert the data from lat/long to UTM (or another coordinate system with linear units) prior to gridding. Having linear, or proportional, XY units could have a slight effect on the search around the grid nodes, or the weight of the data points in the grid node calculation. For instructions on converting the coordinate system of raw data in Surfer, please see our KB article: How can I convert the coordinate system of raw data, such as from UTM to Lat/Long?
For more information about gridding data, please see these additional resources:
Thanks for sharing this example. I had one comment that users may want to consider or keep in mind. Generally I had the understanding that it is best practice to convert your data from a geographic projection (e.g. coordinates in lat/lon in degrees) into some projected coordinate system selected to preserves distance in all directions locally (e.g. a UTM grid) BEFORE you do your gridding and interpolation. The reason for this is that the latitude/longitude grid system is not a uniform grid in terms of distance and the distance represented by a degree of latitude and a degree of longitude are not necessarily the same and are dependent on latitude as the meridians converge at the poles. In your example it appears that you have left the data in lat/lon and the interpolation methods are treating them as if they were in consistent units of length, rather than in degrees. For your example data set from Cuprite, with a latitude of ~37.54 deg, one degree of latitude would be equal to approximately 111 kilometers whereas one degree of longitude would be equal to 88 kilometers. So users should be cognizant of this when applying these techniques to data sets.
You are correct that one degree of latitude is not the same as one degree of longitude.
As you say, you could do the conversion to UTM (or any other coordinate system with linear units) before you grid the data. That may affect the search or the weight of the data points around a node. I'm not sure it would have a significant effect though, but it could. That is a good idea.
Alternatively, for visual appearance you could grid the data in lat/long and rescale the maps. By default, Surfer scales maps proportionally so 1x = 1y. This is how I left the final map in that article (the script scales the maps by size, so that is even more different). However, you can change the scaling of the map easily in Surfer to be more appropriate to lat/long data by multiplying the longitude by the cosine of the latitude (see: http://www.goldensoftware.com/knowledge-base/surfer/1067-how-can-i-scale-my-lat-long-maps-correctly).
I have updated the article with this information! Thank you!
Great Post! Very helpful as I am gridding with similarly acquired ground data. It seems like this whole workflow could be automated into one script for ease of use.