Improving agricultural microinsurance by applying universal kriging and generalised additive models for interpolation of mean daily temperature

Agricultural microinsurance has the potential to protect farmers against crop loss caused by extreme adverse weather conditions. Microinsurance policies for smallholder farmers are often designed on the basis of weather indices, whereby weather insurance variables are measured at ground weather stations and then interpolated to the location of the farm. However, a low density of weather stations causes interpolation error, which contributes to basis risk. The objective of this paper is to investigate whether agricultural microinsurance can be improved by reducing interpolation error through advanced interpolation methods, including universal kriging (UK) and generalised additive models (GAM) used with land surface temperature, elevation, and other covariates. Results indicate that for areas with a lower density of weather stations, UK with elevation substantially improves air temperature interpolation accuracy. The approach developed in this paper may help to improve interpolation and could therefore reduce basis risk for agricultural microinsurance in regions with a low density of weather stations, such as in developing countries.


Introduction
Agricultural risk management and insurance are important to the growth of the agricultural sector worldwide because they help producers to stabilise their income and respond to large and unexpected weather risks. A strong agricultural risk management and insurance system in a global context is important for meeting food demand now and in the future and enhancing world food security (FAO 2009). In developing countries farmers often face poverty and constraints such as a lack of understanding of insurance, unaffordable insurance premiums and inadequate access to credit. Agricultural microinsurance has the potential to reduce poverty and stabilise the income of smallholder farmers. Crop microinsurance tends to use an index-based approach, whereby payouts are based on an easily observable index that is correlated to crop yield outcome. Compared to traditional crop insurance, index-based weather insurance approaches have the potential to address some of the difficulties facing smallholder farmers. For example, index-based weather insurance involves relatively low administration costs, and savings on administration could be passed on to farmers through reduced premiums, resulting in more affordable insurance. However, the problem of basis risk, which refers to the mismatch between the indemnity computed from the index and the actual loss experienced on the farm, may partially explain why index-based insurance has not achieved widespread success.
Most of the index-based crop insurance products around the world, including microinsurance products for smallholder farmers, are based on weather indices. Ground weather station measurements (e.g., temperature and rainfall) are the most common source of information for developing index-based weather insurance (Turvey and McLaurin 2012;Mahul and Skees 2007;Makaudze and Miranda 2010). To a lesser extent, satellite-based vegetation or biophysical parameter indices are used to estimate biomass (Baret et al. 2005;Meroni et al. 2014). The underlying hypothesis is that the weather index is highly correlated with crop yield. In the case of ground weather station data, a low density of ground stations can introduce significant basis risk in situations where a farm is located a considerable distance from the next weather station (Collier et al. 2009; Barnett and Mahul 2007). Some alternative indices that overcome the problem of a low density of weather stations are satellitebased vegetation or biophysical indices that do not rely on weather station infrastructure, but consensus has not been reached regarding the extent to which crop yield can be explained by these methods. Other research has examined the combination of vegetation indices, such as the normalised difference vegetation index (NDVI) with ground weather station variables such as temperature (Turvey and McLaurin 2012). To improve the accuracy of agricultural index-based insurance and index-based microinsurance, it is important to help to overcome the issue of basis risk. In the case of a straightforward index-based weather insurance product based on one weather variable such as temperature, relatively simple interpolation methods are often used in practice. A common example of a simple interpolation method is the nearest neighbour (NN) approach, whereby the weather measurement from the station located at the shortest distance from the farm is used in the insurance index. Since measurements are taken at weather stations and not at the insured farm, there is an interpolation error, which contributes to basis risk.
The objective of this paper is to investigate whether agricultural microinsurance can be improved by reducing interpolation error through the use of advanced interpolation methods that integrate satellite information with observations every 1 km, and ground weather station information with observations that vary in geographic distance. The interpolation methods considered include universal kriging (UK) and generalised additive models (GAM) with ancillary information such as satellite land surface temperature, elevation, and other covariates. Performance is evaluated based on each method's ability to develop daily fine-grain 1 km resolution temperature layers over Alberta, Canada.
UK improves air temperature interpolation at a farm location by optimally weighting weather station air temperature information from nearby weather stations and incorporating other variables to capture weather trends. GAMs provide a flexible regression framework to model non-linearity in the relationship between air temperature and covariates by allowing covariate domains to be split into subdomains to allow for flexibility. Air temperature interpolation error, measured by root-mean-square error (RMSE) is tested on two weather station density scenarios designed to replicate expected conditions in both developed and developing countries. Case A has a higher density of weather stations and is intended to represent developed countries with a more developed infrastructure. Case B has a randomly sampled lower density of weather stations and is intended to represent developing countries with less developed infrastructure. Case A is empirically modelled using the existing weather station network in Alberta, Canada. Case B is empirically modelled using a randomly generated subset of the Alberta network with stations over twice the average distance apart.
The covariates are tested for both cases to determine whether they improve interpolation evaluated by out-of-sample RMSE. The UK and GAMs performance are also compared and evaluated. Model performance is compared to the simpler methods commonly used in existing index-based weather insurance programmes as a benchmark; these include NN and inverse distance weighting (IDW) methods.
Open source data was collected for ten growing seasons from 2007 to 2016 including elevation, land surface temperatures, soil type, and land cover. Remote sensing data was gathered from the National Aeronautics and Space Administration's (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS) project for land surface temperatures (MOD11A1), transformed wetness index (TWI) (transformed MCD43A4), and land cover (MCD12C1).
The remainder of this paper is organised as follows. First, a background on the experience of smallholder farmers is given with an emphasis on the importance of agriculture and the potential role of agricultural insurance. Next, the benefits and challenges of implementing index-based weather insurance are outlined. Some common weather insurance interpolation methods used in agricultural insurance are discussed, followed by an overview of some methods that could improve interpolation accuracy. The data section introduces each data source used in the analysis. Subsequently, the methods are explained with a focus on the methodology of kriging, and the results are presented. The final section gives a summary and discusses possible implications.

