HMCA-Contour: A Visual Basic Program Based on Surfer Automation for Soil Heavy Metal Spatial Distribution and Contamination Assessment Mapping

: Soil heavy metal contamination has become a major concern in many parts of the world. Grasping the spatial distribution patterns of heavy metals and evaluating soil heavy metal contamination are of great signiﬁcance to policy makers for land use planning and pollution identi-ﬁcation. However, these processes usually involve repetitive calculation and plotting, which are time-consuming and laborious, especially in a geochemical survey. In this paper, HMCA-Contour, which is a program written in Visual Basic based on Surfer 13.0, has been designed to simplify the pollution indices recalculation and plotting routines. This software is comprised of three functions, i.e., data processing, plotting, and template setting. With user-friendly interfaces, the user can easily and efﬁciently calculate pollution indices and generate batch publication-quality maps of heavy metal spatial distribution and pollution levels. To demonstrate the advantage of HMCA-Contour, a case study from the Gaerqin Mine area is used in this paper. HMCA-Contour, as a brand new software, is promising and will become a useful tool for environmental management.


Introduction
Soil is usually regarded as both a sink and source for heavy metals (HMs) and other pollutants [1]. An excess level of heavy metals in soil poses a potential threat to the health of plants, animals, and even humans via the food chain because of their persistence, high toxicity, non-biodegradable, and accumulative properties [2][3][4][5][6]. Consequently, soil heavy metal pollution has become a major concern in environmental management. Numerous studies have focused on HMs' spatial distribution and pollution assessment during the last few decades [2,[7][8][9][10][11]. In these studies, many pollution indices have been applied to evaluate soil contamination assessment. Among these indices, the geo-accumulation index [12], pollution load index [13], enrichment factor index [14], risk index [15], and Nemero index are the most commonly used methods. Additionally, the contour map is also a powerful tool for visualizing the contamination degree and grasping the spatial characteristics of HMs' concentration.
Generally, researchers usually rely on interpolation methods, including kriging [16,17] and inverse distance [18,19] approaches, to estimate the HM concentrations and pollution indices of the unsampled area and therefore construct continuous surface distribution maps [20,21]. In most existing studies, Surfer TM (Golden Software, Inc., Golden, CO, USA) and ArcGIS software (ESRI, Sacramento, CA, USA) can provide a way to realize interpolation algorithms in the form of a contour map [22][23][24][25]. Compared with ArcGIS, Surfer software is more efficient, being faster in data processing and plotting. However, when using this software, several time-consuming and labor-intensive problems may still arise due to a flood of geochemical data caused by multiple elements or too many sampling sites. For example, to assess soil contamination, researchers attempt to choose a proper method to calculate the pollution index several times in Microsoft Excel TM , before generating contour maps via Surfer software. Then, they repeatedly plot contour maps based on the above calculation results or spatial distribution contour maps based on sampling data sets for each single heavy metal. Therefore, it is of great significance to eliminate the above routines and tedious manual operations. Gong developed the software EGAPI, which was built with the above pollution indices methods, representing one of the few attempts to improve the efficiency of pollution index calculation [26]. Nevertheless, this is far from enough, as there are still other problems that need to be solved, such as the need for a user-friendly interface, the lack of batch plotting HM contamination assessment or spatial distribution contour maps, and the requirement of the batch output for publication-quality contour maps.
To solve the problems outlined above and develop a more scientific, faster in data handling and contour map drawing software, a combination of Visual Basic (VB) programming language and Surfer automation technology would be a good solution. The Surfer automation technique allows scripter or other high-level programming languages, such as VB or C++, to call surfer's drawing function and is widely used in soil moisture [27], exploration geochemistry [28], meteorology [29], and hydrology [30]. Similarly, an integrated software HMCA-Contour written in VB and based on Surfer automation was designed to eliminate routine and tedious operations involving datasets of HMs in soil topics. The main functions of HMCA-Contour include the following: (1) Batch calculating pollution indices for HMs; (2) batch plotting maps of the spatial distribution and soil contamination assessment of HMs; and (3) template setting for contour maps. This paper presents the details of HMCA-Contour through a case study.

