- Email: [email protected]
An alternative method for calculating variogram surfaces using polar coordinates
ARTICLE IN PRESS
Computers & Geosciences 31 (2005) 801–803 www.elsevier.com/locate/cageo
Note
An alternative method for calculating variogram surfaces using polar coordinates$ F. Iwashita, R.C. Monteiro, P.M.B. Landim Geomathematic Laboratory, Department of Applied Geology, Sa˜o Paulo State University/Rio Claro campus, Av 24-A, 1515, Rio Claro 13506-900, Brazil Received 12 August 2004; received in revised form 15 December 2004; accepted 7 January 2005
1. Introduction
2. Variogram surface
Geostatistics deals with spatial and temporal data from natural phenomena and involves estimating values of regionalized variables. In this analysis the fundamental tool is the variogram that is a variability measure; for example, geological variability related to distance. The variogram shows the degree of spatial dependence between samples for a specific support. The squared differences between pairs of values are estimated assuming stationarity in increments. The semivariogram is calculated as
A variogram is a measure of spatial variability with respect to distance, but such variability can be different when taking into account different directions. This occurs with anisotropic phenomena. In this situation it is important to calculate variograms in different directions. This procedure may be carried out by constructing a variogram surface that contemplates all the possible directions of variability. Examples of variogram surfaces are shown in Pannatier (1996) and Deutsch and Journel (1998). This note presents an alternative approach to obtaining a variogram surface using polar coordinates. To implement this routine it is necessary to convert gðhÞ from Cartesian to polar coordinates to simplify programming: ( x ¼ h cos ðyÞ; y ¼ h sen ðyÞ;
gðhÞ ¼
1 XN h ½Zðxi þ hÞ Zðxi Þ2 , i¼1 2N h
where h ¼ basic lag, Nh ¼ number of pairs. Variogram analysis deals with the experimental variogram calculated from the data and the variogram model fitted to the data. The variogram model is chosen from a set of mathematical functions that describe the spatial relationship. This step involves the knowledge and researcher’s experience concerning the phenomenon and the study area. A good adjustment is necessary for achieving a good result.
$ Code on server at http://www.iamg.org/CGEditor/ index.htm. Corresponding author. Tel.: +55 19 3526 2811; fax: +55 19 3524 2445. E-mail address: [email protected] (F. Iwashita).
where x and y are Cartesian coordinates, h (lag) is the radius and the angle y is the interval (angular tolerance) to each direction between 901 and 901.
3. Program description The routine polarmap.m, running in MATLAB (Mathworks, 2000), is part of a software package for exploratory and geostatistical analysis developed at the
0098-3004/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2005.01.009
ARTICLE IN PRESS 802
F. Iwashita et al. / Computers & Geosciences 31 (2005) 801–803
Geomathematics Laboratory in the Sa˜o Paulo State University/Rio Claro campus. The whole package has other routines for calculating histograms, scatter plots, box-whisker plots, normal probability plots and linear regression. The calculations for the experimental variogram and ordinary kriging are made by the functions vario.m and kriging.m from the Bayesian Maximum Entropy Library/BMELib package (Christakos et al., 2002). These routines are the basis for calculating and plotting indicator variograms and indicator kriging. The package allows calculation and display of standard deviation maps for ordinary and indicator kriging and a specific routine for crossvalidation. An index for experimental variogram fitness to the theoretical variogram is proposed, based on the w2 (chi-square) test: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN ðgðhÞ g ðhÞÞ2 g ðhÞ j¼1 , i¼ N where i is the calculated adjusted index; smaller the index, better the fitting; N is the number of elements in g domain; gðhÞ are the experimental variogram values; g ðhÞ are the theoretical model values for the fitted variogram. The BMELIB routine vario.m produces directional variograms (901 a 901) with angular tolerance defined by the user. The polarmap() routine calculates n (angular tolerance /180) directional variograms. The gðhÞ values calculated by each directional variogram are stored in a matrix (escalaC) in vector form. A polar mesh is created based on angular tolerance defined by the user, and the routine pcolor( ) plots the matrix values in escalaC as an index for a color scale map.
Table 1 Extract from example input file in GeoEAS format
3.1. Input data
Fig. 1. Variogram map by GSLIB (Deutsch and Journel, 1998).
Example. dat—Geostatistical Assessment Software 5 Easting Northing Arsenic Cadmium Lead 288.0 311.0 .850 285.6 288.0 .630 273.6 269.0 1.02 280.8 249.0 1.02 273.6 231.0 1.01 ^ ^ ^ 465.6 216.0 .930 492.0 216.0 .750 345.6 216.0 1.45
Environmental
11.5 8.50 7.00 10.7 11.2 ^ 11.6 6.90 9.90
18.25 30.25 20.00 19.25 151.5 ^ 25.00 33.00 40.75
Cadmium: Variogram map 162.500
Coordinates N-S
30.000 28.000 26.000 24.000 22.000 20.000 18.000 16.000 14.000 12.000 10.000 8.000 6.000 4.000 2.000 0.0
-162.500 -162.500
Coordinates E-W
162.500
The function readGeoeas.m (BMELib package) reads and transforms the database to the vector and matrix structures used by the analysis. The datafile is in GeoEAS format (ASCII) to be read, preferentially using a ‘.dat’ extension. The chosen example comes from Englund and Sparks (1988), using the variable cadmium (Table 1). 3.2. Example To compare the conversion by the polarmap.m routine, three variogram surfaces are presented. The first one was calculated by GSLib (Fig. 1), the second by Variowin (Fig. 2) and the third obtained by polarmap.m routine (Fig. 3). The datasets are the same, but there are some differences in the variogram maps due to inherent differences between Cartesian and polar coordinates and slightly different scales.
Fig. 2. Variogram map by Variowin (Pannatier, 1996).
ARTICLE IN PRESS F. Iwashita et al. / Computers & Geosciences 31 (2005) 801–803
803
involved, but the routine developed here for polar coordinates in MATLAB has just 49 lines. Its simplicity is related to the use of MATLAB toolboxes and their simple and versatile commands for matrix operations. The programming is straightforward, because directional variograms are calculated based on angles and not on Cartesian coordinates, and are more easily implemented than the calculations using Cartesian coordinates. The resulting variogram surfaces are effective representations of the spatial variability and anisotropic properties. This routine and the whole package can be downloaded from the IAMG server or obtained from the author at E-mail to: [email protected].
Acknowledgements To ‘‘The State of Sa˜o Paulo Research Foundation/ FAPESP’’ for a scholarship to the first author of this paper.
References
Fig. 3. (a) Variogram map by polarmap.m. (b) Variogram map by polarmap.m (with flipud command).
4. Discussion Software used to plot variogram maps normally use Cartesian coordinates. Such codes are often lengthy and
Christakos, G., Bogaert, P., Serre, M.L., 2002. Temporal GIS. Advanced Functions for Field-based Applications. Springer, New York, NY 217pp. Deutsch, C.V., Journel, A.G., 1998. GSLIB. Geostatistical Software Library and User0 s Guide. Oxford University Press, New York, NY 369pp. Englund, E., Sparks, A., 1988. GEO-EAS (Geostatistical Environmental Assessment Software) User0 s Guide. EPA/ 600/4-88/033a, v. 1.1. Mathworks, 2000. MATLAB User Guide. The MathWorks, Inc., Natick, MA. Pannatier, Y., 1996. VARIOWIN. Software for Spatial Data Analysis in 2D. Springer, New York, NY 91pp.
Computers & Geosciences 31 (2005) 801–803 www.elsevier.com/locate/cageo
Note
An alternative method for calculating variogram surfaces using polar coordinates$ F. Iwashita, R.C. Monteiro, P.M.B. Landim Geomathematic Laboratory, Department of Applied Geology, Sa˜o Paulo State University/Rio Claro campus, Av 24-A, 1515, Rio Claro 13506-900, Brazil Received 12 August 2004; received in revised form 15 December 2004; accepted 7 January 2005
1. Introduction
2. Variogram surface
Geostatistics deals with spatial and temporal data from natural phenomena and involves estimating values of regionalized variables. In this analysis the fundamental tool is the variogram that is a variability measure; for example, geological variability related to distance. The variogram shows the degree of spatial dependence between samples for a specific support. The squared differences between pairs of values are estimated assuming stationarity in increments. The semivariogram is calculated as
A variogram is a measure of spatial variability with respect to distance, but such variability can be different when taking into account different directions. This occurs with anisotropic phenomena. In this situation it is important to calculate variograms in different directions. This procedure may be carried out by constructing a variogram surface that contemplates all the possible directions of variability. Examples of variogram surfaces are shown in Pannatier (1996) and Deutsch and Journel (1998). This note presents an alternative approach to obtaining a variogram surface using polar coordinates. To implement this routine it is necessary to convert gðhÞ from Cartesian to polar coordinates to simplify programming: ( x ¼ h cos ðyÞ; y ¼ h sen ðyÞ;
gðhÞ ¼
1 XN h ½Zðxi þ hÞ Zðxi Þ2 , i¼1 2N h
where h ¼ basic lag, Nh ¼ number of pairs. Variogram analysis deals with the experimental variogram calculated from the data and the variogram model fitted to the data. The variogram model is chosen from a set of mathematical functions that describe the spatial relationship. This step involves the knowledge and researcher’s experience concerning the phenomenon and the study area. A good adjustment is necessary for achieving a good result.
$ Code on server at http://www.iamg.org/CGEditor/ index.htm. Corresponding author. Tel.: +55 19 3526 2811; fax: +55 19 3524 2445. E-mail address: [email protected] (F. Iwashita).
where x and y are Cartesian coordinates, h (lag) is the radius and the angle y is the interval (angular tolerance) to each direction between 901 and 901.
3. Program description The routine polarmap.m, running in MATLAB (Mathworks, 2000), is part of a software package for exploratory and geostatistical analysis developed at the
0098-3004/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2005.01.009
ARTICLE IN PRESS 802
F. Iwashita et al. / Computers & Geosciences 31 (2005) 801–803
Geomathematics Laboratory in the Sa˜o Paulo State University/Rio Claro campus. The whole package has other routines for calculating histograms, scatter plots, box-whisker plots, normal probability plots and linear regression. The calculations for the experimental variogram and ordinary kriging are made by the functions vario.m and kriging.m from the Bayesian Maximum Entropy Library/BMELib package (Christakos et al., 2002). These routines are the basis for calculating and plotting indicator variograms and indicator kriging. The package allows calculation and display of standard deviation maps for ordinary and indicator kriging and a specific routine for crossvalidation. An index for experimental variogram fitness to the theoretical variogram is proposed, based on the w2 (chi-square) test: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN ðgðhÞ g ðhÞÞ2 g ðhÞ j¼1 , i¼ N where i is the calculated adjusted index; smaller the index, better the fitting; N is the number of elements in g domain; gðhÞ are the experimental variogram values; g ðhÞ are the theoretical model values for the fitted variogram. The BMELIB routine vario.m produces directional variograms (901 a 901) with angular tolerance defined by the user. The polarmap() routine calculates n (angular tolerance /180) directional variograms. The gðhÞ values calculated by each directional variogram are stored in a matrix (escalaC) in vector form. A polar mesh is created based on angular tolerance defined by the user, and the routine pcolor( ) plots the matrix values in escalaC as an index for a color scale map.
Table 1 Extract from example input file in GeoEAS format
3.1. Input data
Fig. 1. Variogram map by GSLIB (Deutsch and Journel, 1998).
Example. dat—Geostatistical Assessment Software 5 Easting Northing Arsenic Cadmium Lead 288.0 311.0 .850 285.6 288.0 .630 273.6 269.0 1.02 280.8 249.0 1.02 273.6 231.0 1.01 ^ ^ ^ 465.6 216.0 .930 492.0 216.0 .750 345.6 216.0 1.45
Environmental
11.5 8.50 7.00 10.7 11.2 ^ 11.6 6.90 9.90
18.25 30.25 20.00 19.25 151.5 ^ 25.00 33.00 40.75
Cadmium: Variogram map 162.500
Coordinates N-S
30.000 28.000 26.000 24.000 22.000 20.000 18.000 16.000 14.000 12.000 10.000 8.000 6.000 4.000 2.000 0.0
-162.500 -162.500
Coordinates E-W
162.500
The function readGeoeas.m (BMELib package) reads and transforms the database to the vector and matrix structures used by the analysis. The datafile is in GeoEAS format (ASCII) to be read, preferentially using a ‘.dat’ extension. The chosen example comes from Englund and Sparks (1988), using the variable cadmium (Table 1). 3.2. Example To compare the conversion by the polarmap.m routine, three variogram surfaces are presented. The first one was calculated by GSLib (Fig. 1), the second by Variowin (Fig. 2) and the third obtained by polarmap.m routine (Fig. 3). The datasets are the same, but there are some differences in the variogram maps due to inherent differences between Cartesian and polar coordinates and slightly different scales.
Fig. 2. Variogram map by Variowin (Pannatier, 1996).
ARTICLE IN PRESS F. Iwashita et al. / Computers & Geosciences 31 (2005) 801–803
803
involved, but the routine developed here for polar coordinates in MATLAB has just 49 lines. Its simplicity is related to the use of MATLAB toolboxes and their simple and versatile commands for matrix operations. The programming is straightforward, because directional variograms are calculated based on angles and not on Cartesian coordinates, and are more easily implemented than the calculations using Cartesian coordinates. The resulting variogram surfaces are effective representations of the spatial variability and anisotropic properties. This routine and the whole package can be downloaded from the IAMG server or obtained from the author at E-mail to: [email protected].
Acknowledgements To ‘‘The State of Sa˜o Paulo Research Foundation/ FAPESP’’ for a scholarship to the first author of this paper.
References
Fig. 3. (a) Variogram map by polarmap.m. (b) Variogram map by polarmap.m (with flipud command).
4. Discussion Software used to plot variogram maps normally use Cartesian coordinates. Such codes are often lengthy and
Christakos, G., Bogaert, P., Serre, M.L., 2002. Temporal GIS. Advanced Functions for Field-based Applications. Springer, New York, NY 217pp. Deutsch, C.V., Journel, A.G., 1998. GSLIB. Geostatistical Software Library and User0 s Guide. Oxford University Press, New York, NY 369pp. Englund, E., Sparks, A., 1988. GEO-EAS (Geostatistical Environmental Assessment Software) User0 s Guide. EPA/ 600/4-88/033a, v. 1.1. Mathworks, 2000. MATLAB User Guide. The MathWorks, Inc., Natick, MA. Pannatier, Y., 1996. VARIOWIN. Software for Spatial Data Analysis in 2D. Springer, New York, NY 91pp.