Background and literature review
There is evidence to suggest that smallholder farming households take significant steps to self-insure their crops and protect themselves from income shocks. These methods of self-insurance can be costly and may not even provide sufficient levels of support during extreme and widespread events such as drought (Dercon et al. 2005;Keyzer et al. 2007;Rosenzweig and Binswanger 1992). Some strategies limit short-term risk but reduce long-term wealth creation. For example, if farmers plant crops with lower but more reliable yields, or if they limit input expenses to reduce the negative shock of crop failure, they are in fact reducing their income. The economic impacts of uninsured weather events may limit households' ability to accumulate wealth (Jensen 2000). Without investment in human capital and wealth accumulation it is difficult to exit the cycle of poverty. Agricultural insurance and other risk management approaches can contribute to improving agricultural productivity by helping producers to invest in more productive agricultural practices (Porth and Tan 2015).
It is well known that agriculture is a risky industry, and farmers may avoid taking risks in the face of extreme loss (Rosenzweig and Wolpin 1993). When farmers limit risk-taking, they often avoid innovations that have the potential to increase productivity. Farmers often rely on their own savings to cushion potential income shortfalls, and investing in technology often requires them to spend money without being certain of a positive return (Dercon 1996). Further, low-income farmers who have borrowed from others may be discouraged from investing funds in higher-reward activities out of fear of the increased default risk; this is described by Boucher et al. (2008) as risk rationing. Indirectly, lenders often charge higher costs for credit to smallholders, which further discourages technology adoption. Also, some banks do not give credit to farmers until they are sure that the conditions are favourable enough to avoid risk to crops already growing, and this creates an obvious shortcoming that is unsuitable for smallholders who need credit to buy inputs at planting time. Lack of access to credit therefore restricts access to agricultural inputs and technologies, such as improved seeds and fertilisers.
Despite the advantages that insurance offers to smallholders, traditional crop insurance has often been difficult to implement outside developed countries. This may be due to moral hazard, whereby farmers who have insured their crop may intentionally reduce the level of crop management or inputs without concern for reduced yield. The problem of adverse selection is also a concern for individual indemnity-based crop insurance, as it may lead to mispricing and the potential for cross-subsidisation of premiums. Further, the high administration costs associated with crop loss adjustment can be difficult to manage, given that developing countries tend to have a large number of small farms.

Indexed-based weather insurance
An alternative to traditional crop insurance is index-based weather insurance, which avoids the critical issues of moral hazard and adverse selection, because the indemnity depends on the performance of an underlying index that is independent of the policyholder's behaviour and cannot be manipulated by the insured customer. Thus index-based insurance is verifiable and transparent. Index-based weather insurance is the most prevalent type, whereby losses are modelled on an underlying index, such as cumulative rainfall or maximum temperature.
In addition to overcoming the challenges associated with moral hazard and adverse selection, index-based insurance may theoretically be more cost efficient, because claims are paid directly based on the index, removing the need for loss assessments at individual farms and saving administration costs. All this makes index-based weather insurance especially useful in the context of many developing countries where crop loss adjustment is difficult and costly, since the majority of farmers are smallholders with less than two acres of land. Although index insurance offers advantages compared to traditional indemnity-based insurance, a new problem of basis risk arises. Basis risk refers to the situation where the loss computed by the index differs from the loss experienced on the farm. This can include the situation in which a farmer experiences a loss but does not receive a payment, or alternatively where a farmer does not experience an actual loss but still receives a payment. Therefore, reducing basis risk is likely to be necessary in order to successfully implement index-based weather insurance (Carter et al. 2014;Clarke 2011;Elabed et al. 2013;Jensen et al. 2018;Lin et al. 2015;Skees 2008).
Generally, basis risk can be categorised as spatial, temporal or variable. Variable basis risk refers to the situation of incorrect variables, or insufficient variables being used to model the relationship with yield. Temporal basis risk refers to the mismatch between the assumptions used to weight the index over the growing period compared to the actual crop growth over different stages of the growing season. The focus of this paper, spatial basis in the context of index-based weather insurance, arises from differences between the weather variables measured at the ground weather station and the actual weather observed at the farm (Norton et al. 2012).
Input interpolation error increases with the distance from the ground weather station to the farm, and this error leads to basis risk. Therefore, the lower density of weather stations in the developing world, together with the problem of incomplete data collection, presents a major challenge to the successful implementation of index-based weather insurance and microinsurance (Hazell et al. 2010). Developing new infrastructure to reduce gaps in the observational network primarily in developing countries is a priority of the World Meteorological Organization (WMO 2016). However, the development of new infrastructure to increase the density of weather stations can be a costly and slow process. Even when the infrastructure has been put in place, long time series of consistent measurements are often required to design insurance products and compute actuarially fair premiums. Within the constraints of existing infrastructure, reducing spatial basis risk can be challenging. One approach to reducing spatial basis risk is through various interpolation techniques, which can be used to predict measurements at unknown locations. In addition, with the increasing availability of free data, a second approach would be to use advanced interpolation techniques that integrate satellite imagery data and other possible variables to improve the weather estimate. With medium-resolution satellite imagery data, weather estimates can be obtained every 1000 m, compared to ground weather stations that are often located more than 20 km apart (Wan et al. 2015).
In commercial applications it is common for input variable interpolation to be computed based on relatively simple methods such as measuring the weather variable at the nearest weather station or weighting weather measurements from nearby stations. For example, in the Corn Heat Unit programme in Alberta, Canada, air temperature (the main input) is measured either from the nearest weather station or from nearby weather stations weighted by the producer (Alberta Financial Services Corporation 2018). These relatively simple methods may have higher interpolation error than more advanced methods, and the additional interpolation error may contribute to spatial basis risk. Using more complex methods and considering additional variables that account for the spatial autocorrelation structure may help improve the accuracy of the index. In developing countries with a lower density of weather stations the interpolation error based on simple approaches is likely to be greater than in countries with a higher density of weather stations. Therefore, more advanced methods, such as UK, may be helpful in reducing basis risk.