Methodology
Several common contamination assessment indexes of soil HMs, including the geoaccumulation index (I geo ), pollution load index (PLI), enrichment factor (EF), potential ecological risk index (RI), single pollution index (SI), and Nemero index (NI), were selected and integrated into HMCA-Contour software.

Geo-Accumulation Index
The geo-accumulation index (I geo ) was put forward by Müller [12]. Originally used in river bottom sediments, it can also be applied to the assessment of heavy metal contamination in soil [31]. I geo is computed by the following formula: where C n is the measured content of element n in soil, B n is the geochemical background value of element n in soil, and the constant 1.5 is the background matrix correction factor due to lithogenic effects in different regions. According to Müller [12], the geo-accumulation index consists of 7 grades or classes (Table 1), whereby the highest class 6 is an open class comprising all values of the index higher than class 5 and may reflect a 100-fold enrichment above the geochemical background value.

Pollution Load Index
The pollution load index (PLI), proposed by Tomlinson et al. [13], was employed to determine the magnitude of heavy metal contamination in the sediment. PLI was determined as the nth root of the product of the n contamination factor (CF). It is expressed as follows: where n is the number of metals and CF is the heavy metal concentration in the soil/background value of the metal. Interpretation PLI > 1 indicates polluted, whereas PLI < 1 indicates no pollution.

Enrichment Factor Index
The enrichment factor (EF) was first proposed by Zoller et al. [32] based on a study of the enrichment degree and source of elements in the atmosphere over Antarctica. Originally applied to the atmosphere and seawater, EF was gradually introduced to environment media in continental environments, such as soils [14], stream sediments [33], and river suspended solids [34] and it has been widely used in environmental geochemical heavy metal pollution assessment. The EF is defined mathematically as where the concentration of the heavy element i is C i and C re f refers to the concentration of the reference element. The most common reference elements are Fe, Al, Ti, Mn, and Sc [34][35][36][37][38][39]. According to Sutherland [35], EF is divided into five categories: Minimal pollution (EF < 2); moderate pollution (EF = 2-5); significant pollution (EF = 5-20); strong pollution (EF = 20-40); and extreme pollution (EF > 40).

Ecological Risk Index
The potential ecological risk index (RI) method is a sediment evaluation method proposed by Swedish scientist Hakanson [15] in 1980. It is a relatively fast, simple, and standard method for evaluating the sediment pollution degree and potential ecological risk index from the perspective of sedimentology.
The potential ecological risk index (RI) was suggested by Hakanson [15] and conducted to evaluate the impact of heavy metals in soil or sediment [40][41][42]. Three aspects were considered in this method: The toxicity of various heavy metals; the synergistic action of multiple elements; and the sensitivity of various heavy metal contamination in soils or sediments. Three Equations (4)-(6) for computing RI are as follows: where T i r refers to the toxic-response factor for a given element i, with values of 10, 2, 5, 30, 5, 1, and 40 for As, Cr, Cu, Cd, Pb, Zn, and Hg, respectively; E i f represents an ecological risk factor for heavy metal i; C i f , C i s , and B i n are the contamination factor, the concentration in the sediment, and the background reference value for element i, respectively; and the potential ecological index (RI) is defined as the sum of risk factors.
According to Hakanson [15], the Er and RI can be divided into four grades and these are described in Table 2. High potential risk 360 ≤ RI < 720 Very strong potential ≥320 Significantly very high ≥720 Highly strong potential

Single and Nemero Pollution Index
The single pollution index (P i ) quantifies the pollution of one single heavy metal [43]. It was defined as the ratio between C i and S i : where C i identifies the concentration of each heavy metal in each sample and S i refers to the soil standard of each heavy metal. The Nemerow pollution index (P N ) provided by Nemerow [44] was initially employed to water quality assessment and gradually extended to soil contamination [45][46][47]. This method takes both the mean and maximum of single pollution levels of all the evaluation factors into consideration. It was adopted and calculated as where P iave and P imax represent the maximum value and average value of the single factor pollution index in each sample, respectively. P i and P N were respectively divided into four and five grades [48,49] and their criteria are listed in Table 3. Table 3. Grading standard for the single pollution index (P i ) and Nemerow pollution index (P N ).

