GSJ Open-File Report, no. 449 Jannuary, 2007

Software system for Aeromagnetic data processing, Grid data manipulation, and Reduction and quantitative interpretation of magnetic anomaly data

Tadashi Nakatsuka (Institute of Geology and Geoinformation, GSJ/AIST)

1. Introduction

I have been engaged in the research on aeromagnetic survey, and developed various kinds of programs for data processing, analysis and interpretation. The outlines of them have already been reported on the scientific journals and at the scientific meetings [Nakatsuka and Okuma, 2005a, 2005b, 2006a, 2006b; etc.]. This report opens their source codes to the public.

In this report, I deal with three groups of programs for the aeromagnetic survey, as follows.
    1) DPAM: programs for aeromagnetic survey data processing,
    2) GDMP: programs for grid data manipulation,
    3) ANAM: programs for analysis of magnetic anomaly data.
This software system is in a style of program libraries, consisting of many individual programs of rather simple functions. Then the actual process to some target/objective would require a series of program executions. The data formats used in this system are customized for our own original style. However, as for its design, all data are handled in forms of ASCII text, so that the user can easily convert data to/from another style of formats.

2. DPAM programs

DPAM group consists of programs for aeromagnetic (AM) survey data processing: log2asc, xldam, xldpn, xldhg, despike, dvcorr, ecomp, fcomp, ggrid, ggrids, pframe, pltrk, pchkdv, pchkmag, pchkres, pchkcomp, xslin and xslina.

The program alog2asc converts AM binary raw data (specialized for our survey system) into ASCII file format, xldam and xldpn generates line data (DPAM line data) file from AM ASCII raw data, making use of Real-time GPS and PostNav processed DGPS position data, respectively, and xldhg generates similar line data (HGAM line data) file from helicopter-borne gradiometer AM acquisition data. In addition, xldam, xldpn and xldhg also calculate the IGRF residuals, using the subprogram "igrf" by Nakatsuka [2005].

The program despike eliminates spike noises in magnetic field values (total magnetic field and IGRF residual) of DPAM line data, dvcorr performs the diurnal correction of magnetic field data making use of Ground-Station magnetic field data, and ecomp and fcomp execute the compensation of aircraft's magnetic effect for compensation flight DPAM line data and for survey flight DPAM line data, respectively, using the method described by Nakatsuka and Okuma [2005b].

The programs ggrid and ggrids generate grid data from random points data (i.e., DPAM line data or equivalent) by using the method of "continuous curvature splines in tension" as developed by Smith and Wessel [1990]. Only the magnetic field (IGRF residual) data are converted into grid by the program ggrid, while ggrids applies the method also on the altitude data.

The program pframe draws an illustration of coordinates framework of survey area, for easier setting of grid location parameters, and pltrk plots trackline paths from DPAM line data. The programs pchkmag, pchkres and pchkcomp are used to plot magnetic field data in order to inspect the data quality. Total magnetic field data is plotted in black by pchkmag, and IGRF residual data is plotted in blue by pchkres, while pchkcomp plot IGRF residual magnetic data both before and after the compensation of aircraft's magnetic noise in red and blue colors.

The programs xslin and xslina extract simple standard line (StdLIN) data from random point data. Various formats of random point (line) data are supported by xslin as input data, while xslina deals with only DPAM line data but with the function of averaged re-sampling.

3. GDMP programs

GDMP group consists of programs for grid data manipulation: sel, seldb, altx, adjlv, gadd, gsub, govlay, gojoin, gmerge, txproj, gtrf, plmap, plmapc, plmapg, plmaps, plmapcs, shade, plmapl, plmapcl, xplmap and xplmapc.

The programs sel and seldb generate new grid data from an existing grid data with re-gridding. For sel program, the existing grid data should be assigned as the input filename parameter, while AMDB grid data (Nakatsuka et al., [2005]) is selected by the code name input for seldb. The standard format magnetic anomaly grid data includes both information of magnetic anomalies and the altitudes of observation surface. The program altx extracts the altitude information into a standard grid data format. The program adjlv adds a constant value to all grid data in order to adjust DC level, gadd adds 2 sets of grid data, and gsub subtracts the 2nd set grid data values from the 1st set grids, generating a new grid data.