Climatology and meteorology methods
Climatology and meteorology literature have proposed more complex interpolation methods to reduce interpolation error. Currently, approaches such as UK, which is the focus of this paper, have not been applied in the context of developing agricultural microinsurance. In climatology, interpolation methods with coordinate information (latitude and longitude) and elevation have been effective in improving air temperature interpolation accuracy (Hijmans et al. 2005).
In addition to elevation, researchers have also found that land surface temperature values (L ST ) are closely correlated to ground air temperature measurements from weather stations (Mildrexler et al. 2011;Vancutsem et al. 2010). Remotely sensed L ST provides daily measurements of each pixel of a 1 km by 1 km global grid if quality issues do not prohibit the creation of the L ST value. Since L ST and weather station air temperature are closely correlated and L ST has a higher resolution compared to the density of ground weather stations, researchers have investigated the use of L ST for improving air temperature interpolation (Mildrexler et al. 2011;Vancutsem et al. 2010;Kilibarda et al. 2014;Parmentier et al. 2015Parmentier et al. , 2014Hengl et al. 2012).
Three common methods for daily air temperature interpolation are kriging, GAMs, and geographically weighted regression (GWR) (Parmentier et al. 2015). For a more thorough review of common interpolation methods used in climatology and meteorology see Hengl (2009). Traditional techniques for daily interpolation of air temperature are often referred to as single timescale methods, which involve using input air temperature data from local weather stations, while multi-timescale methods such as climatology-aided interpolation add temporal climate data to improve prediction accuracy. Both of these techniques can facilitate other covariates such as L ST values to improve accuracy.
In Parmentier et al. (2015) both single timescale and multi-timescale climatology-aided interpolation techniques were tested using kriging, GAM, and GWR methods, and multiple covariates including elevation (Elev), L ST and land cover (L C ) were tested in each model. Employing both-single timescale and multi-timescale methods-they found that adding elevation improved interpolation accuracy and adding L ST improved interpolation accuracy for the summer months. Additionally, variations between the three methods-kriging, GAM, and GWR-were relatively low, although kriging had the lowest RMSE for the multi-timescale approach and GAM had the lowest RMSE for the single timescale approach. However, GAM and GWR may be more effective with a lower density of weather stations. The study area was the state of Oregon in the U.S. Kilibarda et al. (2014) investigated air temperature interpolation globally for 2011 by comparing global RMSE produced by spatio-temporal kriging methods with the inclusion of L ST , elevation information and a topographic wetness layer (TWI). These results were compared to less advanced methods with and without the additional covariates. Spatio-temporal kriging methods are multi-timescale techniques that incorporate temporal information in a kriging framework. Results of this analysis show that global mean, minimum and maximum air temperatures can be predicted with an accuracy between ± 2 °C for areas with a high density of weather stations and from 2 to 4 °C for areas with a low density of weather stations.
Currently, the techniques of climatology and meteorology are not commonly used for index-based weather insurance applications; however, they may significantly reduce basis risk through minimising interpolation error.

Weather station daily temperature
Daily historical mean air temperature measurements from weather stations (T mean ) were collected from the government of Canada's climate and weather database. Measurements were selected for the growing season from 1 April to 1 November for the period 2007-2016 for a total of ten growing seasons. Stations located within Alberta were considered for analysis; air temperature was measured in °C (Government of Canada 2018). T mean measurements were calculated as the average of the minimum and maximum temperatures measured within the daytime interval. To maximise the air temperature sample size, weather stations with continuous daily mean observations for each growing season were selected for the analysis (see Table 1).

Remote sensing land surface temperature images (L ST )
MODIS, operating on the Terra and Aqua satellites, acquires remote sensing information on 36 spectral bands. This information is used to develop many data products including MOD11 land surface temperature (L ST ) and emissivity products. The MOD11 L ST data is retrieved at 1 km resolution by a split-window algorithm and at 6 km resolution by the day/night algorithm. Several L ST products are available in varying resolution and time increments. MOD11A1 daily images were used for analysis (Wan et al. 2015). The MOD11A1 L ST images were converted to °C from Kelvin using the conversion method outlined in the MOD11A1 user manual. Missing values present a challenge for satellite remote sensing data, as cloud cover can obstruct data collection. Inverse distance weighting (IDW) was used to interpolate these missing values. MOD11A1 data was collected for the 2007 to 2016 growing seasons corresponding to the weather station air temperature time series.

Remote sensing transformed wetness index (TWI)
The TWI is a soil moisture index derived from transformations of the MODIS MCD43A4 version 6 Nadir Bidirectional reflectance distribution function/Albedo (BRDF/Albedo) product. The MCD43A4 product is generated daily and available globally. It uses a 16 day retrieval period and a weighting algorithm to generate the daily pixel value at a pixel size of 500 m by 500 m (Schaaf and Wang 2015). The TWI as outlined by Gumbricht (2016) uses two linear transformed components derived from MCD43A4 image bands and creates an index similar to a normalised difference vegetation index (NDVI). First, the components are made from the MCD43A4 bands; then they are used to make the TWI index. The first component (Eq. 1) represents the soil line by arranging dark soil reflectance to light soil reflectance, and the second component (Eq. 2) represents open water. The TWI retains the reflectance value scale factors of the MCD43A4 product.
Once the soil line and water components have been created they can be input as variables in the TWI algorithm. The TWI algorithm is defined by a reference (1) iso-wetness line, and the algorithm is applied using a scale preserving trigonometric rotation. This is combined with a rescaling factor and a calibration factor.
where the reference iso-wetness line is represented as a, the slope of the iso-wetness line as β, the rescaling factor as R, and the calibration factor as C. The values used for the TWI algorithm are: iso-wetness line a (2080), iso-wetness line slope β (− 57), calibration factor C (7000), and rescaling factor R (5942). For the purpose of this analysis, the reference values used in Gumbricht (2016) were used as the default values. The calibration factor C, reference iso-wetness line a and its slope β can be sampled using reference sites.