P i Appraisal Result
Seriously polluted domain

Software Operation Environment Requirements
HMCA-Contour is temporarily only provided for mainstream Microsoft Windows (since 7), and it requires Excel 2007 or later and Surfer 13.0 to be installed on the user's computer, which enables Visual Basic to invoke the automation method provided by Surfer software (version 13.0). The download link can provide access to the installation file from the website (https://github.com/HMCAHome/HMCA-Contour (accessed on 10 October 2020)). After downloading the installation file, the user can then run HMCA-Contour setup.exe and follow the installation directions on the screen. Once the installation October 2020)). After downloading the installation file, the user can then run HMCA-Contour setup.exe and follow the installation directions on the screen. Once the installation procedure is completed, the user can run executable file HMCA-Contour. The flow chart of the operation of this software can be seen in Figure 1.

Data Input
HMCA-Contour accepts two types of input files, i.e., .xlsx or .xls files. With the sheet named "Sheet1", the first two columns of the template input file must be X and Y coordinates, respectively, and followed by a series of elements ( Figure 2) that can be in any order. The user can enter the "Data Processing" interface by clicking the "Data processing" tab in the toolbar of the HMCA-Contour main interface ( Figure 3a). In this interface, three subpanels are given, including "File Location", "Data Extraction", and "Data Processing".
In the "File Location" subframe, the user loads the data file by clicking the "Open" button to choose the input file. After selection of the input file, this program automatically extracts the first row data in the left box of the "Data Extraction" frame, and the user is capable of choosing target items by clicking direction buttons (buttons ">", "<", ">>", and "<<" represent add, reduce, select all, and deselect all, respectively) ( Figure 3b).

Data Input
HMCA-Contour accepts two types of input files, i.e., .xlsx or .xls files. With the sheet named "Sheet1", the first two columns of the template input file must be X and Y coordinates, respectively, and followed by a series of elements ( Figure 2) that can be in any order. The user can enter the "Data Processing" interface by clicking the "Data processing" tab in the toolbar of the HMCA-Contour main interface ( Figure 3a). In this interface, three subpanels are given, including "File Location", "Data Extraction", and "Data Processing".  In the "File Location" subframe, the user loads the data file by clicking the "Open" button to choose the input file. After selection of the input file, this program automatically extracts the first row data in the left box of the "Data Extraction" frame, and the user is capable of choosing target items by clicking direction buttons (buttons ">", "<", ">>", and "<<" represent add, reduce, select all, and deselect all, respectively) ( Figure 3b).

Data Processing
In the "Data Processing" module, six HMCA methods, as mentioned above, are provided in this software (Figure 3b). The user can select the target method to work with data. It is noted that EF requires additional parameters, including the reference element (Cref element) and background value (Cref Background) (Figure 3b). Once the evaluation method is selected, the user can click the "Calculate" button to start calculating and is prompted to input the background value for each item according to the default format in a message box (Figure 3(C1)). After inputting the background value in the default form for each element, the user clicks the "OK" button, and the whole data handling process is then completed in seconds (Figure 3(C2)). The resulting filename is in the form of "source file name"_"selected method".xls and saved in the "Output File" path directory, which is the same as the "Source File" path in default. It is noted that the units of heavy metal concentrations are required to be uniform (mg/kg).

Data Input
In order to enter the plotting interface (Figure 3d), the user simply clicks the "Next Step" button in the "Data Processing" interface or the plotting tab in the main interface of this software. In this function, MS Excel XLSX, XLS, and Comma Separated Value (CSV) files can be used as original data files for HMCA-Contour. Like "Data Extraction" in "Data Processing" capability, the program can extract the data automatically from the original file in the left textbox when the user clicks the "open" button. At the same time, the grid geometry parameters (Minimum, Maximum, Spacing, and Nodes) are calculated automatically. Additionally, the default options for plotting are loaded, for example, "File name", "Output Path", "Image Type", and "Image Quality" in the output filename description frame, and the "Font Face", "Font Size", "Symbol Size", "Symbol Type" in the "Postmap" frame. The program provides four kinds of image quality, from low to high.