The program govlay overlays multiple grid data sequentially with placing slit zone, gojoin joins multiple grid data sequentially with placing transient zone surrounding the overlaid data, and gmerge merges multiple grid data sequentially with placing transient zone along overlapping zone, to create a new grid data. The program txproj translates grid data into another map projection with re-gridding, and gtrf translates IGRF residual grid data into another reference model.

The programs plmap, plmapc, plmapg, plmaps, plmapcs and shade produce illustrations of 3D surface given as altitude distribution of grid data, by line contour map (plmap), color-graded contour map (plmapc), gray-scale grading map (plmapg), shaded-relief contour map (plmaps), shaded-relief color-grading map (plmapcs) and shaded-relief map (shade), all on an A4 sheet. The program plmapl generates line contour map with drawing trackline paths, and plmapcl generates color-graded contour map with trackline paths, both on an A4 sheet.

The programs xplmap and xplmapc create line contour map and color-graded contour map of grid data, respectively, with placing various supplementary items on a selectable size of sheet up to B0.

4. ANAM programs

ANAM group consists of programs for the reduction and quantitative interpretation of magnetic anomaly data: tmcorr, emag, emagf, amag, amagc, galtf, galts, lcecorr, aaptdp, emeq, ameq, ameqc, cmeq, emeqs, ameqs, ameqsc, cmeqs, rpmeqs, edeq, adeq, adeqc, cdeq, edeqs, adeqs, adeqsc, cdeqs and rpdeqs.

The program tmcorr executes the correction of topographic effect on the magnetic anomaly grid data, assuming a uniform terrain magnetization. The programs emag, emagf, amag and amagc are for the magnetization intensity mapping [Okuma et al., 1994; Nakatsuka, 1995]. First, emag or emagf calculates the contribution coefficients matrix COEF, and then amag or amagc executes inversion process of magnetization intensity mapping by the method of conjugate gradients (CG). The source model is approximated into rectangular blocks in emag, while the surface undulation is considered in emagf. The iteration of CG method is controlled by the number of loop count (amag) or by convergency tolerance (amagc).

The program lcecorr estimates the railway loop-current effect in the observed magnetic anomaly grid data, and correct for it. The program aaptdp is for the semi-automatic modeling by point-dipole sources for observed magnetic anomaly grid data. These processes were used in the data interpretation of AM data in the Kobe-Kyoto area [Nakatsuka et al., 2004].

The program galtf interpolates observation altitudes of StdLIN data into grid data (for the purpose of illustrating observation altitudes), and galts generates grid data of smoothed observation altitude from StdLIN data (for the purpose of using as the reference in altitude reduction).

The programs emeq, ameq, ameqc and cmeq are for the altitude reduction from grid data, and the programs edeq, adeq, adeqc and cdeq are for the altitude reduction from random point data, by equivalent anomaly method [Nakatsuka and Okuma, 2006a, 2006b]. First, emeq or edeq calculates the contribution coefficients matrix CMUP or CFUP, and ameq/ameqc or adeq/adeqc executes inversion process of equivalent anomaly derivation by CG method. The iteration of CG method is controlled by the number of loop count in ameq/adeq, or by convergency tolerance in ameqc/adeqc. Then the reduction calculation is performed by cmeq or cdeq program.

The programs emeqs, ameqs, ameqsc, cmeqs and rpmeqs are for the altitude reduction from grid data, and the programs edeqs, adeqs, adeqsc, cdeqs and rpdeqs are for the altitude reduction from random point data, by equivalent source magnetization method [Nakatsuka and Okuma, 2006a, 2006b]. First, emeqs or edeqs calculates the contribution coefficients matrix CMUPS or CFUPS, and ameqs/ameqsc or adeqs/adeqsc executes inversion process of deriving equivalent source magnetization by CG method. The iteration of CG method is controlled by the number of loop count in ameqs/adeqs, or by convergency tolerance in ameqsc/adeqsc. Then the reduction calculation is performed by cmeqs or cdeqs program, and the reduction-to-pole anomaly on that surface is calculated by rpmeqs or rpdeqs program.

