IJERD – International Journal of Environmental and Rural Development (2021) 12-1

Research article
erd

Spatial Variability and Mapping of Soil Properties Using GISBased Geostatistic in Myanmar
TUN TUN HLAING*
Dept. of Soil and Water Science, Yezin Agricultural University, Myanmar
Email: thlaing83@gmail.com

SWE SWE MAR
Dept. of Soil and Water Science, Yezin Agricultural University, Myanmar

CHAW SU LWIN
Dept. of Soil and Water Science, Yezin Agricultural University, Myanmar

HTAY HTAY OO
Dept. of Agronomy, Yezin Agricultural University, Myanmar

KYAW NGWE
Dept. of Soil and Water Science, Yezin Agricultural University, Myanmar
Received 29 December 2020

Accepted 16 May 2021(*Corresponding Author)

Abstract Precise information on the spatial variability of soil is a crucial component for
productive intensive agriculture, sustainable development, and the management of natural
resources. The primary purpose of the study was to investigate the spatial variability of soil
properties of the study site at Yezin Agricultural University Field, Myanmar using geostatistics.
A total of 94 composite soil samples were collected from a depth of 0 to 20 cm, in a systematic
grid (50 x 50 m2) at the site in May 2019. Soil pH, electrical conductivity (EC), soil organic
matter (SOM), total soil nitrogen (TSN), available phosphorus (Ava-P), and available
potassium (Ava-K) were measured using standard analytical methods. Data were analyzed
geostatistically based on semivariogram. The exponential model best fitted the semivariogramfor pH, EC, Ava-P, and Ava-K; SOM was adapted from the Gaussian model while TSN was
adapted from the spherical model. The nugget/sill ratio showed a strong spatial dependence
exists for EC, Ava-P, and Ava-K and a moderate spatial dependence for pH, SOM, and TSN.
Most of the soil was found to be strongly acidic. It was also found that EC, SOM, and Ava-P
are very low in most of the study area. Most of the study area was found to have low TSN
levels, while Ava-K content was low over the entire area. With such an analysis, it is possible toplan better nutrient management practices for agricultural production and environmental
protection. Therefore, geostatistical analysis with ordinary kriging is a useful tool for studying
the spatial variability of soil properties.
Keywords spatial variability, soil properties, geostatistic, semivariogram, kriging

INTRODUCTION
An understanding of the distribution of soil properties is essential for ecological modelling,
environmental predictions, precise agriculture, and management of natural resources (Wang et al.,2009). However, soil properties vary spatially from a small to a larger regional area, and are affectedby intrinsic (parent materials and climate) and extrinsic factors (soil management practices, fertility
status, crop rotation) (Cambardella and Karlen, 1999). Therefore, demands for more accurateⒸ ISERD

15

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

information on spatial variability of soils are significant for intensive agriculture, sustainable
development, and natural resource management (Karlen et al., 2011).
Many studies have used a classical statistical method to quantify spatial characteristics in soil
properties (Salehi et al., 2013). However, physico-chemical characteristics of soil often exhibit spatial
dependency, which cannot be recorded with classical statistical methods (Lin et al., 2005). To
overcome this problem, many researchers apply geostatistical interpolation methods to estimate thespatial variability of soil properties (Cambardella et al., 1994; Webster and Oliver, 2007).
Geostatistics is a set of statistical tools that can be used to investigate and predict the spatial
structure of georeference variables and generate soil property maps (Patil et al., 2011). Based on thegeostatistical analysis, several studies have been carried out to characterize the spatial variability ofvarious soil properties (Weindorf and Zhu, 2010). Among the various geostatistical methods, ordinary
kriging is widely used to map spatial variations in soil fertility because it offers a higher level ofpredictive accuracy (Song et al., 2013). Therefore, it is important to extend the availability of soil
resource information maps to allow the planning of appropriate soil management practices, including
fertilization for agricultural production and environmental protection.
OBJECTIVES
The objective of this research is to investigate the spatial variability and mapping of spatial distributionof selected soil properties status in the study area using GIS-based geostatistical analysis.
METHODOLOGY
Study Area, Soil Sampling and Laboratory Analysis
The study area was Yezin Agricultural University Field (1949'47''1950'21''N and 9615'32''9616'15''E), Zeyarthiri Township, Nay Pyi Taw Union Territory, in central Myanmar. It has anelevation ranging from 121.547 m to 125.205 m above sea level (Fig 1). The area of the research site is
18.19 ha. The study area has an average temperature of 26.8°C and a mean annual rainfall of about
1420 mm. Summer and monsoon rice are the main crops in the study area and these include both rainfed and irrigated rice cultivations.