Data Processing
In the "Data Processing" module, six HMCA methods, as mentioned above, are provided in this software (Figure 3b). The user can select the target method to work with data. It is noted that EF requires additional parameters, including the reference element (Cref element) and background value (Cref Background) (Figure 3b). Once the evaluation method is selected, the user can click the "Calculate" button to start calculating and is prompted to input the background value for each item according to the default format in a message box (Figure 3(C1)). After inputting the background value in the default form for each element, the user clicks the "OK" button, and the whole data handling process is then completed in seconds (Figure 3(C2)). The resulting filename is in the form of "source file name"_"selected method".xls and saved in the "Output File" path directory, which is the same as the "Source File" path in default. It is noted that the units of heavy metal concentrations are required to be uniform (mg/kg).

Data Input
In order to enter the plotting interface (Figure 3d), the user simply clicks the "Next Step" button in the "Data Processing" interface or the plotting tab in the main interface of this software. In this function, MS Excel XLSX, XLS, and Comma Separated Value (CSV) files can be used as original data files for HMCA-Contour. Like "Data Extraction" in "Data Processing" capability, the program can extract the data automatically from the original file in the left textbox when the user clicks the "open" button. At the same time, the grid geometry parameters (Minimum, Maximum, Spacing, and Nodes) are calculated automatically. Additionally, the default options for plotting are loaded, for example, "File name", "Output Path", "Image Type", and "Image Quality" in the output filename description frame, and the "Font Face", "Font Size", "Symbol Size", "Symbol Type" in the "Postmap" frame. The program provides four kinds of image quality, from low to After reading data from the source file, the user may choose the target items from the left text box to the right text box in the "Choosing Items" subpanel. The buttons for selecting target elements are the same as those mentioned above (Figure 3d). It is noted that two items (x, y) must be left in the left text box. The user selects a gridding algorithm from the drop-down menu of "Gridding Method". Additionally, two other kinds of data files are available in this function-the boundary data (.BLN) and the geographic data (.dat, .txt, .xls, and .xlsx)-with the former one being employed for boundary determination and the latter one being employed for the geographic location. Both "blank" and "post map" files are loaded by clicking the "open" button. In the post map frame, the symbol size and type, and the symbol label font face, can be set with the corresponding drop button. Plotting capabilities are based on Surfer automation. The user does not need to worry about the running process, whose flow chart can be seen in Figure 4.

Plotting
After reading data from the source file, the user may choose the target items from t left text box to the right text box in the "Choosing Items" subpanel. The buttons for sele ing target elements are the same as those mentioned above (Figure 3d). It is noted th two items (x, y) must be left in the left text box. The user selects a gridding algorithm fro the drop-down menu of "Gridding Method". Additionally, two other kinds of data fi are available in this function-the boundary data (.BLN) and the geographic data (.d .txt, .xls, and .xlsx)-with the former one being employed for boundary determinati and the latter one being employed for the geographic location. Both "blank" and "p map" files are loaded by clicking the "open" button. In the post map frame, the symb size and type, and the symbol label font face, can be set with the corresponding drop b ton. Plotting capabilities are based on Surfer automation. The user does not need to wor about the running process, whose flow chart can be seen in Figure 4.

Gridding Method Descriptions
Surfer software provides twelve algorithms built in this program [50], including Inverse Distance to a power, kriging, Minimum Curvature, Modified Shepard's Method, Natural Neighbor, Nearest Neighbor, Polynomial Regression, Radial Basis Function, Triangulation with Linear Interpolation, Moving Average, Data Metrics, and Local Polynomial. According to the different characteristics of various types of data and different requirements, the user can scientifically choose the right interpolation through the HMCA-Contour interface to call Surfer software to grid discrete data without worrying about the gridding process. Since these interpolation algorithms are default in Surfer, the contour maps produced by HMCA-Contour and Surfer software can lack model accuracy in comparison to Arcgis under certain circumstances.

