Spatial Prediction of Maximum Temperature Across Sri Lanka Using Kriging Interpolation
Author
H.L.D.H. Sandamini (AS2022636)
Published
September 6, 2026
1 Introduction
This project develops a spatial prediction model to estimate maximum temperature across Sri Lanka using kriging interpolation. Daily maximum temperature records from 24 weather stations for June 2026 were first explored and quality-checked, then aggregated to a single monthly value per station, before being used as input to the spatial interpolation model. Because the objective is spatial interpolation, the temporal dimension is treated as fixed by analysing only the monthly mean maximum temperature for June 2026. The aim of this project is to compare kriging-based predictions with alternative interpolation approaches such as Inverse Distance Weighting (IDW) method and evaluate their predictive performance.
2 Data Sources
The two original source files for this analysis are provided alongside this report and can be viewed or downloaded directly:
Daily maximum temperature values for June 2026 were obtained from the daily weather bulletin images published by the Department of Meteorology, Sri Lanka, on their official Facebook page. One bulletin image was posted per day, giving 30 images in total covering the full month. These images are compiled together in the supplementary PDF daily_bulletins_june2026.pdf.
Manual data entry note
Because the source bulletins were published as images rather than in a machine-readable format, the daily maximum temperature value for each station was manually read off each image and entered into a spreadsheet (maxtemp.xlsx). Every value was cross-checked against the original bulletin image at the time of entry.
3 Exploratory Data Analysis
3.1 Descriptive statistics by station
Daily maximum temperature was summarised by station to check data quality and inspect the shape of each station’s distribution before aggregation.
Table 1: Descriptive statistics of daily maximum temperature by station, June 2026
Location
n_days
Mean
Median
SD
Min
Q1
Q3
Max
IQR
Range
Skewness
Kurtosis
Anuradhapura
30
32.74
32.70
0.82
31.0
32.42
33.48
33.9
1.05
2.9
-0.48
-0.41
Badulla
30
31.46
31.75
1.16
28.4
31.20
32.20
33.2
1.00
4.8
-1.10
0.63
Bandarawela
30
27.49
27.60
1.28
23.0
26.80
28.48
29.8
1.68
6.8
-1.31
3.04
Batticaloa
30
34.34
34.05
1.04
32.3
33.80
34.60
36.8
0.80
4.5
0.81
0.44
Colombo
30
30.91
31.05
0.74
28.7
30.72
31.40
31.8
0.68
3.1
-1.37
1.40
Galle
30
29.94
30.00
0.37
29.2
29.70
30.28
30.6
0.58
1.4
-0.23
-0.98
Hambanthota
30
32.52
31.90
1.56
30.1
31.35
33.67
35.8
2.32
5.7
0.54
-0.90
Jaffna
30
32.91
32.70
0.77
31.8
32.30
33.50
34.2
1.20
2.4
0.26
-1.41
Katugasthota
30
29.12
29.20
1.44
24.8
28.33
29.98
31.8
1.65
7.0
-0.57
0.96
Katunayake
30
31.04
31.15
0.63
28.8
30.83
31.40
31.7
0.57
2.9
-1.70
3.32
Kurunagala
30
31.46
31.70
1.14
28.7
30.80
32.15
33.3
1.35
4.6
-0.66
-0.16
Maha Illuppallama
30
31.89
31.80
0.81
30.1
31.60
32.55
33.6
0.95
3.5
-0.28
-0.11
Mannar
30
31.49
31.50
0.28
30.8
31.40
31.70
32.2
0.30
1.4
-0.09
0.52
Mattala
30
33.78
34.05
1.30
31.0
32.73
34.88
35.8
2.15
4.8
-0.35
-1.03
Monaragala
30
34.87
34.70
1.01
32.7
34.40
35.68
36.8
1.28
4.1
-0.14
-0.63
Mullativu
30
36.01
35.90
0.62
34.9
35.50
36.60
37.1
1.10
2.2
0.04
-1.30
Nuwara Eliya
30
19.88
19.30
1.65
16.3
18.80
21.45
23.1
2.65
6.8
0.13
-0.91
Polonnaruwa
30
35.76
35.65
0.91
34.2
35.20
36.58
37.2
1.38
3.0
0.00
-1.13
Pothuvil
30
35.93
36.15
1.24
32.9
34.95
36.75
38.4
1.80
5.5
-0.31
-0.46
Puttalam
30
31.71
31.75
0.50
30.1
31.50
32.00
32.6
0.50
2.5
-1.00
1.68
Rathmalana
30
31.95
32.15
1.12
27.7
31.60
32.58
33.3
0.98
5.6
-1.96
4.71
Rathnapura
30
31.35
31.60
1.40
27.5
31.02
32.20
33.4
1.18
5.9
-1.19
1.12
Trincomalee
30
35.33
35.50
0.86
33.0
34.80
35.88
36.8
1.08
3.8
-0.55
0.04
Vavuniya
30
34.18
34.45
1.72
26.0
33.82
35.08
35.8
1.25
9.8
-3.50
13.99
All 24 stations report a complete 30-day record with no missing values. Some stations show noticeable skewness (e.g. Vavuniya, Rathmalana, Kurunagala), reflecting isolated unusually hot or cool days rather than data errors, since all recorded values fall within a physically plausible range for Sri Lanka (16.3-38.4°C).
3.2 Distribution of daily maximum temperature by station
daily_data <- daily_data %>%mutate(Location =reorder(Location, Max_Temperature, FUN = median))ggplot(daily_data, aes(x = Location, y = Max_Temperature)) +geom_boxplot(fill ="lightblue", outlier.color ="red", outlier.shape =16)+coord_flip() +labs(title ="Distribution of Daily Maximum Temperature by Station (June 2026)",x ="Station",y ="Max Temperature (°C)" ) +theme_minimal(base_size =12)
Figure 1: Distribution of daily maximum temperature by station, ordered by median (June 2026)
Figure 1 shows a clear spatial variation in daily maximum temperature: hill country stations such as Nuwara Eliya and Bandarawela are consistently coolest, while dry-zone stations such as Mullativu, Pothuvil, and Polonnaruwa are consistently hottest. Several stations (e.g. Batticaloa, Rathnapura, Katugasthota, Bandarawela) show one or more outlier days, consistent with the skewness values reported in Table 1 above.
4 Temporal Aggregation: Daily to Monthly
Daily maximum temperature values were aggregated to a single monthly value per station using the arithmetic mean, consistent with standard climatological convention (e.g. WMO, CRU, WorldClim), where “monthly mean maximum temperature” is defined as the average of daily Tmax observations over the month.
Although the station-level descriptive statistics above revealed some skewness and boxplot outliers within individual stations, these deviations reflect genuine day-to-day weather variability (e.g. isolated hot or rainy days) rather than data entry errors, and are physically plausible given the observed island-wide range of daily maximum temperatures. Since kriging operates on the spatial distribution of the monthly station values rather than the within-station daily distribution, station-level daily skewness does not violate the assumptions underlying the spatial interpolation model. The arithmetic mean was therefore kept as it provides a standard, physically interpretable summary of average thermal conditions at each station, and comparison against the median confirmed the two measures were reasonably close for most stations, supporting the mean as a representative choice.
4.1 Converting daily data to a monthly value
Since the daily mean and median agree closely for most stations, the daily mean was used to derive each station’s monthly value.
monthly_station_data <- station_summary %>%select(Location, Mean_Tmax = Mean)kable(monthly_station_data, digits =2,caption ="Monthly mean maximum temperature by station, June 2026")
Monthly mean maximum temperature by station, June 2026
The aim of this study is to examine the spatial distribution of maximum temperature across Sri Lanka rather than changes over time. Therefore, the analysis uses the monthly mean maximum temperature for June 2026, with one value representing each weather station. This means that all observations are assumed to belong to the same time period. By averaging the daily temperatures into a monthly value, day-to-day variations are removed, allowing the analysis to focus only on the spatial variation of maximum temperature across the country.
5 Station Coordinates
Kriging requires a spatial coordinate (latitude/longitude) for every station, in addition to the monthly temperature value derived above.
5.1 Source
Station coordinates were obtained from the WMO OSCAR/Surface database (https://oscar.wmo.int/surface/), the World Meteorological Organization’s official repository of WIGOS metadata for all internationally registered surface weather stations. Sri Lanka was selected as the country filter, and the resulting station list was exported as a CSV (StationSearchResults.csv), which is provided alongside this report:
The raw OSCAR export contained 26 records for Sri Lanka, two more than the 24 stations present in the temperature dataset. The following cleaning steps were applied to produce a coordinate file that lines up exactly with the 24 stations in maxtemp.xlsx:
Removed 2 non-matching records: KANKASANTURAI (a separate station not present in the temperature dataset) and an underwater ARGO float record (2902365), which is not relevant to this analysis.
Renamed KANDY to Katugasthota: the OSCAR station registered under the name “Kandy” is physically sited in Katugasthota (a suburb of Kandy), which matches the station name used in the temperature dataset.
Station names were matched to the exact spelling used in maxtemp.xlsx (e.g. Hambanthota, Rathmalana, Rathnapura, Mullativu) so the two datasets can be joined directly by Location.
The cleaned result is a 24-row coordinate table, saved as cleaned_station_coordinates.xlsx, provided alongside this report. In addition to latitude and longitude, this file also includes each station’s elevation (also sourced from the OSCAR export), which is used later as a covariate for Regression Kriging.
This section fits Inverse Distance Weighting (IDW) approach and several kriging interpolation approaches to the monthly station data, then compares their predictive performance using cross-validation.
6.1 Converting Station Coordinates to Spatial Objects
stars object with 2 dimensions and 2 attributes
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
var1.pred 19.97645 31.0146 32.20156 32.03726 33.36496 36.00989 1697
var1.var NA NA NA NaN NA NA 4350
dimension(s):
from to offset delta refsys x/y
x 1 50 352448 5000 WGS 84 / UTM zone 44N [x]
y 1 87 1086668 -5000 WGS 84 / UTM zone 44N [y]
ggplot() +geom_stars(data = idw.result, aes(fill = var1.pred, x = x, y = y)) +geom_sf(data =st_cast(map, "MULTILINESTRING")) +geom_sf(data = stations.sf) +scale_fill_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()
Figure 3: IDW prediction of monthly mean maximum temperature
6.4.1 Comparison
library(patchwork)p1 <-ggplot(map) +geom_sf() +geom_sf(data = stations.sf, aes(color = Mean_Tmax), size =3) +scale_color_distiller(palette ="RdBu", direction =-1) +theme_minimal()p2 <-ggplot() +geom_stars(data = idw.result, aes(fill = var1.pred, x = x, y = y)) +geom_sf(data =st_cast(map, "MULTILINESTRING")) +geom_sf(data = stations.sf) +scale_fill_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()p1|p2
Figure 4: Empirical semivariogram of monthly mean maximum temperature
6.5.1 Model Selection Criteria
To determine the most appropriate model for the spatial analysis, several theoretical models available in the R gstat package were compared (Figure 5).
show.vgms(par.strip.text =list(cex =0.75))
Figure 5: Theoretical variogram models available within the gstat package framework.
Although the gstat package provides many theoretical variogram models, this study considered only the Exponential ("Exp"), Spherical ("Sph"), and Gaussian ("Gau") models. These three models were selected based on the following reasons:
Realistic representation: Linear ("Lin") and Power ("Pow") models assume that spatial variation increases continuously with distance, which is not suitable for Sri Lanka’s limited geographical area. Similarly, the Wave ("Wav") model represents cyclic spatial patterns that do not match the distribution of climatic variables in this study.
Limited sample size: The Matérn ("Mat") model requires estimating an additional smoothness parameter. With data from only 24 weather stations, estimating this parameter can lead to overfitting and difficulties in model convergence.
Avoiding unnecessary models: The Circular ("Cir") and Pentaspherical ("Pen") models produce results that are very similar to the Spherical model. Therefore, including them would provide little additional benefit.
Isotropic Assumption: We assume that temperature differences only depend on the distance between weather stations, not the direction they are facing. Even though Sri Lanka’s monsoons and coastlines can cause weather to move in specific directions, we only have 24 weather stations in total. If we tried to build a complicated model that accounts for direction (anisotropic model), we would have to split our tiny dataset into even smaller directional groups. This would create a lot of statistical noise, make the code unstable, and likely crash our model fitting. Keeping the model isotropic is the safest and most reliable choice for a small sample size like ours.
v.m <-fit.variogram(v1, vgm(psill =10, model =c("Exp", "Sph", "Gau"), range =50000, nugget =1))plot(v1, v.m)
The optimized parameter values obtained for the automatically selected model are presented below:
# Optimized theoretical variogram parameters selected by gstatv.m
model psill range
1 Gau 15.1274 20847.33
The Gau model was selected as the best-performing model because it produced the smallest weighted residual sum of squares during variogram fitting.
Since only 24 weather stations were available, the empirical variogram may contain some variability. Therefore, the fitted variogram model should be interpreted with caution, as this limitation is mainly due to the small sample size rather than the model fitting procedure.
ggplot() +geom_stars(data = krigSK, aes(fill = var1.pred, x = x, y = y)) +geom_sf(data =st_cast(map, "MULTILINESTRING")) +geom_sf(data = stations.sf) +scale_fill_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()
Figure 6: Simple kriging prediction of monthly mean maximum temperature
Simple kriging assumes a known, fixed mean (beta), rather than estimating it from the data.
6.7 Ordinary Kriging
krigOK <-krige(Mean_Tmax ~1, stations.sf, grid, v.m, debug.level =0)ggplot() +geom_stars(data = krigOK, aes(fill = var1.pred, x = x, y = y)) +xlab(NULL) +ylab(NULL) +geom_sf(data =st_cast(map, "MULTILINESTRING")) +geom_sf(data = stations.sf) +coord_sf(lims_method ="geometry_bbox") +scale_fill_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()
Figure 7: Ordinary kriging prediction of monthly mean maximum temperature
Ordinary kriging estimates the mean from the data itself, rather than assuming it is known in advance, and is the standard baseline kriging approach used for comparison in this project.
6.8 Universal Kriging
v.uk <-variogram(Mean_Tmax ~ X + Y, stations.sf)v.m.uk <-fit.variogram(v.uk, vgm(psill =10, model =c("Exp", "Sph", "Gau"),range =50000, nugget =1))grid_xy <- grid |>st_as_sf(as_points =TRUE) %>%mutate(X =st_coordinates(.)[,1], Y =st_coordinates(.)[,2])krigUK <-krige(Mean_Tmax ~ X + Y, stations.sf, grid_xy, v.m.uk, debug.level =0)
ggplot() +geom_sf(data = map) +geom_sf(data = krigUK, aes(color = var1.pred)) +geom_sf(data = stations.sf) +scale_color_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()
Figure 8: Universal kriging prediction of monthly mean maximum temperature
Universal kriging models a spatial trend directly from the coordinates (here, a first-order/linear trend surface) rather than assuming a constant mean across the island.
6.9 Regression Kriging
v.rk <-variogram(Mean_Tmax ~ Elevation, stations.sf)v.m.rk <-fit.variogram(v.rk, vgm(psill =10, model =c("Exp", "Sph", "Gau"), range =50000, nugget =1))plot(v.rk, v.m.rk)
Figure 9: Variogram of monthly mean maximum temperature using elevation as a covariate
The optimized parameter values obtained for the automatically selected model are presented below:
# Optimized theoretical variogram parameters selected by gstatv.m.rk
model psill range
1 Gau 2.688574 19933.85
The optimization procedure selected the Gau model as the best-performing model.
Regression kriging needs the elevation covariate at every prediction grid cell, not only at the 24 stations, so a digital elevation model (DEM) covering Sri Lanka was downloaded to provide this.
Figure 10: Regression kriging prediction of monthly mean maximum temperature, using elevation as a covariate
Regression kriging incorporates elevation as a covariate, consistent with the physical expectation that temperature decreases with elevation (visible in the hill-country stations Nuwara Eliya and Bandarawela being the coolest in the dataset).
6.10 Stratified Kriging
Sri Lanka is conventionally divided into three climatic zones (Wet, Dry, and Intermediate). Stratified kriging fits a separate variogram and kriging model within each zone, rather than a single model for the whole island.
6.10.1 Basis for the zone classification
The three-zone classification used in this study follows the standard climatic zones of Sri Lanka, which are mainly defined based on average annual rainfall (Panabokke, 1996; Department of Meteorology, Sri Lanka). The Wet Zone, located mainly in the south-western region and central highlands, receives more than 2,500 mm of annual rainfall and does not experience a clear dry season. The Dry Zone, covering the northern, eastern, and south-eastern lowlands, receives less than 1,750 mm of annual rainfall and has a distinct dry period from May to September. The Intermediate Zone lies between these two regions and receives approximately 1,750–2,500 mm of annual rainfall.
In this study, the 24 weather stations were categorized into these climatic zones based on their known geographical locations. For example, stations such as Colombo and Galle were assigned to the Wet Zone due to their location in the south-western region, while stations such as Jaffna, Mannar, and Anuradhapura were classified under the Dry Zone. The classification was based on general geographic positioning of the stations rather than using exact station coordinates with an official digital climatic zone boundary map.
Limitation
A few districts in Sri Lanka extend across more than one climatic zone, meaning they do not belong entirely to a single zone. Previous classifications have identified areas such as Badulla, Kurunegala, and Hambantota as districts with mixed climatic characteristics, including combinations of dry, intermediate, and wet conditions. Therefore, the assignments used in this study, where Badulla was classified as Intermediate Zone, Kurunegala as Intermediate Zone, and Hambantota as Dry Zone, should be considered reasonable approximations rather than exact classifications. This uncertainty may have influenced the Stratified Kriging results, together with the limitation caused by the relatively small number of stations available within each climatic zone.
With only 24 stations split across 3 zones (roughly 5–13 stations per zone per zone), each zone’s variogram is fitted on very few pairs and should be treated as illustrative rather than reliable. This section demonstrates the method, but the resulting map is more uncertain than the island-wide methods above.
To produce a single map, each grid cell is first assigned to the zone of its nearest station, then predicted using that zone’s own kriging model.
ggplot() +geom_stars(data = krigBK, aes(fill = var1.pred, x = x, y = y)) +geom_sf(data =st_cast(map, "MULTILINESTRING")) +geom_sf(data = stations.sf) +scale_fill_distiller(palette ="RdBu", direction =-1, name ="Mean_Tmax.pred (°C)")+theme_minimal()
Figure 12: Block kriging prediction of monthly mean maximum temperature (5km x 5km block averages)
Block kriging predicts the average value over a 5km x 5km block centred at each grid point, rather than a point estimate, typically producing smoother predictions with lower prediction variance than point kriging.
Table 3: Leave-one-out cross-validation performance by method
Leave-one-out cross-validation performance by method
Method
RMSE
MAE
Sample Size
Simple Kriging
3.113
2.125
24
Ordinary Kriging
3.207
2.202
24
Universal Kriging
2.528
1.609
24
Regression Kriging
1.990
1.720
24
Block Kriging
3.205
2.201
24
Stratified Kriging
2.728
1.717
24
IDW
2.988
2.029
24
The cross-validation results show that Regression Kriging and Universal Kriging performed better than the other interpolation methods. Regression Kriging obtained the lowest RMSE value (1.990), which indicates that it had the best overall prediction accuracy among the tested models. Meanwhile, Universal Kriging recorded the lowest MAE value (1.609), meaning that its average prediction errors were slightly smaller across the weather stations. These results suggest that including additional information, such as spatial trends in Universal Kriging or elevation as a predictor in Regression Kriging, can improve temperature predictions compared with Simple Kriging, Ordinary Kriging, Block Kriging, Stratified Kriging, and IDW.
Based on the results, Regression Kriging can be considered the most suitable method for this dataset because it achieved the lowest RMSE and also accounts for elevation, which is an important factor affecting temperature differences across Sri Lanka. However, the difference between Regression Kriging and Universal Kriging should be interpreted carefully because the analysis was performed using only 24 weather stations. A limited number of stations can affect variogram estimation and may cause the model parameters to be influenced by individual observations. Therefore, these cross-validation results provide a useful comparison between methods, but further data from more weather stations would be needed to confirm the accuracy of the interpolation results across the entire country.