5. Common features of the programs

I have developed these programs on a Linux machine, on which the GNU Compiler Collection [] is installed, and all programs are coded in FORTRAN language, except alog2asc program and a few subprograms (for the common use from FORTRAN programs) are coded in C language. The programs were developed step by step during the years of research activities on aeromagnetic survey at the Geological Survey of Japan, and they were verified to function well by some practical data. The program codes, however, have been revised further for the improvement on reasonable data handling or easier operation. So there is a little possibility that some unaware bugs are left. But, I convince the user will easily correct bugs, as the source codes are now presented.

In these programs, the subroutine library software by Nakatsuka [2005] is utilized quite frequently. So the installation of the subroutine library in advance is required to compile these programs.

Most programs use some array variables, and the necessary size (dimension) of the array is dependent on the actual data to be dealt with. (Although too large size definition does not harm the correct result, the memory resource may be wasted.) The definitions of such dimensions are often given by using 'parameter' statement in source codes, to enable easier adjustment. Actually the definitions of such parameters may not be consistent among programs even in the same group. The user is requested to customize them according to his/her requirement.

When executing these programs, a few to several or more parameters have to be specified, along with selecting working directory path. In these programs, the LWKDIR (Assist to set-up working dir. and process parameters, etc.) function with 'opnpin' mechanism of the subroutine library mentioned above is implemented. The use of this function provides us two merits. The one is that the parameter input is prompted with informative message or valid choice, which is very useful to learn how to operate the program. The other is that the same program is used to the non-interactive process for predetermined parameters with the help of 'opnpin' mechanism. The latter style of execution of program is realized by creating 'Job Control File' and submitting job/job1 command of the Utility program appended to the subroutine library. Refer to the section of Job Control File, Misc. Utility Program of GSJ Open-file Report no.422 [Nakatsuka, 2005], with respect to the detail of executing a series of job steps with the help of 'opnpin' mechanism.

6. Files included in this Report

Source programs coded in FORTRAN language have the filenames with .f (specifier) extensions, and those in C language with .c (specifier) extensions. To compile these programs under UNIX environment, @mkall script in each directory can be used, provided that the subroutine library by Nakatsuka [2005] is installed as the archive file libgm.a under the directory of /home/SHARE/lib. (Here, g77 and gcc are the commands to execute FORTRAN and C language compilers, as we use the GNU compilers.)

There are also HTML documentations for each of three program groups, DPAM, GDMP and ANAM. The explanation of the function of individual program is given briefly, and the parameters required to execute the program are described.

