Computers
PII: SOO!l8-3004(%)00077-5
& Geoscienws Vol. 23, No. I, pp. 127-131, 1997 k? 1997 Elsevier Science Ltd Printed in GreatBritain. All rights reserved 0098-3004/97 $17.00 + 0.00
SHORT NOTE PROFILE: A MICROSOFT QUICK BASIC PROGRAM FOR RETRIEVING DATA ALONG A GIVEN PROFILE FROM GRIDDED DATA Y. SATYANARAYANA,
D. CHANCHAL,
and G. R. K. MURTY
Government of India, Ministry of Defence, Defence Research and Development Org., Naval Physical and Oceanographic Laboratory, BMC Post, Thrikkakara, Kochi021, India (e-mail:
[email protected]) (Received
19 December
1995; revised 12 August 1996)
INTRODUCTION A popular operation performed on any geoscientific contour map (such as gravity, magnetic, bathymetry or oceanographic data) is the construction of twodimensional (2-D) profiles between two given points for further qualitative or quantitative analysis. Generally the process of interpolation is done manually by digitising the contour values to retrieve the data along the profile. This process is time consuming and tedious. Programs for this type of operation are commercially available as a part of integrated software packages whose source code is not freely available. In this note, a simple program PROFILE written in Microsoft Quick Basic is presented for constructing 2-D profiles. The program utilises a 2-D interpolation technique to retrieve the data without deforming the original trend along any profile from a gridded data set. The significant features of the program are:
either in ASCII or binary the commercial package Software Inc. (1994).
format generated from SURFER of Golden
INTERPOLATION
TECHNIQUE
The symmetrical weighting interpolation formula used in this program is derived as follows: zi, z2, z3, and z4 are the z values at four corners of the grid cell, see Figure 1. The values of z at points ziz and z34 are given by 212 = (ZIXe + zzxW)/& and 234 = (z3-G
+ Z4XwY&v
where g, = (x, + x,) is the grid spacing in x direction. Therefore the value of Z at point z is given by Z =
(ZlZYs
+
ZUY”)/&’
=
be(w,
+
z3yd+
(1) The profile can be defined in two ways: by Xw(Z2Ys + z4Y”MMy~
two end points of the profile, and by the starting point and the bearing of the profile. interpolation technique is (2) An appropriate used to represent the original trend of the data with little deformation. (3) It accepts Surfer grid files both in ASCII and binary formats. (4) It works with files of any size since arrays are not used. (5) The program is user friendly.
where g,, = (_yn+ y,) is the grid spacing in the y direction. The values obtained from this formula retain the original trend of the data with little deformation.
DESCRIPTION
OF THE PROGRAM
The program is written in Microsoft Quick Basic (Microsoft Corporation, 1989a,b). All statements are standard. This program can also be implemented in QBASIC supplied with the DOS 5.0 (or later) version of the operating system. Upon execution, the user will be asked to enter the name of input data file (e.g Surfer grid file) in either ASCII or binary format. The type of format will be ident-
The program takes any regularly spaced data either in square or in rectangular form. This program is helpful in constructing 2-D profiles from regularly spaced data such as those collected in ground geophysical surveys for mineral or geotechnical investigations. It also accepts directly any gridded data file 127
Short Note
128
-I-
+ Yn I t
I
----_-_ xw
‘2 -----------_xe
7
I
Ys
z4
234
1 =3
Figure 1. Plan view of one grid cell: Z, interpolated value at point z; z,, z2, z3 and z4 are z values at four nearest surrounding stations; xW +x, = g,; y, + y, = g!.
ified by a header label in the grid file by the program. A data set SAMPLE.DAT in Surfer grid ASCII format used to illustrate the program PROFILE is given in Appendix A. The function ReadGridValue() retrieves the Z(ij) value from the binary grid file. The function computes the byte position of element (ij), uses the SEEK statement to go there and the GET statement to read it. This is fast and allows the program to work with any size file because arrays are not used. In the situation of an ASCII file, the grid is converted into a temporary binary grid file in order to avoid reading the entire grid into memory. To construct a profile, the user selects one of the two options available for defining the limits of the profile. Under option 1 the user must supply the coordinates of the two end points (A and B in Fig. 2) of the profile, while under option 2, a starting point and the bearing of the profile with respect
10.4
10.5
Figure 2. Input options for constructing 2-D profile from gridded data: ., data point; AB, option 1; CD, option 2 (0 is the bearing of the profile); g,, g,., grid spacings along x and y directions.
to the positive y-axis of the map (the point C and angle 0 in Fig. 2) are required. The boundary values (the limits) of the four corners of the map will be displayed to enable the user to give the correct input. If the user enters a point outside the boundary, the program cautions the user with an error message. Under the first option, the data will be retrieved at a specified sampling interval along the profile by interpolation using the function INTERPOL, without any change in trend of the values. The second option activates the subroutine INTERSEC before it calls the INTERPOL subroutine. The subroutine INTERSEC computes the end point of the profile at the map boundary. It is suggested that selected the sampling interval should always be less than the station spacing in the grid in order to obtain a smooth curve. In this program a maximum of 5 samples in a grid cell (i.e. ndrs = 5) is used. The user can change the value of ndrs.
10.6
10.7
10.8
10.9
Figure 3. Sample contour map using data given in Appendix A; AB is section along which profile is required.
I29
Short Note
0.5
0 Distance
Figure 4. Graphic
display
(units)
of profile AB of Figure 3.
The subroutine PLOT enables instant graphic display of the profile. A sample contour map using the gridded data SAMPLE.DAT in the Surfer ASCII format (Appendix A) is shown in Figure 3. The values of to illustrate the prox, y, z are taken arbitrarily gram. The input data are uniformly spaced in a projected normal reference frame (i.e. the axes are orthogonal). The graphic display of the profile along AB on the contour map (Fig. 3) is shown in Figure 4. The distances from the starting point along the profile and z values are plotted along the x-axis and y-axis, respectively. The diagonal length of the map is selected as the maximum value of x. The values of Z,in and z,,, are taken as y-axis limits in the current display. The user can change the axis limits by changing the values of the window limits in the program.
After each run, the program asks for a new profile. The user can save the profile data in a file by selecting the output option. The sequence of display on execution of the program is given below for option 1 in Appendix B. The program PROFILE is available by anonymous FTP from the server IAMG.ORG. The same program can be used, with minor modifications, for retrieving data and constructing stacked 2-D profiles of different parameters for the area such as bathymetry, gravity and magnetic simultaneously. Acknowledgments-The authors would like to thank the Director of NPOL for approving this work for publication and John A. Grant,* -Airbo& Geophysics Section, Geological Survey of Canada, for critical review and valuable suggestions in improving this manuscript.
REFERENCES Microsoft Corporation, 1989a, Microsoft BASIC, BASIC Language Reference, Ver. 7.0: Microsoft Corporation, Redmond, Washington, variously paged. BASIC, 1989b, Microsoft Microsoft Corporation, Programmer’s Guide, Ver. 7.0: Microsoft Corporation, _ Redmond, Washington, variously paged. Golden Software Inc., 1994, SURFER for Windows Guide: Golden Software Inc., Golden, Users’ Colorado, variously paged.
Short Note
130
APPENDIX
A
Test Dataset in Surfer ASCII Format D6M 15 17 10.4
1O.B
2 2.4 10
500
14.2168
40.1126
6a.BS64
S1.2WlS
S76.6a6
Ha.721
44B.704
417.4
46.21@5
54.4976
70.262S
S6.6la7
47.2172
46.a66
a62.206
414.062
455.766
456.624
464.655
76.6013
61.2623
65.6477
71.2745
76.6261
4la.612
4aS.lla
449.01a
461.644
46a.662
66.6164
lOa.
390.641
414.507
140.422
0
a78.655
a66.81@
116.63 411.169
lob.276
67.0925
440.541
124.017
108.665
104.22s
177.66
107.508
@@.05@6 lla.624
177.702
155.626
a40.264
a47.56S
101.004 a06.229
445.16a
aSO.a4@ la6.946
162.al6
256.658
2a6.4a7
217.161
21O.BOa
160.057
aa1.646
266.601
aaO.603
a3a.BO4
2S6.261
266.646
265.117
274.246
215.BS6
aa6.696
b76.6711
aS4.741
346.657
SS1.001
261.646
260.666
264.+6
a61.aO2
361.5
274.661
2711.266 a66.674
266.18@
a06.626
145.224
asa.
266.87 331.16a aO6.016
all.67 826.16a ala.B22
265.76
824.09
a25.66
S24.866
[email protected]
a41.715 aos.166 a28.626
2BS.46B aa7.724 aO7.a51
aO6.64)
a25.@71
aa)b.aSO
a2a.666
alS.161
816.137
aPa.
306.aa6
S14.12
a41.662
alS.076
a26.222 a14.652
al6.196
a66.a26 all.046
alS.a7a
a66.616
401.766
a44.166
aS7.65
a17.62
a46.441
a62.6a7
267.666
261.767
a42.0114
a6a.S6a
2a6.6a6
275.666
aa6.504
140.14
16s.066
206.462
246.542
267.a27
176.a611
160.411
206.664
254.1s
167.666
al6.1al
ala.6
211.441
26S.22a
26a.S6S
ala.@47
aaO.211
al2.611
al4.OBl
406.656
420.74s
[email protected]
2Sa.Sal
a26.6a6
267.2a7
266.146
26S.a66
271.626
267.267
262.la6
276.64a
26a.47
272.a64
266.617
264.127
a40.066
PS7.1a2
a25.486
307.66
26a.061
207.01)s
a56.464
8a7.#4
a41.262
a2a.766
a60.66
a21.2S6
ao5.146 a40.674
aO4.045
a26.027
676.061
all.665
26B.lOa 327.15
aos.84 a20.12
226.W7
a44.176
aa1.61
281.661
aa4.aO6
Sas.406
a60.601
a46.7)
a60.7#1
446.6s
a80.022
110.741
a46.716
26a.746
274.216
146.604
180.966
126.144
272.415
419.253
a44.611
aaS.
432.12a
a62.111
164.047
262.276
lsO.172
361.066
a66.21
a44.762
4aa.486
lB1.601
SO0 42a.626
26a.764
121.7aP
114.662
461.221
104.161
64.2296
124.174
491.a66
a16.266
272.166
812.27a
275.876
all.a66
276.65
271.72
271.446
265.417
aal.
265.07
827.226
264.6a2
916.622
alS.106 a2#.6M a18.064
271.a66
131
Short Note APPENDIX
B
Example Dialogue
PROBILE.BAs
>
1 rlbo End Pointst Option 2 tstarting Point 8 Bearing
Option Enter
c A:\SAMPLE.DAT
I 3
Enter Surfer grid file Name
'1' or '2' 3
After entering
Cl>
1 or 2, the minimum & maximum values of x and
y will be displayed
on the screen.
If the option 1 is selected then display will be I ?
< 10.4,2.1
>
(x,y) : ?
c 10.8,2.2
>
Enter Starting Point Enter End Point
The output plot of the profile will be displayed the comments
(Fig.4):
Enter cT>ry again
: cS>ave
t
uit
along with