Remote sensing land classification (L C )
The MODIS MCD12C1 land cover (L C ) product 1 provides land classification maps of three classification schemes: the International Geosphere-Biosphere Programme (IGBP), University of Maryland (UMD), and the Leaf Area Index (LAI). In this analysis we use the IGBP classification scheme, which classifies land cover into 17 categories including cropland, and is of particular interest for agriculture (Friedl and Sulla-Menashe 2015). We consider the effect of cropland on mean air temperatures (T mean ) and discard other land covers because they were not well represented in the sample.

Digital elevation model (Elev)
The Canadian Digital Elevation Model (CDEM) provides an array of vertical elevation values in a common measurement, metres, and vertical accuracy over the Alberta region of approximately 0-5 m (Natural Resources Canada 2016).

EcoDistrict soil type (S T )
Soil type (S T ) is grouped spatially within a national ecological framework. The framework's most specific identifier, the EcoDistrict, is composed of several different attribute classifications including soil type. The attributes are spatially located, and each attribute is derived in unique ways (Agriculture and Agri-Food Canada 2013). The EcoDistrict soil development classification system is based on the Canadian System of Soil Classification. There are a total of 33 classes composed of 31 soil types, a water/rock class, and an urban development class. Soil type was transformed into a binary variable representing Chernozem soil, a fertile agriculture soil, and other. Short descriptions of each covariate are listed in Table 2.

Weather station density: Cases A and B
Weather station data sets were selected to test interpolation improvements in two cases: Case A (full-density) uses Alberta's original weather station network with an average distance between stations of approximately 20 km; Case B (half-density) uses a random subset of 40 weather stations selected from the same network but with the weather stations over twice the distance apart at an average of 40 km. The lower-density network of Case B was selected to simulate weather station density in developing countries (Figs. 1, 2).

Data processing and preparation
Before air temperature prediction can be made and accuracy assessed, each variable must be processed and projected on the same coordinate reference system. Analysis was conducted in the NAD83(CSRS)/Alberta 10-TM(Resource) coordinate reference system. Each MODIS data type comes in science data sets that are extracted using the Geospatial Data Abstraction Library (GDAL) tools (GDAL/OGR contributors 2018). One of the conditions required for air temperature interpolation is that all data are projected in the same coordinate space, otherwise the layers would not match interpolation values and an error would result. MODIS data types use a custom sinusoidal grid projection coordinate reference system, and because the weather station data is unprojected, both data types must be projected to match the Alberta grid. A similar transformation must be done for EcoDistrict derived soil type information. These transformations are made using GIS tools available in R statistical software.
Once the target information has been projected on to the specified coordinate reference system, a data set can be extracted from the desired weather station locations for each growing season day for the period 2007-2016. Each data source presents unique challenges, as different formats are used for MOD11A1 land surface temperature (L ST ), MCD43A4 derived TWI, MCD12Q1 land cover (L C ), soil type (S T ), elevation (Elev), and weather station mean air temperature (T mean ). The MOD11A1 derived data sets are represented by gridded layers with land surface temperature values in each corresponding to a 1 km by 1 km grid and NA values, with pixel values removed due to obstructed measurement or viewing angle error. Similarly, in MCD43A4 the input bands used for the construction of the TWI values are represented by gridded layers but with a different pixel size of 500 m by 500 m. The land cover data set slightly differs from the other MODIS products; MCD12Q1 comes in a gridded data set with 5600 m by 5600 m pixel size and does not have the same issues with missing values, because land cover is composed annually and the method used to derive land cover uses a larger temporal time window, resulting in less missing values. The EcoDistrict soil type values are stored in spatial polygons with the percentage of soil type coverage and soil type identifiers. In total, 33 soil development types are measured in Canada; however, soil types were aggregated to represent fertile agriculture soil type and other soil types. The soil type for each EcoDistrict was selected as the soil type value. Finally, the digital elevation model presents Elev values in a gridded layer with an accuracy of approximately 2 to 6 m. These values are transformed on to the Alberta grid and extracted at the weather station locations.
Before the corresponding data values for each data type can be extracted from the weather station locations, the missing values in the MOD11A1 L ST grid layer and the processed MCD43A4 TWI grid layer are filled by IDW. Next, the value for each variable at the intersection of the weather station points is extracted into a spatial point data set. These spatial points are then used for analysis, and the five methods outlined below are used to predict mean air temperature (T mean ): NN, IDW, OK, UK, and GAMs.

Nearest neighbour (NN)
The simplest approach to index-based weather insurance is to measure the index variables at the weather stations and base farmers' claims on the nearest weather station to their location. This method is easy to use and transparent for both the insurer and the producer. However, it ignores additional available data, and this may lead to a higher interpolation error than when using methods that incorporate additional information.

Inverse distance weighting (IDW)
IDW is similar to the NN approach but uses multiple weather stations instead of only one to make a weighted interpolation of air temperature. Weighting is determined based on distance to nearby weather stations and by an a priori decay rate that is defaulted to p = 2. The predicted air temperature at a given location is then a weighted average of the air temperatures recorded at nearby weather stations.
where λ i is the weight assigned to the air temperature at weather station i, and the sum of the weather station weights must equal one.
Using the IDW method, weights are determined by distance, and all available weather stations are used as input unless a maximum search radius is pre-determined. A maximum search radius was not set for the analysis.
where p is a coefficient used to adjust weights, and d(s 0 ,s i ) is the distance between the target weather station and a known sampled weather station.

