- Email: [email protected]
A fortran IV program for interpolating irregularly spaced data using the difference equations for minimum curvature
Comp,t,ra & Geoscumcts, rot. 4. p. 209 Pcrlpumon Press Ltd.. 197t. Printed in Great Britain
009g-~04Dlk10~l-0'209r~2.00/0
LETTERS TO THE EDITOR CLARIFICATION OF DAVID'S METHOD FOR ORE ESTIMATION
to deal with data containing large gaps. A rule of thumb for finding a suitable value of L (which must be a power of 2) is:
In Clark (1977), three references are made to David (1976). The problem of estimating a mining block on a bench from the surrounding DDH is discussed. It is suggested that "'David's method" ignores all samples above and below the bench. If this were the situation, we could not agree more with Clark, this would be an oversimplification which certainly would loose information. It should be pointed out however that this is absolutely not the situation and samples above and below the bench are certainly included in all our programs. More details are presented in David (1977). What may not be legitimate is the use of the expression "2-D". A block is collapsed as a square, a sample, the same length as the bench height is collapsed as a point but there is no need for both to be in the same 2D plane! Clark (1977) is an original contribution, why not present it for itself rather than "against" something else. When will geostatisticians stop beating drums and raising flags?
L = (largest gap size) + (grid increment x F),
Clark, 1.. 1977. Practical kriging in three dimensions: Computers & Geosciences. v. 3. no. I. p. 173-180. l~vid. M., 1976,The practice of kfigin8in: Advanced geostatistics in the mining industry: D. Reidel. Dordrecht-Holland,p. 31-.48. David, M.. 1977. Geostatistical ore reserve estimation: Elsevier Scientific Publ. Co.. Amsterdam. 364 p.
Ecole Polytechnique, Montreal, Quebec, Canada.
MICHELDAVID
ERRATUM--A FORTRAN IV PROGRAM FOR INTERPOLATING IRREGULARLY SPACED DATA USING THE DIFFERENCE EQUATIONS FOR MINIMUM CURVATURE On reading Waldo Tobler's letter (Computers & Geosciences v. 3, no. I, p. 181) it struck me how unlucky he was to choose the combination of L = 4, N (or M) = 13 to test the program. This is not to excuse myself for failing to foresee the possibility of it happening, and the first modification he suggests is valid. The second is rather a matter of opinion: the method used in the original is merely to ensure that the program runs even if incorrect parameters are input (at the expense of possibly cutting down the map area slightly). The program was intended for large arrays of grid points. The purpose of introducing the L parameter was
where F is between 5 and 10. So, if N is 13, L would not need to be larger than 2. What actually happens in the boundary zone (outer two rows and columns) is that the "off-grid" equations are not applied but instead the data point is moved to the nearest grid point and this value is fixed and controls the surrounding grid values during iteration. This leads to slight inaccuracies but is not important if a large array is used. To simplify the program DUIJ is extracted from an area which excludes the boundary zone. It should be noted that the spectral content of the data may be modified if the grid spacing is more than about half the minimum data spacing, but this probably will not be serious unless it exceeds the minimum data spacing in which case, of course, some data points will be omitted. A slightly more subtle error in the program has come to light recently. In subroutine SORT the line: IF (R.EQ.M) GO TO 5 should occur after line 25 (MARK(I J) = R). The reason is that FORTRAN compilers (unlike ALGOL compilers) test the DO loop index at the end of the loop so that at • least one execution of the loop always occurs. The effect of this error is that if no data occur in one of the sorting rectangles the data points may be sorted incorrectly and in consequence the message SEARCH RADIUS TOO SMALL may be printed when it does not seem to be justified. A larger search radius would cure the problem--but at the expense of increased execution time. The program now has been modified so that the data are packed into the reference array IU. This has the effect that, if subroutine HEIGHT is modified suitably, then subroutine SORT becomes redundant (the data being effectively sorted in IU) along with all the arrays except U, IU, and MAP, which makes the program easier to use. The contour following routine mentioned in the I~aper was implemented but using a weighted paraboloid algoritbm to interpolate the value within the grid (use of the "off-grid" equations leads to discontinuities near data points). It is possible to include the data points in this interpolation so that contours precisely honor them. The results are satisfactory and could be improved further if the interpolation could somehow be based on the minimum curvature criterion.
209
Department of Geolosy, Leicato Univereity, Leicester,England
C. J. SWATH
009g-~04Dlk10~l-0'209r~2.00/0
LETTERS TO THE EDITOR CLARIFICATION OF DAVID'S METHOD FOR ORE ESTIMATION
to deal with data containing large gaps. A rule of thumb for finding a suitable value of L (which must be a power of 2) is:
In Clark (1977), three references are made to David (1976). The problem of estimating a mining block on a bench from the surrounding DDH is discussed. It is suggested that "'David's method" ignores all samples above and below the bench. If this were the situation, we could not agree more with Clark, this would be an oversimplification which certainly would loose information. It should be pointed out however that this is absolutely not the situation and samples above and below the bench are certainly included in all our programs. More details are presented in David (1977). What may not be legitimate is the use of the expression "2-D". A block is collapsed as a square, a sample, the same length as the bench height is collapsed as a point but there is no need for both to be in the same 2D plane! Clark (1977) is an original contribution, why not present it for itself rather than "against" something else. When will geostatisticians stop beating drums and raising flags?
L = (largest gap size) + (grid increment x F),
Clark, 1.. 1977. Practical kriging in three dimensions: Computers & Geosciences. v. 3. no. I. p. 173-180. l~vid. M., 1976,The practice of kfigin8in: Advanced geostatistics in the mining industry: D. Reidel. Dordrecht-Holland,p. 31-.48. David, M.. 1977. Geostatistical ore reserve estimation: Elsevier Scientific Publ. Co.. Amsterdam. 364 p.
Ecole Polytechnique, Montreal, Quebec, Canada.
MICHELDAVID
ERRATUM--A FORTRAN IV PROGRAM FOR INTERPOLATING IRREGULARLY SPACED DATA USING THE DIFFERENCE EQUATIONS FOR MINIMUM CURVATURE On reading Waldo Tobler's letter (Computers & Geosciences v. 3, no. I, p. 181) it struck me how unlucky he was to choose the combination of L = 4, N (or M) = 13 to test the program. This is not to excuse myself for failing to foresee the possibility of it happening, and the first modification he suggests is valid. The second is rather a matter of opinion: the method used in the original is merely to ensure that the program runs even if incorrect parameters are input (at the expense of possibly cutting down the map area slightly). The program was intended for large arrays of grid points. The purpose of introducing the L parameter was
where F is between 5 and 10. So, if N is 13, L would not need to be larger than 2. What actually happens in the boundary zone (outer two rows and columns) is that the "off-grid" equations are not applied but instead the data point is moved to the nearest grid point and this value is fixed and controls the surrounding grid values during iteration. This leads to slight inaccuracies but is not important if a large array is used. To simplify the program DUIJ is extracted from an area which excludes the boundary zone. It should be noted that the spectral content of the data may be modified if the grid spacing is more than about half the minimum data spacing, but this probably will not be serious unless it exceeds the minimum data spacing in which case, of course, some data points will be omitted. A slightly more subtle error in the program has come to light recently. In subroutine SORT the line: IF (R.EQ.M) GO TO 5 should occur after line 25 (MARK(I J) = R). The reason is that FORTRAN compilers (unlike ALGOL compilers) test the DO loop index at the end of the loop so that at • least one execution of the loop always occurs. The effect of this error is that if no data occur in one of the sorting rectangles the data points may be sorted incorrectly and in consequence the message SEARCH RADIUS TOO SMALL may be printed when it does not seem to be justified. A larger search radius would cure the problem--but at the expense of increased execution time. The program now has been modified so that the data are packed into the reference array IU. This has the effect that, if subroutine HEIGHT is modified suitably, then subroutine SORT becomes redundant (the data being effectively sorted in IU) along with all the arrays except U, IU, and MAP, which makes the program easier to use. The contour following routine mentioned in the I~aper was implemented but using a weighted paraboloid algoritbm to interpolate the value within the grid (use of the "off-grid" equations leads to discontinuities near data points). It is possible to include the data points in this interpolation so that contours precisely honor them. The results are satisfactory and could be improved further if the interpolation could somehow be based on the minimum curvature criterion.
209
Department of Geolosy, Leicato Univereity, Leicester,England
C. J. SWATH