Fig. 1 Georeferenced sampling sites of the research area
A total of 94 composite soil samples were collected at a depth of 0 to 20 cm based on a systematic
grid (50 x 50 m2) determined with the help of a hand-held GPS device, in May 2019. Soil samples were
air-dried and ground so as to pass through a 2-mm sieve. Soil pH and electrical conductivity (EC) wereⒸ ISERD

16

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

measured in a 1:5 soil/water extract, soil organic matter (SOM) was determined by Heanes wet
oxidation method, total soil nitrogen (TSN) was analyzed by the Semi-micro Kjeldahl steam
distillation method, available phosphorus (Ava-P) was measured by the Olsen-P method, and available
potassium (Ava-K) was determined by extraction with 1M ammonium chloride (Rayment and Lyons,2011).
Geostatistical Analysis
The geostatistical analysis, consisting of semivariogram calculation, cross-validation, and mapping,
was performed with the Geostatistical Analyst Extension Tool of ArcGIS 10.7. The spatial variation ofsoil properties was analyzed using geographical semivariogram to quantify the spatial variation of aregionalized variables which derives important parameters used for ordinary kriging (OK) spatial
interpolation (Krige, 1951). The semivariogram is used as a fundamental tool to study the spatial
distribution structure of soil properties. The semivariogram analyzes were performed before theapplication of OK interpolation, as the semivariogram model determines the interpolation function(Goovaerts, 1997), defined as:2

 ( h) =1 N (h) z ( xi ) − z ( xi + h)
2 N (h) i =1

(1)

where,γ (h) is the experimental semivariogram value in a distance interval h (in meters or km), N (h)is the number of sample pairs that are located by a particular distance (h) from each other. z (xi) and z(xi+ h) are the values of a regionalized variable at location xi and xi+ h, respectively (Wang and Shao2013).
Theoretical semivariogram models fit the empirical semivariogram obtained from the data to
generate geostatistical parameters, including nugget variance (C0), structured variance (C), sill variance(C0+C), and distance parameters (A). The nugget/ sill ratio, C0/ (C0 + C) is calculated to characterize
the spatial dependency of the values. A nugget/ sill ratio is classified as strongly spatially dependent if
the ratio is less than or equal to 0.25, moderately spatially dependent if the ratio is between greater than
0.25 and less than or equal to 0.75 while it is classified as a weak spatial dependent if it is greater than
0.75 (Cambardella et al., 1994).
Several semivariogram models were evaluated to select the best fit with the data. The model,
spherical, Gaussian, or exponential that offers the best fit varies depending on the soil parameters(Ramzan and Wani, 2018). A cross-validation technique was used to evaluate and compare theperformance of the OK interpolation method. The lowest RMSE (Root Mean Square Error) value
indicates the best fit for the variogram model (Panday et al., 2018). The predictive maps of soil
properties are then created using a semivariogram model through OK.
Assessment of Accuracy of Interpolation Map
The effectiveness of interpolation was evaluated based on Goodness-of-Prediction Estimate (G) Eq. (2).
A “G” value of 100% indicates a perfect prediction, positive values (i.e., 0 to 100%) show that thepredictions are more reliable than the use of the sample mean, and negative values indicate that thepredictions are less reliable than using the sample means (Laekemariam et al., 2018).
N2  z ( xi ) − zˆ ( xi ) 
  100
G = 1 − i=1 N
2


ˆ


z
(
x
)
−
y

i


i =1

(2)

Ⓒ ISERD

17

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

Where, z(xi) is the observed value at location i, ž(xi) is the predicted value at location i, N is the sample
size, and ŷ is the sample mean.
RESULTS AND DISCUSSION
Geostatistical Analysis
The semivariogram parameters obtained from the best-fit model are in Table 1. An exponential model
produced the best fit to semivariogram for pH, EC, Ava-P, and Ava-K. This model is one of thestandard models used in the study of soil properties (Cambardella et al., 1994; Reza et al., 2016). The
spherical model was the best suited to the semivariogram of TSN, while a Gaussian model was the best
fit for SOM.
The range for all soil properties varies from 96.76 m to 276.68 m, and therefore the length of thespatial autocorrelation is much longer than the sampling interval of 50 m. According to Goovaerts(1997), the current sample design is appropriate for this study, and it is expected that the interpolated
map will display good spatial structure.
In the present study, the nugget/ sill ratio showed that EC, Ava-P, and Ava-K were strongly
dependent spatially whereas, pH, SOM, and TSN were moderately dependent spatially. The strong
spatial dependency suggests that intrinsic factors, such as climate, parent material, topography, soil
properties, and other natural factors, play important roles in spatial variability. The weak spatial
dependency indicates that the spatial variability is mainly caused by extrinsic factors, such as
fertilization, local farming practice, cropping systems, and other human activities. The moderatespatial dependency shows that spatial variability is caused by a mix of extrinsic and intrinsic factors(Bhunia et al., 2018; Cambardella et al., 1994).
The G-values are greater than zero for all soil parameters. This value indicated that spatial
prediction using semivariogram parameters is better than assuming that the mean of observed values isthe best value for an unsampled location. This result also shows that semivariogram parameters
obtained from fitting experimental semivariogram values describe the spatial variation reasonably(Reza et al., 2010).
Table 1 Geostatistical parameters of the fitted semivariogram models for soil properties
Parameters
Soil pH
EC
SOM
TSN
Ava-P
Ava-K

