GSJ Open-File Report, no. 519 December, 2009

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

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.]. Their source codes are open to public through the GSJ Open-file Reports.

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.

This report is a revised version for the previous GSJ Open-File Report no.449 [Nakatsuka, 2007]. The points of revision are:

  1. Programs gtrim and gtopo were added to DPAM group,
  2. Programs tmcfix, cmag, cmagf, plamag, plamagc and calmas were added to ANAM group,
  3. All programs are updated so as to record the time-stamp of execution on the log file,
  4. Response to the change of specification of the library subprogram 'libgm', and improved effective use of the library,
  5. Modification of coding for minor improvements,
  6. Japanese documentation was added.
Be aware that the revised version of 'libgm' (GSJ Open-file Report no.518 [Nakatsuka, 2009]) shoule be combined with to compile the source codes in this report.

2. DPAM programs

DPAM group consists of programs for aeromagnetic survey data processing as follows.
alog2asc, xldam, xldpn, xldhg, despike, dvcorr, ecomp, fcomp, ggrid, ggrids, pframe, pltrk, pchkdv, pchkmag, pchkres, pchkcomp, xslin and xslina.

The program alog2asc converts aeromagnetic binary raw data (specialized for our survey system) into ASCII file format, the programs xldam and xldpn generate line data (DPAM line data) file from aeromagnetic ASCII raw data, making use of Real-time GPS and PostNav processed DGPS position data, respectively, and the program xldhg generates similar line data (HGAM line data) file from helicopter-borne gradiometer aeromagnetic data acquired. 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, the program dvcorr performs the diurnal correction of magnetic field data making use of Ground-Station magnetic field data, and the programs 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 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 the program 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 plots 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, gtrim, govlay, gojoin, gmerge, txproj, gtopo, 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, the program gadd adds 2 sets of grid data, and the program gsub subtracts the 2nd set grid data values from the 1st set grids, as generating a new grid data. The program gtrim replaces values of input grid data with NULL (undefined) data value for the range where the values of reference grid data are NULL data value, as generating a new grid data.

The program govlay overlays multiple grid data sequentially with placing slit zone, the program gojoin joins multiple grid data sequentially with placing transient zone surrounding the overlaid data, and the program 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, the program gtopo generates topographic altitudes grid data from DEM40** (Digital Elevation Model 40m) data, and the program 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 the program 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.

* AMDB grid files should be recovered from TGZ file and stored in the directory /pub/AMDB/DATA/grd with filenames as "code-name.gd".
** DEM40 is a group of files of terrain elevation data in Japan, extracted from Digital Maps (50m mesh elevation, Japan-I, Japan-II and Japan-III) by the Geographical Survey Institute, Japan (GSI). The source data are subject to GSI's copyright. The procedure how to generate DEM40 files from GSI's Digital Maps is presented in genDEM40j.html (only in Japanese).

4. ANAM programs

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

The programs emag, emagf, amag, amagc, cmag and cmagf are used for the magnetization intensity mapping [Okuma et al., 1994; Nakatsuka, 1995]. First, emag / emagf calculates the contribution coefficients matrix COEF, then amag / amagc executes inversion process of magnetization intensity mapping by the method of conjugate gradients (CG), and cmag / cmagf may be used to calculate the synthetic magnetic anomalies on a desired altitude surface by using the result of amag / amagc. The source model is approximated into rectangular blocks in emag / cmag, while the surface undulation is considered in emagf / cmagf. The iteration of CG method is controlled by the number of loop count (amag) or by convergency tolerance (amagc). The programs plamag and plamagc plot the result of the magnetization intensity mapping in the form of line contour map and color-graded contour map, respectively, with masking the surrounding source region.

The program tmcorr executes the correction of topographic effect on the magnetic anomaly grid data, assuming a uniform terrain magnetization, while the program tmcfix executes the correction of topographic effect with the assigned value of terrain magnetization. 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 aeromagnetic 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 the program 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.

The program calmas calculates theoretical magnetic anomaly distribution on the specified observation surface, caused by a simple source model.

5. Common features of the programs