In the diskette appended to this report, HTML files, FORTRAN and C source files, and templates of Job Control files for non-interactive execution are included, in the following tree structure.

        <<< Tree structure of files included in this report >>>

 +-- openfile0449.html  Cover page HTML (in Japanese/English)
 +-- index.html      Introduction HTML in English (This document)
 +-- dpam.html       Documantation for DPAM programs
 +-- gdmp.html       Documantation for GDMP programs
 +-- anam.html       Documantation for ANAM programs
 +-- dpam/           (Directory containing DPAM source programs)
 |    |
 |    +-- @mkall     (Script to compile all DPAM programs below)
 |    +-- alog2asc.c, xldam.f, xldpn.f, xldhg.f, despike.f,
 |         dvcorr.f, sml5.c, ecomp.f, xldam.f, gsurf.c,
 |          ggrid.f, ggrids.f, pframe.f, pltrk.f,
 |           pchkdv.f, pchkmag.f, pchkres.f, pchkcomp.f,
 |            xslin.f, xslina.f
 +-- gdmp/           (Directory containing GDMP source programs)
 |    |
 |    +-- @mkall     (Script to compile all GDMP programs below)
 |    +-- sel.f, seldb.f, altx.f, adjlv.f, gadd.f, gsub.f,
 |         govlay.f, gojoin.f, gmerge.f, txproj.f, gtrf.f,
 |          plmap.f, plmapc.f, plmapg.f, plmaps.f,
 |           plmapcs.f, shade.f, plmapl.f, plmapcl.f,
 |            xplgobj.c, xplmap.f, xplmapc.f
 +-- anam/           (Directory containing ANAM source programs)
 |    |
 |    +-- @mkall     (Script to compile all ANAM programs below)
 |    +-- tmcorr.f, emagf.f, emag.f, amag.f, amagc.f,
 |         galtf.f, galts.f, lcecorr.f, aaptdp.f,
 |          emeq.f, ameq.f, ameqc.f, cmeq.f, emeqs.f,
 |           ameqs.f, ameqsc.f, cmeqs.f, rpmeqs.f,
 |            edeq.f, adeq.f, adeqc.f, cdeq.f, edeqs.f,
 |             adeqs.f, adeqsc.f, cdeqs.f, rpdeqs.f
 +-- dpam_tp/        (Directory containing templates for DPAM programs)
 |    |
 |    +--,,,,,
 | ,,,
 +-- gdmp_tp/        (Directory containing templates for GDMP programs)
 |    |
 |    +--,,,,,,
 | ,,,,
 |  ,
 +-- anam_tp/        (Directory containing templates for ANAM programs)


  1. Nakatsuka, T. [1995] Minimum norm inversion of magnetic anomalies with application to aeromagnetic data in the Tanna area, Central Japan. J. Geomag. Geoelectr., 47, 295-311.
  2. Nakatsuka, T. [2005] Calculation of the International Geomagnetic Reference Field (4). GSJ Open-file Report, no.423, 39p.+ 1 Diskette, Geol. Surv. Japan, AIST. (in Japanese with English abstract)
  3. Nakatsuka, T. [2006] Library software for geophysical data processing and representation (2). GSJ Open-File Report, no. 442, 102p.+ 1 CD-ROM, Geol. Surv. Japan, AIST.
  4. Nakatsuka, T., and S. Okuma [2005a] Data processing system for helicopter-borne aeromagnetic survey. Proc. 113th SEGJ Conference, 239-242. (in Japanese with English abstract)
  5. Nakatsuka, T., and S. Okuma [2005b] Development of a helicopter-borne magnetic survey system in use of a nose boom magnetic sensor and its passive compensation for aircraft's magnetic field. Butsuri-Tansa (Geophys. Explor.), 58, 451-459. (in Japanese with English abstract)
  6. Nakatsuka, T., and S. Okuma [2006a] Reduction of geomagnetic anomalies from the observation at varying elevations with helicopter survey. Intl. Symp. Airborne Geophysics 2006, Tsukuba, January, 2006.
  7. Nakatsuka, T., and S. Okuma [2006b] Reduction of geomagnetic anomaly observations from helicopter surveys at varying elevations. Explor. Geophys., 37, 121-128; Butsuri-Tansa (Geophys. Explor.), 59, 121-128; Mulli-Tamsa (Geophys. Explor.), 9, 121-128.
  8. Nakatsuka, T., S. Okuma, R. Morijiri, and M. Makino [2004] High-resolution aeromagnetic anomaly map of the Kobe Kyoto area (1:100,000). Total Intensity Aeromagnetic Maps, no.42, Geol. Surv. Japan, AIST.
  9. Nakatsuka, T., S. Okuma, M. Makino, and R. Morijiri [2005] Aeromagnetic Survey Database of Japan. Digital Geoscience Maps, P-6 (2 CD-ROM disks), Geol. Surv. Japan, AIST.
  10. Okuma, S., M. Makino, and T. Nakatsuka [1994] Magnetization intensity mapping in and around Izu-Oshima Volcano, Japan. J. Geomag. Geoelectr., 46, 541-556.
  11. Smith, W.H.F., and P. Wessel [1990] Gridding with continuous curvature splines in tension. Geophysics, 55, 293-305.