Study Area
Borneo island is in the center of maritime southeast Asia, east of Sumatra, north of Java, and west of Sulawesi. It is located at 109° to 119° east longitude and -4° to 7° south latitude (Fig 1). This study divided Borneo island into 8 regions with each region consisting of 9 subregions. Subregions were defined such that the centers were located at latitudes and longitudes of widths 150 pixels (95 km) apart in order to cover all parts of the island. There were 72 subregions, each one comprising 49 pixels in a 7×7 array.
As shown in Figure 1, region A represents Sabah & Brunei, region B represents Sarawak, region C represents North-Kalimantan, region D represents West-Kalimantan, region E represents East-Kalimantan, region F represents West-central Kalimantan, region G represents Central-east Kalimantan region and region H represents South-Kalimantan.
Data
The MODIS Land Surface Temperature (LST) database, available online, was used in this analysis. The website covers all areas of the world and provides LST data during the daytime and nighttime. The data contains average temperatures every 8 days, clear skies permitting, for areas of size 0.859 km2 (ORNL DAAC, 2018; Phan et al. 2008). A sinusoidal projection with tiles of size 10×10 latitude degrees was used to ensure area equality of all pixels, with each tile in turn divided into 1200×1200 pixels. The downloaded LST data was based on the center of the subregion on the island to avoid missing data. Missing values were deleted from the analysis. Natural disasters that can cause unexpected transformations of data behavior were excluded. To maintain comprehensiveness of the LST data, all outliers were kept in the data set. Temperature data, originally stored in the Kelvin scale, was converted to the Celsius scale before analyzing.
Methods
A constant seasonal pattern in LST was assumed to be the same for every year. A cubic spline was used to model the continuous seasonal pattern. The model takes the form:
where S(t) is the Spline function at time t in Julian calendar, and the defined knots are t1 < t2 < ... < tp and (t – tk)+ is (t – tk)>0 for t > x and 0 otherwise. S(t) for t < t1 equals S(t) for t > tp and a, b, ck are the coefficients of the combination between the linear and cubic spline model.
Selecting the position and the number of knots in the spline is important for smoothing. The LST variation in the different regions in the world can be affected by inter-seasonal variation (Singh, Grover, and Zhan 2014). During the rainy season, the LST will be lower in places that have dry and rainy seasons (Jesus and Santana 2017). Variations in LST may be related to heatwaves in tropical areas (in April and May) and rainfall (in June-September) (Gogoi et al. 2019). LST data in the tropical area will use the time during the rainy and dry season to determine the position of the knots (Lukas et al. 2010; Wongsai et al. 2017). Based on the seasonal characteristics of tropical regions, we chose 8 knots, placing 4 at the beginning of the year and the remaining 4 at the end of the year.
A seasonally adjusted time series was used to minimize the seasonal effect on the LST per day for 18 years using a vector of spline fitted values that we estimated from the cubic splines and the average LST per year in a linear model.
A second order autoregressive model AR(2) was used to fit the LST seasonally adjusted. The model is given by:
where Yat is the seasonally adjusted LST at time t, and Yat-1 is the LST at time t-1, t = 1, …, 365 days, α1 and α2 are unknown parameters to be estimated and εt is the random error with zero mean and finite variance (Venables and Ripley 2002).
A multivariate regression model (Mardia et al. 1979) was then used to analyse the seasonally adjusted LST data to detect the spatial correlation. The model is given by:
Y = XB + U [4]
where Y is the outcome matrix of variables with dimension n × m, n is the number of observations, m is the number of subregions, X is a matrix of independent variables of dimension n×q, where q is the number of independent variables, B is a regression parameter matrix with dimension q×m, and U is an unobserved random disturbance matrix.
All analyses and graphical displays were done using R (R Core Team 2018).