I have developed these programs on a Linux machine, on which the GNU Compiler Collection [http://www.gnu.org/software/gcc/] 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 'libgm' by Nakatsuka [2009] is utilized quite frequently. So the installation of the subroutine library 'libgm' 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 'libgm' 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.518 [Nakatsuka, 2009], 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 'libgm' by Nakatsuka [2005] is installed as the archive file libgm.a under the directory of /home/SHARE/lib, and the command names 'fort' and 'gcc' is aliased to (or the name of) the desired FORTRAN and C 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 CD-ROM 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 >>>

(root)
 |
 +-- 0519index.html  Cover page HTML
 +-- index.html      Introduction HTML in English (This document)
 +-- dpam.html       English Documantation for DPAM programs
 +-- gdmp.html       English Documantation for GDMP programs
 +-- anam.html       English Documantation for ANAM programs
 |
 +-- indexj.html     Introduction HTML in Japanese
 +-- dpamj.html      Japanese Documantation for DPAM programs
 +-- gdmpj.html      Japanese Documantation for GDMP programs
 +-- anamj.html      Japanese Documantation for ANAM programs
 +-- genDEM40j.html  (in Japanese) How to generate DEM40 file
 |
 +-- dpam.tgz        Gzipped tar archive of the directory 'dpam'
 +-- gdmp.tgz        Gzipped tar archive of the directory 'gdmp'
 +-- anam.tgz        Gzipped tar archive of the directory 'anam'
 +-- dpam_tp.tgz     Gzipped tar archive of the directory 'dpam_tp'
 +-- gdmp_tp.tgz     Gzipped tar archive of the directory 'gdmp_tp'
 +-- anam_tp.tgz     Gzipped tar archive of the directory 'anam_tp'
 |
 +-- 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, fcomp.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,
 |         gtrim.f, govlay.f, gojoin.f, gmerge.f, txproj.f,
 |          gtopo.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)
 |    +-- emag.f, emagf.f, amag.f, amagc.f, cmag.f, cmagf.f,
 |         plamag.f, plamagc.f, tmcorr.f, tmcfix.f, lcecorr.f,
 |          aaptdp.f, galtf.f, galts.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 calmas.f
 |
 +-- dpam_tp/        (Directory containing templates for DPAM programs)
 |    |
 |    +-- alog2asc.tp, xldam.tp, xldpn.tp, xldhg.tp, despike.tp,
 |         dvcorr.tp, ecomp.tp, fcomp.tp, ggrid.tp, ggrids.tp,
 |          pframe.tp, pltrk.tp, pchkdv.tp, pchkmag.tp,
 |           pchkres.tp, pchkcomp.tp, xslin.tp, xslina.tp
 |
 +-- gdmp_tp/        (Directory containing templates for GDMP programs)
 |    |
 |    +-- sel.tp, seldb.tp, altx.tp, adjlv.tp, gadd.tp, gsub.tp,
 |         gtrim.tp, govlay.tp, gojoin.tp, gmerge.tp, txproj.tp,
 |          gtopo.tp, gtrf.tp, plmap.tp, plmapc.tp, plmapg.tp,
 |           plmaps.tp, plmapcs.tp, shade.tp, plmapl.tp,
 |            plmapcl.tp, xplmap.tp, xplmapc.tp
 |
 +-- anam_tp/        (Directory containing templates for ANAM programs)
      |
      +-- emag.tp, emagf.tp, amag.tp, amagc.tp, cmag.tp, cmagf.tp,
           plamag.tp, plamagc.tp, tmcorr.tp, tmcfix.tp, lcecorr.tp,
            aaptdp.tp, galtf.tp, galts.tp, emeq.tp, ameq.tp,
             ameqc.tp, cmeq.tp, emeqs.tp, ameqs.tp, ameqsc.tp,
              cmeqs.tp, rpmeqs.tp, edeq.tp, adeq.tp, adeqc.tp,
               cdeq.tp, edeqs.tp, adeqs.tp, adeqsc.tp, cdeqs.tp,
                rpdeqs.tp, calmas.tp

References

  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. [2007] Software system for Aeromagnetic data processing, Grid data manipulation, and Reduction and quantitative interpretation of magnetic anomaly data. GSJ Open-file Report, no.449, 29p. + 1 Diskette, Geol. Surv. Japan, AIST.
  4. Nakatsuka, T. [2009] Library software for geophysical data processing and representation (3). GSJ Open-File Report, no.518, 107p.+ 1 CD-ROM, Geol. Surv. Japan, AIST.
  5. 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)
  6. 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)
  7. 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.
  8. 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.
  9. 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.
  10. 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.
  11. Okuma, S., M. Makino, and T. Nakatsuka [1994] Magnetization intensity mapping in and around Izu-Oshima Volcano, Japan. J. Geomag. Geoelectr., 46, 541-556.
  12. Smith, W.H.F., and P. Wessel [1990] Gridding with continuous curvature splines in tension. Geophysics, 55, 293-305.