Model

NuggetSill

Exponential
Exponential
Gaussian
Spherical
Exponential
Exponential

0.10370.00000.09000.00050.08560.00000.1822
0.48220.28690.00080.6445214.0108Range(m)
253.3496.76
276.68227.30178.09108.58NuggetDSD/ Sill0.57 Moderate
0 Strong
0.31 Moderate
0.38 Moderate
0.13 Strong
0 Strong

RMSE

G (%)

0.41
0.03
0.34
0.03
6.71
12.26

6.11
13.25
47.43
25.08
25.66
39.90

DSD: Degree of Spatial Dependence, RMSE: Root Mean Square Error, G: Goodness of Prediction

Spatial Distribution of Soil Properties
The parameter of the fitted semivariogram models is used for OK to produce a spatial distribution mapof soil properties in the study area. The spatial distribution of soil properties such as pH, EC, SOM,
TSN, Ava-P, and Ava-K are shown in (Fig. 2a-e). The distribution of the predicted soil pH map (Fig.
2a) shows that 0.77%, 80.92%, and 18.31% of the study area were very strongly acid, strongly acid,and moderately acid, respectively. Most of the soils were strongly acid which may be caused by a
mixture of the nature of the soil mineralogy, the use of acidic fertilizers, low input of organic materials,Ⓒ ISERD

18

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

and removal of base nutrients (Rawal et al., 2018). The predicted map of EC (Fig. 2b) shows that89.22% and 10.78% of the soil in the study area, can be described as very low and low, respectively.
According to the soil guide (Moore, 2001); low EC levels only have a minimal impact on plant growth.
(a)

(b)

(c)

(d)

(e)

(f)

Fig. 2 Spatial distribution maps for (a) pH, (b) EC, (c) SOM, (d) TSN, (e) Ava-P and (f) Ava-K
The distribution of SOM (Fig. 2c) ranged from very low (62.73%) to low (37.27%), but low levels
were most prevalent. The lower organic matter content in these soils can be attributed to the poor
management practices such as intensive cropping, the complete removal of crop residues, and lack ofaddition of organic fertilizer sources (Gebreselassie, 2002). The majority of the soils were low(97.75%) in total nitrogen content, whereas 2.25% of the study area is rated at a medium level (Fig. 2d).
The low nitrogen levels in most areas could be as a result of continuous cereal-based cropping, lower
external organic-N inputs (like plant residues, animal manures), N (nitrate ions) leaching problem, andaddition of a low amounts of SOM (Patil et al., 2011). The Ava-P map (Fig. 2e) shows that, in terms ofⒸ ISERD

19

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

area coverage, 68.83%, 28.42%, and 2.75% of the study area has very low, low, and medium levels
respectively. This result indicates that, in general, the majority (97.25%) of the study area is deficientin phosphorus. The reason for low soil phosphorus levels may be due to the intensive cropping system,
low pH (acidic) soils, the imbalanced use of fertilizer, and nutrient mining (Sertsu and Ali, 1983). The
whole of the study area can be classified with a low rating for Ava-K content (Fig. 2f). The lowest
Ava-K in the study area might be due to the lowest SOM and the continuous removal of potassium by
cereal crops, as the field has been intensively cultivated for a long period.
CONCLUSION
Most of the soil was strongly acidic. The selected soil properties had a distribution showing low
content levels in most of the study area. The distribution of variability of the soil properties across thelandscape might be sufficient to construct fertility maps. The generated spatial distribution and fertility
maps can serve as a powerful tool for farmers, decision-makers, and planners to understand the existingsoil conditions and make sensible decisions to better manage the soil for sustainability and productivity.
These results show that geostatistical analysis with kriging is an effective tool for studying the spatial
variability of soil properties and that it will be useful technique for future soil sampling campaigns in
Myanmar.
ACKNOWLEDGEMENTS
It is a pleasure to express special gratitude to the supervisory committee and Australia Centre forInternational Agricultural Research (ACIAR) which supports research grant and technical support forthis study under the project title of “Management of nutrients for improved profitability andsustainability of crop production in central Burma (SMCN– 2014 -044)”.
REFERENCES
Bhunia, G.S., Shit, P.K. and Chattopadhyay, R. 2018. Assessment of spatial variability of soil properties using
geostatistical approach of lateritic soil (West Bengal, India). Annals of Agrarian Science, 16 (4), 436-443.
Cambardella, C.A. and Karlen, D.L. 1999. Spatial analysis of soil fertility parameters. Precision Agriculture, 1(1), 5-14.
Cambardella, C.A., Moorman, T.B., Parkin, T.B., Karlen, D.L., Novak, J.M., Turco, R.F. and Konopka, A.E.
1994. Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal,
58, 1501-1511.
Gebreselassie, Y. 2002. Selected chemical and physical characteristic of soils of Adet research center and its
testing sites in north-western Ethipoia. Ethiopian Journal of Natural Resources, 4 (2), 199-215.
Goovaerts, P. 1997. Geostatistics for natural resources evaluation. Oxford University Press.
Karlen, D.L., Wienhold, B.J., Kang, S., Zobeck, T.M. and Andrews, S.S. 2011. Indices for soil management
decisions, USDA-ARS/UNL Faculty, Paper, 1381.
Krige, D.G. 1951. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of
the Southern African, Institute of Mining and Metallurgy, 52 (6), 119-139.
Laekemariam, F., Kibret, K., Mamo, T. and Shiferaw, H. 2018. Accounting spatial variability of soil properties
and mapping fertilizer types using geostatistics in southern Ethiopia. Communications in Soil Science andPlant Analysis, 49 (1), 124–137.
Lin, H., Wheeler, D., Bell, J. and Wilding, L. 2005. Assessment of soil spatial variability at multiple scales.
Ecological Modelling, 182 (3), 271-290.
Moore, G.A. 2001. Soilguide (soil guide) : A handbook for understanding and managing agricultural soils.
Department of Agriculture and Food, Western Australia, Perth. Bulletin 4343.

