Kriging of road surface temperatures: data preparation

28 May, 2025

Kriging for Road Temperature Interpolation

As discussed in my previous post, I am exploring kriging methods as a primary approach to interpolate road temperatures to other locations.

Using the simple application available here, I obtained the following data points:

999901,dummy station,10.056845696219929,56.30947843428772
999902,dummy station,9.889097380698168,56.37031342981681
999903,dummy station,9.953001500896933,56.43503151016691
999904,dummy station,9.94022067685718,56.494572144089
999905,dummy station,10.197434760657213,56.521753737836036
999906,dummy station,10.154299479523047,56.48162852801897
999907,dummy station,10.044064872180176,56.374196514637816
999908,dummy station,10.310864574010024,56.34830928249778
999909,dummy station,10.162287494547893,56.45833001909294
999910,dummy station,10.013710415085761,56.52563682265704
999911,dummy station,10.117554610408757,56.45703565748594
999912,dummy station,10.031284048140423,56.4272653405249

These points, marked by "X"s, are plotted in the following figure. The red dots indicate locations where temperatures are known. selected_sample

Kriging is a robust spatial interpolation technique commonly used in geosciences to predict values at unmeasured locations based on the spatial autocorrelation of measured points. For a comprehensive introduction and relevant references, refer to this Medium article. A significant advantage of kriging is its ability to quantify the uncertainty associated with predictions.

The simplest form of kriging is Ordinary Kriging (SK), but a more complex algorithm is Universal Kriging (UK). Simple Kriging assumes a known and constant mean across the spatial domain, incorporated directly into the kriging equations. Universal Kriging, however, assumes an unknown mean that varies spatially, modeled as a function of spatial coordinates (often a polynomial). Universal Kriging can incorporate covariates such as elevation data from a Digital Surface Model (DSM).

Simple Kriging is suitable when there is a strong understanding of the underlying process and a constant mean can be confidently specified. Universal Kriging is more appropriate when there is a trend or drift in the data. Given that road temperatures vary across the domain, Universal Kriging appears to be the better option. Furthermore, we have access to a detailed Digital Surface Model (DSM) of Denmark, allowing us to incorporate terrain orography into the UK modeling.

Before performing the kriging calculation, we need to compute geometric quantities that characterize the terrain's orography, which will be used as covariates in the UK model. We will use GIS utilities to read the TIF data from the DSM and calculate terrain elevation, slope, and slope aspect at the station locations. The procedure involves:

Converting the latitude and longitude coordinates to UTM coordinates (since the DSM data is in UTM coordinates).
Identifying the TIF tiles containing the coordinates of interest.
Extracting terrain elevation, slope, and slope aspect at the selected locations.

I utilized the following GIS libraries to accomplish these tasks:

rasterio for loading the TIF files from the DSM.
scipy.ndimage for calculating geometrical features on the data points.
skimage.morphology to calculate minima and maxima around the geographical features.
gdal for extracting boundaries and extent of domain from the TIF files.

These steps are implemented using the following scripts:

extract_metadata_tif_files.sh extracts metadata from the TIF files.
identify_tiles_needed.py identifies the specific TIF files needed for each station.
terrain_metrics.py calculates the orographic features from the terrain at each station location.

The following figure illustrates the orography around station 999905, showing a slice of the terrain height in the North-South direction:

terrain_profile_station_999905

The line orientation and terrain features are depicted in the image below. The 500 m line provides a sense of scale. Each TIF tile from the Danish Surface Model covers a 1 km x 1 km square. The profile extends along the same distance but may appear truncated if the station is near the edge of the domain.

terrain_height_station_999905

With this information, we can derive the following orographic features for the selected stations:

station_id,easting,northing,lat,lon,elev_m,slope_deg,aspect_deg
999901,554229.127877,6249754.167336,56.38929079720315,9.878367933790281,57.95357894897461,37.905677795410156,152.46356201171875
999902,557883.625126,6251682.875208,56.406183385014394,9.93797835390155,47.71885681152344,1.6183518171310425,334.412353515625
999903,557754.877344,6249284.037186,56.384649976373,9.93536386178335,55.55642318725586,1.624045491218567,79.27023315429688
999904,564597.521676,6249506.576003,56.38576377336779,10.046218327254515,62.16939163208008,4.255456924438477,334.4957580566406
999905,551707.752393,6247346.371866,56.36794302138657,9.8370589584469,67.3304443359375,1.3737967014312744,115.48489379882812
999906,551499.63955,6251187.33407,56.40247072835509,9.83444446634275,51.787410736083984,85.89315795898438,105.11801147460938
999907,553826.561777,6256030.973102,56.44572317848233,9.873138949562492,22.089370727539062,3.1103343963623047,1.6468505859375
999908,551633.54944,6261376.38577,56.49398771510929,9.838627653717548,,,
999909,556726.193386,6260986.985054,56.48990379278139,9.921245604401985,,,
999910,567487.021357,6260690.299263,56.485819870454804,10.095893677351198,32.54896545410156,3.8217124938964844,222.7812957763672
999911,562642.918649,6259396.56108,56.47486753329134,10.01693601562108,50.14634323120117,5.993014812469482,43.314117431640625
999912,566698.315706,6257598.602895,56.45816057830541,10.082298318371421,1.242899775505066,2.4549121856689453,29.056732177734375
999913,568328.847189,6254421.393458,56.429387489166395,10.107920341062103,23.738506317138672,2.722937822341919,33.678924560546875
999914,575223.095911,6254910.08896,56.43272888016291,10.219820603374616,15.142535209655762,3.186929702758789,263.4203186035156
999915,568399.788773,6258018.387564,56.46168760213609,10.110011934740477,0.0,0.0,0.0
999916,571885.425034,6259626.087648,56.47561006462536,10.167007862751875,-1.0328165292739868,10.073546409606934,286.3743591308594
999917,574470.977044,6260993.452016,56.48749056594994,10.209362634936879,,,
999918,576575.101763,6251482.95542,56.401728197021164,10.24073654024968,9.114789009094238,1.5538591146469116,261.32293701171875
999919,565367.403446,6251936.194187,56.4074828148552,10.059290787802523,53.06315612792969,2.274454116821289,59.5836181640625
999920,570579.016566,6249064.445447,56.38093731970368,10.142954535337202,43.3277473449707,6.066507816314697,245.16543579101562
999921,574738.02608,6247647.605625,56.367571755714515,10.209885533347508,40.998775482177734,2.257882833480835,332.2108154296875
999922,566007.669174,6246035.824937,56.35439182456608,10.068180060987647,77.24187469482422,13.15386962890625,321.3916931152344
999923,562096.545795,6248291.377861,56.37518270187392,10.005432250347262,77.78323364257812,86.92675018310547,329.6190185546875
999924,562564.218622,6251625.299414,56.40506958801595,10.013798625087073,78.13107299804688,31.939407348632812,12.27001953125
999925,560363.432189,6254321.19025,56.42957312199534,9.97876443080986,9.917754173278809,60.57387924194336,84.58279418945312