Ordinary kriging (OK)
OK was initially developed for mining applications (Krige 1951;Cressie 2015). It is similar to the previous methods in that it uses weather station location information, but it does not use other covariates in the model and does not include information regarding spatial autocorrelation structure to improve interpolation estimates. Air temperature predictions are based on the following model: where µ is the constant global mean and ε is the spatially correlated stochastic error.
OK predictions are made in the same way as IDW predictions, as shown in Formula (2). Weights are determined to reflect the true spatial autocorrelation structure by analysing the semivariances, where semivariance is defined as: where z(s i ) is air temperature at a select weather station, and z(s i +h) is air temperature at a neighbouring weather station at a distance of s i + h.
If there are n unique weather station observations, then n*(n − 1)/2 point pairs are available to compute semivariances. If the semivariances for these point pairs are plotted by distance, we have a variogram cloud that can be binned by distance and further examined. The result of this process is the experimental variogram, which displays the structure of the spatial autocorrelation. Common variogram models include spherical, exponential, and Gaussian models (Cressie 2015;Bivand et al. 2008). Variograms are usually fitted by re-weighted least squares methods, where the weights are determined by the number of point pairs or based on distance (Hengl 2009;Cressie 2015).
Once the variogram is estimated, semivariances are calculated for all the weather station air temperature measurements, and the kriging weights can be solved by: where C is a covariance matrix for n × n, and C 0 is a vector of covariances at the new weather station location.
A disadvantage of OK is that it requires all observations to follow an identically and independently distributed stationary mean. This requirement will likely not hold for the Alberta region, because the geographic area is extensive and temperature variation may be large due to factors such as elevation.

Universal kriging (UK)
UK is an extension of OK that adds explanatory variables to the model to address the stationarity concerns associated with OK. For kriging prediction, variability is separated into global and local components. UK explains the variation in the global component with a trend surface that may include covariates, resulting in residuals that are more likely to follow an individual identically distributed stationary mean. UK requires the fitting of variograms for each day over the time period. The automatic fitting of variograms using the weighted least squares approach was implemented in the automap R package (Hiemstra et al. 2009).

Generalised additive models (GAM)
GAMs are generalised linear models with a linear predictor composed of a sum of smooth functions of covariates (Hastie and Tibshirani 1990). They are frequently used for interpolation purposes and provide flexibility for applying many different relationships. GAMs require the choice of covariate basis functions and knot parameters. A common basis function choice is regression splines, specifically thin-plate splines, which are often used for interpolation purposes (Parmentier et al. 2014;Hijmans et al. 2005). In this paper they are chosen for basis functions. The choice of knot parameters is determined through a cross-validation approach. For more information regarding GAMs, refer to Wood (2006). GAMs are fitted using the mgcv package available in R using the GCV cross-validation method (Wood 2006).

Model assessment
Model interpolation performance is assessed using out-of-sample RMSE, mean absolute error (MAE), and mean error (ME) determined by tenfold cross-validation. The model interpolation accuracy is assessed based on RMSE and MAE measured in °C, and the average bias is assessed by ME. As an example, RMSE is outlined below where ẑ s i is the predicted value at the weather station i, z(s i ) is the observed air temperature value at weather station i, and cross-validation is performed on the data set for each day over the period 2007-2016. Interpolation performance is evaluated for Case A (full-density) and Case B (half-density). The analysis follows three steps for each case. First, model performance is evaluated relative to the simpler interpolation methods, including NN and IDW. In this step each covariate is added in a UK and a GAM framework, and model accuracy is compared relative to the model performance over the aggregated time period. Second, model performance is evaluated for the more complicated models that add additional covariates to the models tested in the first step. These models are evaluated based on their performance relative to the best performing model in the first step. Third, the temporal performance of the models is evaluated. The temporal performance of the simpler interpolation methods, which include the NN and IDW approaches, and the top performing advanced interpolation methods, including UK and GAM, are then compared over the aggregated growing season months April to October.