Template Setting
The user may enter the plotting setting interface by either clicking the "Next" button in the "Data processing" interface or clicking the function "Template Setting" tab in the main interface of HMCA-Contour. The Template Setting interface has five tabs: "Color scale"; "Axis Ticks"; "Axis Title"; "Axis Labels"; and "Map Title & Scale bar & North Arrow" ( Figure 5). This function is designed to build a publication-ready template for a series of thematic maps. The details of the tabs are given below.

Color Scale
The program provides three group parameters to customize a color scale (Figure 5a,b). Firstly, the user can type a title for the color scale that consists of a prefix, element, and suffix. The "prefix" and "suffix" text boxes can be activated by their front check boxes. When activated, the "title" text can be formed by elements and activated text boxes. Secondly, two forms of color scale bars are offered: The simple one is a system-based preset, which provides several kinds of color scale bars in the "Preset" drop-down menu, whilst the advanced one requires the user to click the "more" button ( Figure 5a). As soon as the "more" button is clicked, there is a pop window through which the user is capable of adjusting a desired color scale example based on "nums", "zlevel", and "lcolor" parameters ( Figure 5b). The user then clicks the "ok" or "return" button. Note that the user must choose one method to show the color scale bar. Lastly, the color scale bar label is set as the axis label.

Axis Ticks
Four kinds of axis ticks are provided, including bottom, left, right, and top ticks. Each of them has the same parameters (Figure 5c). Take the bottom axis tick as an example, where the user can set parameters in terms of five aspects: (1) A major or minor tick type drop-down menu that enables the user to choose the desired tick type from a list containing "none", "in", "out", and "both"; (2) a major or minor tick length spin box, which allows the user to use an up or down arrow button to respectively increase or decrease the number value of the tick length in the text field by a 0.01 increment; (3) a major or minor gird line drop-down menu that offers two choices ("True and False") for the user to decide whether or not to display it; (4) a minor per major spin button that enables the user to set the number of minor ticks between each two major ticks; and (5) an axis line width spin button makes it possible for the user to adjust the axis line width by clicking the up or down arrow button. The same approach of adjusting parameters can be applied for the remaining kinds of axis ticks.

Axis Titles
The user is capable of customizing the axis title properties (Figure 5d). For example, if the user wants to add a bottom axis title, the user can type text in the text box, and then change the angle of the title by clicking the up or down spin button. Additionally, the user can set the title text font, size, vertical and horizontal distance to the axis, bold, and italic maps produced by HMCA-Contour and Surfer software can lack model accuracy in comparison to Arcgis under certain circumstances.

Template Setting
The user may enter the plotting setting interface by either clicking the "Next" button in the "Data processing" interface or clicking the function "Template Setting" tab in the main interface of HMCA-Contour. The Template Setting interface has five tabs: "Color scale"; "Axis Ticks"; "Axis Title"; "Axis Labels"; and "Map Title & Scale bar & North Arrow" (Figure 5). This function is designed to build a publication-ready template for a series of thematic maps. The details of the tabs are given below.

Color Scale
The program provides three group parameters to customize a color scale (Figure 5a,b). Firstly, the user can type a title for the color scale that consists of a prefix, element, and suffix. The "prefix" and "suffix" text boxes can be activated by their front check boxes. When activated, the "title" text can be formed by elements and activated text boxes. Secondly, two forms of color scale bars are offered: The simple one is a system-based preset, which pro-

Axis Labels
Four check boxes allow the user to decide which axis labels need to be displayed on a map (Figure 5e). Once a check box is clicked, the user can further adjust the corresponding axis label parameters. The user clicks the up or down spin button to set a proper distance value between the axis line and labels. Moreover, by clicking the drop-down buttons, the font face, size, bold, and italic and angle of axis label can be adjusted via their drop list.

Map Title, Scale Bar, and North Arrow
Four sections are divided (Figure 5f). Among them, the method of setting parameters in the map title section occurs in the same way as in the color scale title section. With respect to the scale bar component, the user needs to adjust the scale bar label and appearance separately. A scale bar is a tool that helps the user with the reading of this scale when a scale rule is not available, whose style can be approached from three aspects: The spacing between cycles (Cyclespacing); the increment of the label value (Increment); and the number of cycles (Numcycles). All of these parameters can be set from their corresponding drop-down lists. Moreover, the user is able to change the angle of the scale bar. The north arrow section enables a north arrow to be set up on the map display. The user can use the "symbol type" drop-down list to choose the style, and then specify the location of its center position from the "X position" and "Y position" drop-down lists.