Ⓒ ISERD

20

IJERD – International Journal of Environmental and Rural Development (2021) 12-1

Panday, D., Maharjan, B., Chalise, D., Shrestha, R.K. and Twanabasu, B. 2018. Digital soil mapping in the Bara
district of Nepal using kriging tool in ArcGIS. Plos One, 13 (10), e0206350.
Patil, S.S., Patil, V.C. and Gaadi, K.A.I.A. 2011. Spatial variability in fertility status of surface soils. World
Applied Sciences Journal, 14 (7), 1020-1024.
Ramzan, S. and Wani, M.A. 2018. Geographic information system and geostatistical techniques to characterizespatial variability of soil micronutrients including toxic metals in an agricultural farm. Communications inSoil Science and Plant Analysis, 49 (4), 463-477.
Rawal, N., Acharya, K.K., Bam, C.R. and Acharya, K. 2018. Soil fertility mapping of different VDCs of Sunsari
District, Nepal using GIS. International Journal of Applied Science, Biotechnology, 6 (2), 142-151.
Rayment, G.E. and Lyons, D.J. 2011. Soil chemical methods: Australasia, 3. CSIRO publishing.
Reza, S.K., Baruah, U., Sarkar, D. and Singh, S.K. 2016. Spatial variability of soil properties using geostatisticalmethod: A case study of lower Brahmaputra plains, India. Arabian Journal of Geosciences, 9 (6), 446.
Reza, S.K., Sarkari, D., Daruah, U. and Das, T.H. 2010. Evaluation and comparison of ordinary kriging andinverse distance weighting methods for prediction of spatial variability of some chemical parameters of
Dhalai district. Tripura, 12.
Salehi, M., Safaei, Z., Esfandiarpour-Borujeni, I. and Mohammadi, J. 2013. Generalisation of continuous models
to estimate soil characteristics into similar delineations of a detailed soil map. Soil Research, 51 (4), 350-361.
Sertsu, S. and Ali, A. 1983. Phosphorus sorption characteristics of some ethiopian soils. Ethiopian Journal of
Agricultural Sciences.
Song, G., Zhang, L., Wang, K. and Fang, M. 2013. Spatial simulation of soil attribute based on principles of soil
science. 21st International Conference on Geoinformatics, 20-22 June 2013. Kaifeng, China.
Wang, Y., Zhang, X. and Huang, C. 2009. Spatial variability of soil total nitrogen and soil total phosphorus under
different land uses in a small watershed on the Loess Plateau, China. Geoderma, 150, 141-149.
Wang, Y.Q. and Shao, M.A. 2013. Spatial variability of soil physical properties in a region of the Loess Plateau
of PR China subject to wind and water erosion. Land Degrad Dev., 24 (3), 296-304.
Webster, R. and Oliver, M.A. 2007. Geostatistics for environmental scientists. John Wiley & Sons.
Weindorf, D.C. and Zhu, Y. 2010. Spatial variability of soil properties at Capulin volcano, New Mexico, USA:
Implications for sampling strategy. Pedosphere, 20 (2), 185-197.

Ⓒ ISERD

21