Full-density network: Case A
First, the UK and GAM interpolation methods are evaluated relative to the simpler methods, which include the NN and IDW approaches. NN and IDW are commonly used interpolation approaches in existing index-based weather insurance programmes, and more advanced methods may improve interpolation accuracy. In Table 3 model performance is displayed for each model aggregated over each growing season for the period 2007-2016 based on standard assessment metrics that include RMSE, MAE and ME. In Table 3 in mod4 to mod8 for the advanced interpolation methods that include UK and GAM one covariate is added to each to identify the covariate that most improves interpolation accuracy. The covariates tested include elevation (Elev), land surface temperature (L ST ), TWI, satellite land cover (L C ), and soil type (S T ).
The results show that the more advanced methods, including UK and GAM, make substantial improvements in model performance over the simpler interpolation methods of NN and IDW in terms of reducing out-of-sample RMSE. The average RMSE for the NN approach over the time period was 1.0063 compared to the top performing single covariate advanced model UK with Elev, which had an RMSE of 0.6915. Further, the UK method seems to outperform the GAM method for each covariate combination.
According to the results, ME seems to be slightly negative for the simpler methods of NN and IDW. GAMs have less negative ME than UK models and far less compared to the simpler methods.
The covariate that seems to have the most impact on interpolation performance over the time period is Elev. The addition of Elev in both the UK and GAM framework results in the top performing models with RMSE of 0.6915 and 0.7253, respectively. This is consistent with interpolation literature, as the addition of elevation information is believed to improve temperature interpolation accuracy. The effect of adding covariate information for interpolation can be observed visually as in Parmentier et al. (2014), who show that adding important covariates improves the fine-grain characteristics of the predicted surface. A surface map for 1 June 2016 is shown in Fig. 3 as an example to illustrate the fine-grain characteristics of each model as well as the improved texture that each covariate adds. The addition of Elev improves the surface visual texture, and additionally the L ST and TWI covariates improve the fine-grain surface. Next, two additional covariate models are assessed to determine whether UK with Elev, and GAMs with Elev can be improved by adding one of the additional covariates, including L ST , TWI, L C , and S T . Table 4 shows the performance results of each of the tested UK and GAM models. The UK models continue to outperform a NN, b IDW, c OK, d UK and elevation, e UK and land surface temperature, f UK and transformed wetness index, g UK and land cover, h UK and soil type. The addition of the elevation covariate in the UK model shown in d produces a fine-grain surface that captures the spatial variability. The addition of other covariates does not seem to capture the spatial variability as well as the addition of elevation. Images have been scaled for visual appearance the GAMs, and the addition of L ST , TWI, and S T to the UK model with Elev seems to improve interpolation accuracy. The top performing models according to RMSE are UK with Elev and TWI (mod3), and UK with Elev and L ST (mod2). The improvement of interpolation performance by using the Elev and L ST covariates is consistent with Parmentier et al. (2015), who show that L ST may improve interpolation accuracy in the summer months. TWI also seems to improve interpolation accuracy over the summer months, and similar information may be gained from the TWI index and the L ST index, but further investigation is needed.
Results suggest that the addition of L ST , TWI, and S T to the UK model with Elev may improve interpolation accuracy over the summer months. Further, examining the results in Table 4, UK and Elev (mod1) has an RMSE of 0.6915 and an MAE of 0.5160, and UK with Elev and L ST (mod2) has an RMSE of 0.6892 and an MAE of 0.5173. These results indicate that the addition of L ST reduces the RMSE but results in an increase in the MAE. A possible explanation is that RMSE penalises large interpolation errors more heavily than MAE, and the addition of L ST may reduce the large interpolation misses, but it may also slightly reduce the model accuracy. Similar results are found when adding the covariates TWI and S T ; however, the addition of the satellite L C covariate does not improve the interpolation accuracy. Further, according to the ME results the addition of a second covariate reduces the negative bias over the time period.
The effect of additional covariates can be shown visually, as the simpler methods result in a smooth surface that lacks spatial variability, while advanced methods often capture much of the spatial variability resulting in a fine-grain textured surface. The addition of the covariates on the fine-grain spatial variability of the interpolated surface is shown in Fig. 4. It shows that the improvement for capturing spatial variability is large with the addition of Elev, but changes are marginal with additional covariates and are not easily recognised through visual assessment. Further inspection of the change in spatial variability is necessary to assess the impact of adding additional covariates such as L ST . Figure 5 shows the change in the standard deviation of interpolation accuracy between a. UK model with Elev, and b. UK model with Elev and L ST . The positive changes shown in Fig. 5c, displayed in light grey, indicate the regions in which the addition of L ST reduces the standard deviation of the interpolated surface.
Conversely, the negative changes in Fig. 5c that are displayed in dark grey indicate areas where the addition of L ST increases the standard deviation of the Fig. 4 Simple methods and multiple covariate models: 1 June 2016 surface maps, Alberta, Canada. Case A. a NN, b IDW, c OK, d UK and elevation, e UK with elevation and land surface temperature, f UK with elevation and TWI, g UK with elevation and land cover, h UK with elevation and soil type. The addition of the elevation covariate in the UK model shown in d produces a fine-grain surface that captures the spatial variability. Adding a second covariate to d UK with elevation such as land surface temperature, TWI, land cover, and soil type does not visually improve the fine-grain texture of the surface. Images have been scaled for visual appearance interpolated surface. According to Fig. 5, the addition of L ST may improve interpolation accuracy over most of Alberta, including the growing regions, and may reduce the accuracy of interpolation in the mountainous regions where very little agricultural activity occurs.
Temporal interpolation accuracy of each method may change throughout the growing season, and these changes may affect weather insurance performance. Further investigation of the temporal changes of interpolation accuracy over the growing season may be useful for improving index-based weather insurance and microinsurance programmes. Table 5 displays the monthly out-of-sample interpolation based on RMSE aggregated from 2007 to 2016. Results are for the simpler interpolation methods, which include NN and IDW, and the more advanced methods based on UK, including OK, as well as the four best performing UK models, including UK with Elev, UK with Elev and L ST , UK with Elev and TWI, and UK with Elev and S T .
According to Table 5, the more advanced interpolation methods Models 3 to Model 7 outperform the simpler methods Models 1 and 2 over each of the aggregated growing season months. Additionally, the UK model with Elev and L ST (Model 5) improves the interpolation performance over UK with Elev (Model 4) for the months May, July, August and September. Also, UK with Elev and TWI (Model 6) improves performance over UK with Elev (Model 4) in June, July, August and September. The addition of covariates appears to reduce some of the large interpolation errors. This is best displayed graphically, as shown in Fig. 6. Negative values represent where model b has higher standard deviation. Model b UK with elevation and land surface temperature has lower standard deviation over most of Alberta, Canada including the major growing areas, and it has higher standard deviation in some of the higher elevation regions where agriculture is less prevalent. Images have been scaled for visual appearance Figure 6a shows NN interpolation RMSE throughout the growing season by month for the period 2007-2016. There are concentrations of large interpolation errors, and Fig. 6b shows that using UK with Elev significantly improves the interpolation accuracy throughout the season; however, there is a large interpolation error in the month of June. By adding L ST this large interpolation error is reduced as seen in Fig. 6c. Controlling for large interpolation errors is important for indexbased weather insurance purposes. For example, a single large interpolation error could result in a policy failure if a farmer suffered a loss but no loss was signalled by the index.

