Gridding and Contouring Airborne Geophysical Survey Data

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:

1. Create a classed post map to visualize the data distribution.
2. Run the sample script GridData_Comparison.bas to compare gridding methods and narrow down the gridding methods to just a few.
3. Manually grid the data with the gridding methods previously selected with more specific gridding options.
4. Calculate the residuals to determine the most accurate method (optional).

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.

1. Click Map | New | Classed Post Map, select NV_1152.dat and click Open.
2. Select the Classed Post layer in the Object Manager, and in the Property Manager, on the General page, change the X,  Y and Z columns to C, D and K, respectively. I noticed the data points were very close together in the X direction, but father apart in the Y direction. 
3. Click the Classes tab in the Property Manager.
4. Click the Edit Classes button.
5. In the Classes for Map dialog:
   1. Set the Binning method to Equal Intervals.
   2. Set the Number of classes to 7.  
   3. Note that interval between classes is about 23, the minimum of the first class is about -288, and the minimum of the last class is about -154.
   4. Click the Symbol button. In the Class Symbol Properties dialog, set the Line color to Black with 10% opacity, and change the symbol itself to a square shape. Click OK.
   5. Click the Size button.  Set the Minimum size and Maximum size to the same small value, 0.005, and click OK.
   6. Click each of the symbols for each class and change the Fill color for that class. I chose navy blue, cyan, green, yellow, red, magenta and pink for my classes.
   7. Click Save and save the class definitions to a CLS file (ClassedPostClasses.cls) for future use.
   8. Click OK.
6. Then, I measure the distance of points. I zoom in to the map so I could see the points individually along the lines and click Map | Measure.
   1. In the Y direction, the flight lines were spaced apart about 0.005° on average.
   2. In the X direction, the points were spaced apart about every 0.0001° on average.
7. Press the ESC key to exit measuring mode.
8. I want to make the symbol sizes a little larger now. Click the Edit Classes button again and click the Size button.
9. Set the Minimum size and Maximum size to the same small value, 0.02, and click OK and OK.

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.

1. Open Scripter by clicking Windows Start | All Programs | Golden Software Surfer 13 | Scripter.
2. Once Scripter is opened, click File | Open, navigate to the Surfer 13 SamplesScripts directory (within the installation folder), select GridData_Comparison.bas and click Open.
3. Click Script | Run.
4. When prompted for the data, select NV_1152.dat and click Open.
5. When asked for the X data column, enter 3 and click OK.
6. When asked for the Y data column, enter 4 and click OK.
7. When asked for the Z data column, enter 11 and click OK. The script opens a new Surfer window and creates a map in it for each of the eight gridding 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.

1. Click back to the plot window containing the classed post map previously created, copy it, click back to the window with the contour maps in it, and paste it five times.
2. Overlay each of the classed post maps with each of the remaining contour maps by dragging a Classed Post layer into each of the maps above the Contour layers in the Object Manager. 
3. Adjust the contours to match the classed post map colors. To do this, select one of the contour layers and click on the Levels tab in the Property Manager. On the Levels page:
   1. Change the Level method to Advanced, and click Edit Levels.
   2. Click the Level button, set the Minimum to -288 and the Interval to 23 (just like the classed post classes) and click OK.
   3. Double click on each of the levels under the Fill column and change the color fill for each level to the same as that of the classed post classes (navy blue, cyan, green, yellow, red, magenta, pink).
   4. You may want to fine tune the values for each level to match the minimum value of the classed post map classes exactly.
   5. Click Save and saved the levels to an LVL file (ContourLevels.lvl) and click OK.
   6. Check the check box next to Fill contours.
4. Load the LVL file for each of the four remaining contour layers and fill the contours. For example:
   1. Select one of the other contour layers and click on the Levels tab in the Property Manager.
   2. Change the Level method to Advanced, and click Edit Levels.
   3. Click Load, select the LVL file and click Open and OK. 
   4. Check the check box next to Fill contours.
