Computerized Contour Mapping
By Darrel Dunn, Ph.D., PG, Hydrogeologist - Geologist
(View Résumé 🔳)
Purpose - Describe computerized contouring methods
The purpose of this web page is to describe two computerized contouring methods that are compatible with the FloPy module that develops MODFLOW input files. The contouring methods could also be used to produce maps for petroleum exploration and production applications, or the analysis of any areally distributed data. The contouring methods are (1) gridded averaging, and (2) Matplotlib linear interpolation on a triangular grid (matlotlib.tri.LinearTriInterpolator). This page is technical and assumes the reader is familiar with Python and MODFLOW.
Gridded averaging contouring method
The gridded averaging contouring method has not been described previously. The method uses a set of Python programs that I developed. The programs produce an array of values for cells on a gridded map. Known values (constants) are assigned to grid cells according to their locations, and values are calculated for the remaining cells. The calculated values are the average of the values in the immediately adjacent cells. If the grid used in a MODFLOW model is also used for the contouring, then a Python NumPy array is produced that can be used directly as input to FloPy. One of the programs produces a grid plot of the constant and calculated values and another produces a contour map. The gridded averaging method is capable of introducing directional trend to the contours by using weighted averaging. Introducing trend might be useful in some geologic mapping where trends are known to exist. Since NumPy arrays are involved, there is also an opportunity to manipulate the arrays before contouring. Gridded averaging could be extended to three-dimensional interpolation.
The plotting and contouring method involves the following steps:
Develop a map showing the locations of the constant values.
Overlay the map with a grid that may be the same as a MODFLOW model grid.
Enter the grid locations of the constants and their values into fdcontour_in.py. The fdcontour_in.py script is shown in Attachment 1. The locations are entered in an array that is analogous to the IBOUND array of MODFLOW.
Run final_plot.py with Python3. This run produces the grid plot as a png file. The final_plot.py script is shown in Attachment 2. This step may be skipped if a plot is not wanted.
Run array_contour.py with Python3. This run produces a Numpy array and a contour map of the constant and calculated values as a png file. The array_contour.py script is shown in Attachment 3.
The scripts final_plot.py and array_contour.py each import variables from fdcontour.py. The fdcontour.py script is shown in Attachment 4. The fdcontour.py program calculates the values in the cells that do not contain constants as the average value of the immediately adjacent cells. This is accomplished by solving the set of equations for average cell value simultaneously using successive overrelaxation (SOR).
Figure 1 shows the NumPy array produced for a simple artificial example of the method. In this example the array shape is 5X5. A constant value of 10 is assigned to the center cell, and constant values of 5 are assigned to corner cells and cells in the middle of each side of the grid. Figure 2 shows the plot of the array, and Figure 3 shows a contour map of the array. Inspection of Figure 1 shows that computed values in interior cells are indeed the average of the four surrounding cells. Computed values of cells along the sides of the grid are the average of the three adjacent cells. It is not necessary to put constants in the corner cells. When corner cells are computed values, they are the average of the two adjacent cells. Figure 2 is a plot of the values in Figure 1. Figure 3 is a contour map of the array in Figure 1 produced by the contour function of Matplotlib.
Figure 1. NumPy array.
Figure 2. Grid plot of symmetrical array.
Figure 3. Contour map of symmetrical array produced by gridded averaging.
Figure 4 shows a Matplotlib contour map of the array produced by moving the 10 constant from the center cell to the cell located left of the center, so that the constant locations are not symmetrical.
Figure 4. Contour map of asymmetrical array produced by gridded averaging.
Linear interpolation on triangular grid contouring method
I performed linear interpolation on a triangular grid with Matplotlib using tri.Triangulation and tri.LinearTriInterpolator. The Python script is shown in Attachment 5. I applied this interpolation method to the same constants as in the symmetrical example above. It produced an array similar to that produced by the averaging method. Figure 5 shows the contour map produced from the array. Figure 5 is analogous to Figure 3.
Figure 5. Contour map of symmetrical array produced by Matplotlib triangular linear interpolation.
Figure 6 shows the result when the triangular grid method was applied to the asymmetrical data. Figure 6 is analogous to Figure 4.
Figure 6. Contour map of asymmetrical array produced by Matplotlib triangular linear interpolation.
Application of gridded averaging contouring method
Figure 7 illustrates the use of the gridded averaging contouring method for analyzing real data. It is a map of bedrock elevation beneath surficial deposits. It shows a bedrock surface sloping toward a stream that crosses the southwest corner of the map. The map was constructed by putting bedrock tops reported in the boreholes shown on the map into fdcontour_in.py. I also added the constraint that the bedrock is at least 20 feet below the land surface in areas where there is no borehole data. The grid cells are 1000 X 1000 feet, as shown in Figure 8.
Figure 7. Bedrock elevation contour map.
Figure 8. Bedrock elevation gridded color plot.
I also made a gridded map of the land surface elevation. I did this by putting critical elevations into fdcontour_in.py and running array_contour.py. The result is shown in Figure 9. Some minor offset is caused by the gridding.
Figure 9. Land surface elevation contour map.
Since the numpy arrays for the bedrock elevation and the land surface elevation are the same shape, bedrock may be subtracted from land surface to get an isopach map of the thickness of the surficial deposits (Figure 10).
Figure 10. Isopach map of surficial deposits.
Since these gridded maps are numpy arrays, they are compatible with FloPy development of MODFLOW input files.
Attachment 1 - fdcontour_in.py
Attachment 2 - final_plot.py
Attachment 3 - array_contour.py
Attachment 4 - fdcontour.py
Attachment 5 - pycontour.py
Posted January 17, 2019.
Revised April 8, 2021.