Half-density network: Case B
The half-density network-with twice the average distance between stations as in Case A and intended to simulate developing countries-is examined for each day over the ten growing seasons for the period 2007-2016. The same methods and models are applied to this half-density network as to the full-density network, and model assessment is conducted with the same three-step process described in the section "Model assessment". The overall interpolation accuracy of all models and methods decreases as the distance between weather stations increases, as expected. As seen in Table 6, the NN method's RMSE is measured at 1.1083 for the half-density network, compared to 1.0063 for the full-density network. IDW accuracy decreases but not as significantly as in the NN and the OK approaches. OK has an RMSE change from 0.8092 to 1.1852, which is a fairly large decrease in performance; however, with the addition of covariates using the UK model, performance remains high. Consistent with the full-density network, the top performing single covariate model is UK with Elev (mod4), which has an RMSE of 0.8548. For comparison, the RMSE of the simpler methods on the full-density network is 1.0063 and 0.9613 for NN and IDW, respectively.
According to the ME results, there seems to be a negative bias for most models with the exception of the OK model (mod3), and the UK with Elev model (mod4a). The models with the ME closest to zero are the UK model with L ST (mod5a), GAM model with Elev (mod4b), and the UK model with Elev (mod4a). Table 6 Results: simple methods and single covariate models Case B, n = 40. Despite the lower density of stations in Case B, the top performing model, mod4a, significantly outperforms the simpler models used on the full-density network, Case A. All units are in °C 1 NN is the nearest neighbour method 2 IDW is inverse distance weighting 3 OK is ordinary kriging 4 UK is universal kriging 5 GAM is generalised additive models a NN method uses coordinate information represented by 1, and searches for the nearest station b IDW method uses coordinate information represented by 1, and a decaying search window c mod4a uses Lat, Lon and Elev to interpolate T mean using the UK method, mod4b uses the same variables and the GAM method d mod4a, UK with Elev, is the best performing simple or single covariate interpolation model  The effect of adding the covariates can be visually inspected, which shows that adding important covariates may add a texture element to the fine-grain surface, indicating captured spatial variability. As an example, Fig. 7 shows a surface map of Alberta on 1 June 2016 which displays the surface created for the single covariate kriging models, including OK and the simpler methods NN and IDW, using the half-density network.
In Fig. 7 NN and IDW look similar to the surfaces created with the full-density network; however, the OK model Figs. 7c and 7e-h do not seem to capture the spatial variability. In Figure 7d the UK model with Elev seems to capture much of the spatial variability structure, resulting in a fine-grain surface.
Next, models with two covariates are assessed to determine if UK with Elev, or GAM with Elev can be improved by adding one of the additional covariates, L ST , TWI, L C , and S T . Table 7 shows the performance results of each tested UK and GAM model. The top performing models assessed by RMSE are UK with Elev (mod1a) and UK with Elev and L ST (mod2a). The GAMs do not perform as well as the UK models; however, the model accuracy of the GAMs is closer to the UK methods for the half-density network (Case B) than for the full-density network (Case A).
Results indicate that UK with the Elev covariate is the best performing interpolation method, and the addition of a second covariate does not significantly improve interpolation accuracy. Further inspection of method performance results shows that although UK methods outperform GAMs as they did in the full-density case, the difference in performance is much smaller. In Table 7 the RMSE of UK with Elev (mod1a) is 0.8548, and GAM with Elev (mod1b) is 0.8564. This is compared to the full-density case displayed in Table 6, in which the RMSE of UK with Elev (mod1a) is 0.6915, and GAM with Elev (mod1b) is 0.7253. ME results suggest that the UK models have a small positive bias and the GAMs have a small negative bias-with the exception of the GAMs with Elev and L ST (mod2b). The effect of adding the covariates can be visually inspected. This shows that adding important covariates may add a texture element to the fine-grain surface, indicating captured spatial variability. Figure 8 shows a surface map of Alberta on 1 June 2016 as an example, displaying the surface created for the kriging models with two covariates, including OK and the simpler methods NN and IDW, using the half-density network.
Through visual inspection of Fig. 8, the addition of Elev shown in Fig. 8d improves the fine-grain texture of the surface; however, the impact of the addition of the other covariates on the surface texture displayed in Fig. 8 is not clear.  NN, b IDW, c OK, d UK and elevation, e UK with elevation and land surface temperature, f UK with elevation and TWI, g UK with elevation and land cover, h UK with elevation and soil type. The addition of the elevation covariate in the UK model shown in d produces a fine-grain surface that captures the spatial variability. Adding a second covariate to d UK with elevation such as land surface temperature, transformed wetness index, land cover, and soil type does not visually improve the fine-grain texture of the surface. Images have been scaled for visual appearance Further inspection of the change in spatial texture is necessary to assess the impact of additional covariates such as L ST . Figure 9 shows the change in standard deviation over Alberta for 1 June 2016 as an example, where the positive (light grey) areas show reductions in standard deviation and negative (dark grey) areas show where standard deviation has increased. According to Fig. 9, adding L ST improves the accuracy of interpolation over most of Alberta, including the major agricultural regions, and reduced interpolation accuracy in the high elevation Rocky Mountains.
The temporal interpolation accuracy of the half-density case B may change throughout the growing season, and these changes may have an impact on the performance of weather insurance. Further investigation of the temporal changes of interpolation accuracy over the growing season may be useful for improving indexbased weather insurance and microinsurance programmes. Table 8 displays the monthly out-of-sample interpolation RMSE aggregated for the period 2007-2016 for the simpler interpolation methods including NN and IDW, and the kriging Negative values represent where model b has higher standard deviation. Model b, UK with Elev and L ST has lower standard deviation over most of Alberta, Canada, including the major growing areas, and it has higher standard deviation in some of the higher elevation regions, where agriculture is less prevalent. Images have been scaled for visual appearance  Table 8 shows that the more advanced methods-with the exception of OK (Model 3)-outperform the simpler methods for each month of the growing season. The OK method seems to have difficulty accurately interpolating temperature with the half-density network. Notably for the month of May the OK performance is poor. Further inspection of the interpolation accuracy by month of the more advanced two covariate interpolation models (Models 4 to 7) shows that UK with Elev (Model 4) performs best but not in every month. With the addition of L ST the model UK with Elev and L ST (Model 5) improves interpolation accuracy for the months April, July and August. The addition of TWI, UK with Elev and TWI (Model 6) improves accuracy in the months of July and August. Finally, the addition of S T , UK with Elev and S T (Model 7) improves interpolation accuracy for the month of August. The temporal improvements in some months by adding additional covariates to the UK with Elev model are largely due to the reduction of large interpolation errors throughout the growing season. This effect is best represented graphically and can be seen in Fig. 10. Figure 10a shows NN interpolation RMSE throughout the growing season by month for the period 2007-2016. There are concentrations of large interpolation errors, and Fig. 10b shows that using UK with Elev significantly improves the interpolation accuracy throughout the season; however, there is a large interpolation error in the month of August. By adding L ST this large interpolation error can be controlled, as seen in Fig. 10c. Controlling for large interpolation errors is important for index-based weather insurance purposes.