Study Area and Data
To test the applicability of HMCA-Contour, we examined the soil geochemical data taken from the Gerqin Mine area [51]. This mining area (83 • 23 00 E~83 • 27 00 E, 32 • 47 00 N~32 • 50 00 N) is located in Wuma Township, Gaize County, Tibet. It is located about 100 km northwest of Gaize county, with a height of 4700-5500 m, relative height of 100-700 m, and area of 1700 km 2 . The climate of this area is cold and dry, with a long sunshine time, a large temperature difference between the day and night, and annual precipitation of about 190 mm. There is a large amount of evaporation in the mining area.
The soil in the study area is mainly composed of Quaternary alluvial diluvium and weathering products of the bedrock parent material. From east to west, there are eight mineral deposits (occurrences), including Bolong, Duobuza, Rongna, Naruo, Sena, Tiegelong, Saijiao, and Gaerqin. In this study, we take the southeast of the Duolong mining area, referred to as the Gaerqin mining area, as the study area. Currently, the Duolong area is in the stage of pre-mining, so it is of great significance to know the extent of HM concentration distribution characteristics and assess soil HM contamination.
In this study, a total of 835 surface topsoil samples were collected, with a density of about four samples per km 2 . In the vicinity of the pre-set sampling point, radiating 10-30 m to the surrounding area, 4-5 points were combined equally to form a mixed sample, and attempts were made to ensure the consistency of the sampling depth and soil type. Subsequently, the concentrations of Cu, Pb, Zn, Cd, Cr, Hg, and As in the samples were measured by using X-ray fluorescence (XRF) and Atomic absorption spectrometry (AAS). More detailed information on sample preparation and analysis can be found in [49]. In this study, Pb, Zn, Cr, and Hg were selected to demonstrate the advantage of HMCA-Contour, and a statistic description of their data is shown in Table 4.

Pollution Indices Calculation
As mentioned above, original data are arranged in the form of Figure 2, and the file format is .xlsx. Through the data processing function window, the user loads the test data and chooses the target elements to be processed. After selecting the I geo method, the user clicks the "calculate" button and the I geo calculation message box then immediately pops up, requiring the user to input the corresponding background values (500, 500, 300, and 1.5 mg/kg for Pb, Zn, Cr, and Hg, respectively) for each selected heavy metal in the default form of "{500, 500, 300, 1.5}". After inputting the background values and clicking the "OK" button, the calculated I geo result with its statistic description file (Table 4) is then saved in the output path directory.

Spatial Distribution of Heavy Metal Concentrations
Spatial distribution maps of Pb, Zn, Cr, and Hg were created from the original geochemical data by means of the "Inverse Distance to a power" method in this software. In the plotting function window, the user loads the calculated I geo file and chooses Pb, Zn, Cr, and Hg as the target elements. Additionally, the user sets the parameters of "Axis Ticks", "Axis Labels", and "Axis Titles". The next step is to click the "OK" button to plot, and the bottom progress bar shows specific information about the batch plot. After clicking the prefix and suffix buttons of both the map title and color scale title to activate text boxes, the user enters text into the corresponding text boxes.
The geochemical maps of soil heavy metals produced by HMCA-Contour are shown in Figure 6. The spatial distribution patterns of Pb and Zn concentration hotspots are almost the same in the center of the maps, indicating that they probably have similar geochemical behaviors. In comparison, the Cr concentration is mainly concentrated in the northeast of the study area. Ore and rock weathering may be responsible for the spatial distribution patterns of soil heavy metals. According to the geological facts in the study area, there are ultramafic rock units in the northeast part of the study area, which can be easily eroded and weathered into soil, resulting in the accumulation of Cr in the topsoil. and the bottom progress bar shows specific information about the batch plot. After clicking the prefix and suffix buttons of both the map title and color scale title to activate text boxes, the user enters text into the corresponding text boxes. The geochemical maps of soil heavy metals produced by HMCA-Contour are shown in Figure 6. The spatial distribution patterns of Pb and Zn concentration hotspots are almost the same in the center of the maps, indicating that they probably have similar geochemical behaviors. In comparison, the Cr concentration is mainly concentrated in the northeast of the study area. Ore and rock weathering may be responsible for the spatial distribution patterns of soil heavy metals. According to the geological facts in the study area, there are ultramafic rock units in the northeast part of the study area, which can be easily eroded and weathered into soil, resulting in the accumulation of Cr in the topsoil.