5. All the methods seem relatively accurate to the data, so that is good. But I do not like how Triangulation with Linear Interpolation interpolated the data in between the flight lines. When I zoom in, it appears quite jagged. So I remove that option. That leaves four gridding methods.

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.

1. Click Grid | Data, select NV_1152.dat and click Open.
2. In the Grid Data dialog:
   1. Set the XYZ columns to C, D and K, respectively.
   2. Select Inverse Distance to a Power gridding method
   3. Change the Spacing in the X direction to 0.0001 (the distance between points in the X direction, as measured above).
   4. Change the Spacing in the Y direction to 0.001 (about 1/5 the distance between flight lines).
   5. Check Blank grid outside convex hull of data so that the areas outside the flight lines are not interpolated.
   6. Change the Output Grid File name to include the gridding method (e.g. InverseDistance_HiRes.grd).
   7. Click OK. The grid is created.
3. Now let’s update the map with the new grid. Select the contour layer for Inverse Distance, and in the Property Manager, on the General page, click the Open Grid button, select the new grid file and click Open. The contours are updated.
4. Repeat steps 1-3 for Kriging, Natural Neighbor, and Minimum Curvature.
5. Once all the maps are updated, I can see that gridding with Minimum Curvature at a higher resolution method really smeared the contours out. I don’t like that. Setting an anisotropy value (an advanced option for that gridding method) didn’t seem to help, so I eliminated that method. That left 3 potential gridding methods: Inverse Distance to a Power, Kriging, and Natural Neighbor. 

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.

1. Click Grid | Residuals.
2. Select the high resolution Kriging grid and click Open.
3. Select NV_1152.dat and click Open.
4. Select columns C, D and K as XYZ and store the residuals in column L.
5. Click OK. The worksheet opens with the data and the new Residuals column.
6. Click File | Save to save the worksheet and then File | Close to close it.
7. Click Grid | Residuals again.
8. Select the high resolution Inverse Distance grid and click Open.
9. Select NV_1152.dat and click Open.
10. Select columns C, D and K as XYZ and store the residuals in column M.
11. Click OK. The worksheet opens with the data and the new Residuals column.
12. Click File | Save to save the worksheet and then File | Close to close it.
13. Click Grid | Residuals a last time.
14. Select the high resolution Natural Neighbor grid and click Open.
15. Select NV_1152.dat and click Open.
16. Select columns C, D and K as XYZ and store the residuals in column N.
17. Click OK. Now all three residuals are there in the worksheet.
18. Let’s calculate the square of the residuals.
    1. Select column L by clicking on the column header.
    2. Click Data | Transform, enter the function L=L*L
    3. Click OK.
    4. Click Data | Statistics, make sure Sum is checked and click OK. The sum of the square of the residuals (for the Kriging grid) is displayed. Note this value and click Close.
    5. 
19. Repeat Step 17 for columns M (Inverse Distance) and N (Natural Neighbor).
20. The results show that Kriging has the lowest residuals and so is the most accurate relative to the original data. You can calculate the residuals for the original Kriging grid created by the script with the default settings to compare if decreasing the grid spacing (increasing the resolution) has increased the accuracy, and yes it has.

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:

1. The Help (Help | Contents): On the Contents page, review the topics under the Gridding | Gridding Methods book
2. Webinar Gridding in Surfer
3. Knowledge Base articles:
   1. What types of gridding methods are available?
   2. What gridding method is best for my data file?
   3. What advanced options are available for each of the different gridding methods?
4. Newsletter articles:
   1. What gridding method should I use? Gridding for non-geostatisticians
   2. A Basic Understanding of Surfer Gridding Methods – Part 1
   3. A Basic Understanding of Surfer Gridding Methods – Part 2
5. Paper by Chin-Shung Yang, Szu-Pyng Kao, Fen-Bin Lee, Pen-Shan Hung: TWELVE DIFFERENT INTERPOLATION METHODS: A CASE STUDY OF SURFER 8.0