Summary and implications
Agricultural weather insurance is based on underlying weather variables such as temperature or rainfall. In the case of weather insurance policies for temperature, simple methods are used to interpolate air temperatures observed on the insured farm. Since air temperature measurements are made at weather stations and not at the insured farm itself, there is an interpolation error, which may result in basis risk.
The objective of this paper was to investigate whether agricultural microinsurance could be improved by reducing interpolation error through advanced interpolation methods that integrate satellite information with observations every 1 km, and ground weather station information with observations that vary in geographic distance. The interpolation methods considered include UK and GAM with ancillary information such as satellite land surface temperature, elevation and other covariates. Performance was evaluated based on each method's ability to develop daily fine-grain 1 km resolution temperature layers over Alberta, Canada.
The methods used in this study are UK and GAM interpolation methods. UK improves air temperature interpolation at a farm location by optimally weighting weather station air temperature information from nearby weather stations and incorporating other variables to capture weather trends. GAMs provide a flexible regression framework to model non-linearity in the relationships between air temperature and covariates. Air temperature interpolation error, measured by RMSE, was tested on two different densities of weather stations designed to replicate expected conditions in developed and developing countries. Case A (a full-density network representing developed countries) was empirically modelled using the existing weather station network in Alberta, Canada. Case B (a half-density network representing developing countries) was empirically modelled using a random subset of Alberta's weather station network with the average distance between stations twice that of the full-density network. Data was generated for each day of the analysis.
The covariates were tested for both cases to determine whether they improved interpolation evaluated by out-of-sample RMSE. Also, UK and GAM performances were compared and evaluated. Model performance was compared to the simpler methods commonly used in existing index-based weather insurance programmes that include NN and IDW methods.
Open source data was collected for ten growing seasons from 2007 to 2016, including elevation, land surface temperatures, soil type and land cover. Remote sensing data was gathered from the MODIS project for land surface temperatures (MOD11A1), TWI (transformed MCD43A4), and land cover (MCD12C1).
Results indicate that for both density scenarios the advanced interpolation methods UK and GAM with the addition of covariates substantially improved interpolation accuracy over the simpler methods NN and IDW. For agricultural microinsurance policies in regions where the basis risk caused by interpolation error may otherwise have been too high, these improved interpolation techniques may significantly reduce basis risk.
In Case A (full-density network) UK interpolation with the covariates elevation, and TWI and UK with elevation and land surface temperature resulted in the lowest interpolation RMSE. By using UK methods with covariate information, air temperature interpolation error may be improved for temperature index-based insurance in cases where weather station density is high. Further analysis suggests that the addition of land surface temperature may help to improve interpolation error over the flat growing regions and reduce interpolation accuracy over the high elevation regions where agriculture is less prevalent. UK outperformed the GAM over the full-density network. These advanced methods may be used in air temperature index-based insurance policies to reduce basis risk caused by interpolation error, such as in the Corn Heat Unit insurance programme in Alberta, Canada.
For Case B (half-density network) UK with the covariates elevation resulted in the top performing model over the ten growing seasons. However, the addition of land surface temperatures led to superior air temperature interpolation accuracy in many months and controlled for large interpolation mistakes. The accuracy achieved by using UK with elevation was similar to a much denser weather station network. Despite the low density of the weather station network in Case B, UK with covariates resulted in higher accuracy than the full-density network in Case A using simpler interpolation methods such as NN and IDW. For temperature weather microinsurance policies in developing countries with low-density networks that result in high interpolation error, the basis risk may be significantly reduced by using UK with elevation information and possibly adding land surface temperature information to reduce the large interpolation errors. Additionally, the use of land surface temperatures may improve air temperature interpolation accuracy over growing regions for important growing season months. In general, the UK method showed better performance than the GAM; however, as distance between weather stations increased, the additional performance improvement gained by using UK over GAM decreased and may continue to decrease as the average distance between weather stations increases. This is consistent with Parmentier et al. (2014), who tested the effect of weather station sparsity on GAM and UK methods, and found that GAMs were less sensitive to weather station sparsity. Although findings suggest that UK methods outperform GAMs for index-based weather insurance programmes in regions where weather station density is low, GAMs with covariates should also be considered and empirically tested.
Future research could test whether interpolation error caused by basis risk could be reduced for actual target weather microinsurance areas using the proposed methods. Further, by using actual farm-level data the effect of interpolation error reduction on basis risk could be measured and the effect of basis risk reduction on insurance performance could be simulated. Also, other similar interpolation methods such as GWR or multi-timescale methods such as spatio-temporal kriging could be tested. Additionally, similar studies could investigate the effect of using advanced interpolation methods for reducing basis risk for precipitation-based microinsurance programmes.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 11
Model residuals at weather station locations: 1 June 2016 surface maps, Alberta, Canada. Case A, n = 177. a NN, b IDW, c OK, d UK with elevation (Elev), e UK with Elev and land surface temperature (L ST ), f UK with Elev and TWI, g UK with Elev and land cover (L C ), h UK with Elev and soil type (S T ). Points represent out-of-sample residuals at weather station locations; larger points signify larger residuals. Dark grey points indicate positive residuals and light grey points indicate negative residuals. Visual comparison of model performance is difficult; however, residuals seem to be large in high elevation areas and lower in the growing regions , f UK with Elev and TWI, g UK with Elev and land cover (L C ), h UK with Elev and soil type (S T ). Points represent out-of-sample residuals at weather station locations, larger points signify larger residuals. Dark grey points indicate positive residuals and light grey points indicate negative residuals. Visual comparison of model performance is difficult; however, residuals seem to be large in high elevation areas and lower in the growing regions Mitchell Roznik is a Ph.D. candidate in Actuarial Science and Agricultural Economics at the University of Manitoba, Canada. His research focus is the application of satellite remote sensing data for agriculture insurance. He holds a bachelor's degree in Finance and a master's degree in Agriculture Economics from the University of Manitoba. Milton Boyd is a Professor in the Department of Agribusiness and Agricultural Economics at the University of Manitoba in Winnipeg, Canada. He is also an Adjunct Professor in the area of Risk Management and Actuarial Science at the Asper School of Business of the University of Manitoba. He holds a doctorate degree from Purdue University (USA). Some of his research interests include quantitative finance, risk management, and derivatives.

C. Brock Porth is an Assistant Professor in Biosystems
Katerina Roznik is an ecology and mathematics student at the University of Manitoba, Canada, interested in graduate research related to microinsurance, statistics, and environmental sustainability.