Heavy Metal Contamination Assessment Maps
Heavy metal contamination assessment maps were generated via ( Figure 7). When the value is higher than 0, soils are polluted. Therefore, there is no soil pollution caused by Pb, Zn, and Cd in this study area. However, based on the result map of Cr , the most contaminated areas are located in the northeast part of the Duolong mining area, where the ranged from 0 to 5, indicating that the soil quality ranges from a moderately contaminated level to extremely contaminated level. Therefore, further studies are required to investigate and analyze the contamination source in this area. Knowing the extent of pollution and pollution sources is an essential and important step towards the prevention and control of heavy metal contamination.

Heavy Metal Contamination Assessment Maps
Heavy metal contamination assessment maps were generated via I geo (Figure 7). When the I geo value is higher than 0, soils are polluted. Therefore, there is no soil pollution caused by Pb, Zn, and Cd in this study area. However, based on the result map of Cr I geo , the most contaminated areas are located in the northeast part of the Duolong mining area, where the I geo ranged from 0 to 5, indicating that the soil quality ranges from a moderately contaminated level to extremely contaminated level. Therefore, further studies are required to investigate and analyze the contamination source in this area. Knowing the extent of pollution and pollution sources is an essential and important step towards the prevention and control of heavy metal contamination.

Conclusions
Contour maps are vital for allowing environmental geochemists to visualize the spatial distribution of soil heavy metals and contamination degree. A user-friendly, handy, and high-efficiency software, named HMCA-Contour, was designed based on Visual Basic 6.0 language and Surfer automation technology and a case study was introduced in this paper. HMCA-Contour can automatically calculate the soil pollution indices with their statistic description and plot the spatial distribution of soil heavy metal maps and soil contamination assessment maps. Meanwhile, several details of generated maps can be set according to the user's requirements in the template setting function. Previous repetitive calculation and contour map drawing routines can be shortened by HMCA-Contour software. The real case study tells us that it only takes within a second to calculate the pollution indices and 3~5 s to draw a contour map for a single heavy metal. After a very small manual intervention, the user can easily achieve the expected results. The program can be widely used in pollution indices calculation and plotting maps of soil HMs' spatial distribution and soil pollution levels. It also has implications for the development and improvement of heavy metal contamination assessment software.
Author Contributions: Q.L. designed the study, developed the software, and drafted the manuscript; G.L. and W.C. revised and refined the manuscript; G.C. contributed to the data calculation and plotting. All authors have read and agreed to the published version of the manuscript.

Acknowledgments:
The authors wish to thank Zhifeng Zhang, Zebin Zhang, and Ruiliang Wang for testing this software. We also thank the anonymous reviewers and the associate editor for their comments, which improved the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Conclusions
Contour maps are vital for allowing environmental geochemists to visualize the spatial distribution of soil heavy metals and contamination degree. A user-friendly, handy, and high-efficiency software, named HMCA-Contour, was designed based on Visual Basic 6.0 language and Surfer automation technology and a case study was introduced in this paper. HMCA-Contour can automatically calculate the soil pollution indices with their statistic description and plot the spatial distribution of soil heavy metal maps and soil contamination assessment maps. Meanwhile, several details of generated maps can be set according to the user's requirements in the template setting function. Previous repetitive calculation and contour map drawing routines can be shortened by HMCA-Contour software. The real case study tells us that it only takes within a second to calculate the pollution indices and 3~5 s to draw a contour map for a single heavy metal. After a very small manual intervention, the user can easily achieve the expected results. The program can be widely used in pollution indices calculation and plotting maps of soil HMs' spatial distribution and soil pollution levels. It also has implications for the development and improvement of heavy metal contamination assessment software.