Modiﬁed Linear Scaling and Quantile Mapping Mean Bias Correction of MODIS Land Surface Temperature for Surface Air Temperature Estimation for the Lowland Areas of Peninsular Malaysia

: MODIS land surface temperature data (MODIS T s ) products are quantiﬁed from the earth surface’s reﬂected thermal infrared signal via sensors onboard the Terra and Aqua satellites. MODIS T s products are a great value to many environmental applications but often subject to discrepancies when compared to the air temperature (T a ) data that represent the temperature measured at 2 m above the ground surface. Although they are different in their nature, the relationship between T s and T a has been established by many researchers. Further validation and correction on the relationship between these two has enabled the estimation of T a from MODIS T s products in order to overcome the limitation of T a that can only provide data in a point form with a very limited area coverage. Therefore, this study was conducted with the objective to assess the accuracy of MODIS T s products, i.e., MOD11A1, MOD11A2, MYD11A1, and MYD11A2 against T a and to identify the performance of a modiﬁed Linear Scaling using a constant and monthly correction factor (LS-MBC), and Quantile Mapping Mean Bias Correction (QM-MBC) methods for lowland area of Peninsular Malaysia. Furthermore, the correction factor (CF) values for each MBC were adjusted according to the condition set depending on the different bias levels. Then, the performance of the pre- and post-MBC correction for by stations and regions analysis were evaluated through root mean square error (RMSE), percentage bias (PBIAS), mean absolute error (MAE), and correlation coefﬁcient (r). The region dataset is obtained by stacking the air temperature (T a_r ) and surface temperature (T s_r ) data corresponding to the number of stations within the identiﬁed regions. The assessment of pre-MBC data for both 36 stations and 5 regions demonstrated poor correspondence with high average errors and percentage biases, i.e., RMSE = 3.33–5.42 ◦ C, PBIAS = 1.36–12.07%, MAE = 2.88–4.89 ◦ C, and r = 0.16–0.29. The application of the MBCs has successfully reduced the errors and bias percentages, and slightly increased the r values for all MODIS T s products. All post-MBC depicted good average accuracies (RMSE and MAE < 3 ◦ C and PBIAS between ± 5%) and r between 0.18 and 0.31. In detail, for the station analysis, the LS-MBC using monthly CF recorded better performance than the LS-MBC using constant CF or the QM-MBC. For the regional study, the QM-MBC outperformed the others. This study illustrated that the proposed LS-MBC, in spite of its simplicity, managed to perform well in reducing the error and bias terms of MODIS T s as much as the performance of the more complex QM-MBC method. (DM), and empirical quantile mapping (EQM) onto the Regional Climate Models’ (RCMs) simulated temperature and rainfall data with a reference to the ground meteorological data (1965–2005) for Kaidu River Basin in Xinjiang, China. The ﬁndings indicated that all the bias correction methods performed well and improved the RCMs data with MAE values = 0.76–0.89 ◦ C, PBIAS = 0.20–0.00%, and R 2 = 0.99. A similar research was conducted by Fang et al. [33] in an arid area in northwestern China. Three bias correction methods were applied using daily data from 1975 to 2005 for RCM-simulated temperature data, i.e., LS, VARI, and DM. The monthly scaled bias correction approach for LS, VARI and DM were able to gain PBIAS between 3.04% and 4.78%, R 2 between 0.88 and 0.95, and MAE between 2.35 and 2.52 ◦ C. Another similar research was done by Teutschbein and Seibert [34] in order to remove the bias from RCM-simulated temperature data using four bias correction method, i.e., LS, VARI, DM, and delta change. The evaluation was carried out for ﬁve catchments in Sweden using RCM temperature data from 1961 to 1990. All the bias correction method excepts the delta change resulted in MAE between 0.13 and 1.36 ◦ C. researchers the potential of day T s to estimate T a , weather moderate temperature range between 24 35 C, high high rainfall distribution


Introduction
Air temperature (T a ) is one of the parameters measured at meteorological station. By a standard, T a is measured 2 m above the ground surface to represent the surrounding air temperature. T a has been widely used in various fields of study such as monitoring of temperature trend [1], global warming [2], crop yield modelling [3], and urban planning [4]. Although T a measured at meteorological station is highly accurate and has high temporal frequency, the data provided are dedicated for a point location of the meteorological stations. Hence, T a can only portray the temperature at a local scale and is unable to describe heterogeneous temperature over a large area. This is worsened by the sparse distribution of meteorological stations due to limitations such as topography and operational cost. Since temperature studies for vast and continuous areas are important for many climatic-related applications, numerous studies have been carried out in attempts to solve this limitation, mainly via interpolation of the air temperature data from different localities. However, given the inconsistent distance between weather stations and poorly distributed locations between each retrieval station, the accuracy of the interpolated T a is often compromised [5].
The limitation imposed by the weather stations could be overcome with remotely sensed land surface temperature (LST) data that comes in spatially wide coverages with high temporal resolution. LST measures the skin temperature that is based on the radiations reflected from the earth surface and detected by satellite sensors. The radiation observed are largely controlled by the type of land surface and atmospheric condition such as air temperature, solar radiation, and cloud condition [6,7]. The launch of Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on board of Terra and Aqua satellite has enabled remotely sensed day and night-time skin temperature data to be retrieved daily and every 8-days at 1 km resolution. The Terra satellite provides skin temperature data known as MOD11A1 (daily) and MOD11A2 (every 8-days) at acquisition times approximately 10:30 a.m. and 10:30 p.m. local time, while the Aqua satellite measures skin temperature data known as MYD11A1 (daily) and MYD11A2 (every 8-days) at acquisition times approximately 1:30 a.m. and 1:30 p.m. local time. These skin temperature data are calculated from band 31 (10.78-11.28 µm) and 32 (11.77-12.27 µm), which are designated for skin temperature.
Since these MODIS LST (T s ) products can provide temperature data for mixed elements of targeted earth's land surface [8], they have been a great input in many environmental applications such as hydrology [9], agriculture [10], ecology [11], drought [12], evapotranspiration [13], and soil moisture estimation [14]. Although T s can provide high spatial and temporal temperature data, it suffers from data uncertainty and should be assessed for accuracy; the quality of the day T s is highly dependent on the radiant temperature of the land surface elements, seasonality [15], and the atmospheric condition of particular days [16]. Day T s is more complex than the night T s because the former is acquired from a mix 1 km 2 targeted earth surface that are much affected by solar radiation and other factors [17,18]. Additionally, Zhang et al. [19] revealed that presence of clouds or mixed clouds-earth surface tended to distort minimum and maximum T a estimations.
Prior to utilization of T s to estimate T a values, researchers such as Jin and Dickinson [6], Simó et al. [20], and Sobrino et al. [21] suggested that the complex relationship between both temperatures need to be understood. The vital difference between the T s and T a is that both are measured at certain different heights from the earth surface. While T a is recorded at 2 m above an open cleared land, day T s is acquired from a mix 1 km 2 targeted earth surface, and prone to wind and solar radiation interference. Furthermore, estimation of T a from T s data is also affected by the changes in elevation. Phan et al. [22] revealed that with every increase in elevation by 1000 m, T s value would drop by 3.8 to 6.1 • C for night T s and 1.5 to 5.8 • C for day T s .
Many studies had indicated the performance of MODIS T s to estimate T a in various regions of the world. For instance, Lu et al. [23] conducted a comparison between day-time daily T s with in-situ data for arid area in Northwest China. The finding disclosed that the relationship between both data varied accordingly to seasonality and weather conditions with RMSE between 2.39 and 3.05 • C. Furthermore, Marques da Silva et al. [24] evaluated the performance of accumulated monthly MODIS-derived T s against the accumulated monthly T a , in continental and near sea regions of Portugal. It was discovered that both data were highly linearly related to each other with R 2 > 0.98. In addition, Zhang et al. [25] reported that the MAE between MODIS-derived T s to maximum T a was between 2.8 and 4.1 • C, whereby these values were highly affected given the presence of partial cloud covered pixels. Moreover, El Kenawy et al. [15] in a study to compare day T s against maximum and minimum T a in Egypt reported an overestimation with an average of 5 • C according to different seasons. The difference was the highest during the summer and lowest during the winter, while comparison between the night T s against the min T a demonstrated low over-/underestimations with an average lower than 1.5 • C according to different seasons.
Apart from evaluating the T s products accuracy, many researchers have dedicated their efforts into applying various T s to T a correction methods. Benali et al. [26] applied statistical models and MODIS T s to accurately estimate maximum, minimum, and average T a data of 106 meteorological stations in Portugal for 10 years (2000-2009) period. The study revealed that T a could be accurately estimated from MODIS T s with RMSE between 1.83 and 1.74 • C. Furthermore, Huang et al. [10] estimated and mapped daily mean T a of 23 meteorological stations from T s using a linear calibration method for data range from 2003 to 2011. The proposed model was able to estimate T a with RMSE and MAE of 2.41 and 1.84 • C, respectively. In another study, Zhu et al. [27] incorporated the temperaturevegetation index method into the estimation of T a from day T s over Xiangride river basin at north Tibetan Plateau. The method was able to improve the estimation of maximum T a from day T s with RMSE and MAE of 3.79 and 3.03 • C from 7.45 and 6.21 • C, respectively. Williamson et al. [28] developed an interpolated curve mean daily surface temperature (ICM) method that interpolates MOD11A1 data values by utilizing ground air temperature curves in southwest Yukon, Canada. Validation for the ICM method recorded R 2 and RMSE between 0.72-0.85 and 4.09-4.90 • C, respectively. Further south, Meyer et al. [29] conducted a study to map daily air temperature from spatial MODIS LST data using linear regression and machine learning approaches for Antarctica. The air temperatures from 32 weather stations and the MODIS products, i.e., MOD11A1 and MYD11A1 were used as the absolute and simulated data, respectively. The machine learning method (R 2 = 0.71 and RMSE = 10.51 • C) slightly outperformed the simple linear regression (R 2 = 0.64 and RMSE = 11.02 • C).
The T s to T a correction methods for regional study were also actively carried out by a number of researchers. The analysis at a regional scale has enabled a validated map of estimated T a to be produced. Recently, Hereher and El Kenawy [30] made use of MYD11A2 to produce regional daytime, monthly, and annual surface air temperature map for Egypt using regression correlation and statistical significance analysis. Agreement between both daytime and monthly daytime LST data to air temperature data are R 2 = 0.87-0.95 and 0.77, respectively. It is reported that the produced gridded map with spatial resolution of 1 km (MODIS LST) can act as an effective tool to delineate a temperature variation throughout a region. Meanwhile, Rhee and Im [31] proposed a diurnal air temperature change model for air temperature estimation in South Korea regions with limited ground data. The model considered the sunrise, sunset, and solar noon times as the input parameters. Ground air temperature data retrieved from 60 weather stations and two MODIS LST data were used (MYD11A1 and MYD11A2) as the absolute and the simulated data, respectively. The model produced a monthly corrected air temperature estimation with MAE values between 1.73 and 1.86 • C.
Despite all the mentioned T s to T a correction methods above, mean bias correction method is one of the promising methods that has been widely used to correct Regional Climate Models' (RCMs) simulated temperature data but are still not widely explored for MODIS LST application. Recently, Luo et al. [32] applied five bias correction methods, i.e., linear scaling (LS), daily translation (DT), variance scaling (VARI), distribution mapping (DM), and empirical quantile mapping (EQM) onto the Regional Climate Models' (RCMs) simulated temperature and rainfall data with a reference to the ground meteorological data  for Kaidu River Basin in Xinjiang, China. The findings indicated that all the bias correction methods performed well and improved the RCMs data with MAE values = 0.76-0.89 • C, PBIAS = 0.20-0.00%, and R 2 = 0.99. A similar research was conducted by Fang et al. [33] in an arid area in northwestern China. Three bias correction methods were applied using daily data from 1975 to 2005 for RCM-simulated temperature data, i.e., LS, VARI, and DM. The monthly scaled bias correction approach for LS, VARI and DM were able to gain PBIAS between 3.04% and 4.78%, R 2 between 0.88 and 0.95, and MAE between 2.35 and 2.52 • C. Another similar research was done by Teutschbein and Seibert [34] in order to remove the bias from RCM-simulated temperature data using four bias correction method, i.e., LS, VARI, DM, and delta change. The evaluation was carried out for five catchments in Sweden using RCM temperature data from 1961 to 1990. All the bias correction method excepts the delta change resulted in MAE between 0.13 and 1.36 While many researchers have demonstrated the potential of day T s to estimate T a , most of the studies were conducted in temperate, non-equatorial regions. In light of this, we attempted to (1) evaluate the performance of four types of MODIS T s products, namely, MOD11A1, MOD11A2, MYD11A1, and MYD11A2 for lowland area along east and west coast of Peninsular Malaysia. Our study is the first to conduct such a test for equatorial regions, which extend from 0 • to 10 • in the North as well as the South, and whose climate is characterized by hot and humid weather with moderate temperature range between 24 and 35 • C, high humidity (>80%), high rainfall distribution depending on the monsoonal seasons and high cloud cover throughout the year. Additionally, (2) we comparatively assessed the performance of two mean bias corrections (MBC) approaches, i.e., linear scaling (LS) and quantile mapping (QM), in improving the estimation of T a from the four T s products.

Study Area
Peninsular Malaysia or also known as West Malaysia is located at the south end of Asia continental plate ( Figure 1). It is geographically located between 1 • N and 7 • N and 99 • E to 105 • E. It is characterized by hot and humid tropics with temperature ranging between 23 and 35 • C, daily humidity levels exceeding 80%, and high rainfall and cloud cover throughout the year [35]. The weather characteristic is dependent on two monsoonal wind seasons, i.e., Southwest and Northeast monsoons and elevation. The Southwest monsoon outsets from June to September while the Northeast outsets from November to March, with the latter brings higher rainfall density than the former [35]. The coastal region is characterized by lowland regions, and there is an increase in elevation towards the center of peninsular, where mountain ranges such as the Titiwangsa, Benom and Tahan mountain ranges act as the backbone and separates between the eastern and western part. This geographic condition results in lower temperature towards the center of the peninsular compared to the coastal areas.

Ground Measurement
The ground measurements of air temperature (T a ) were collected from 36 automatic weather stations (AWS) belonging to Department of Environment (DOE) for year 2003 until 2016. The AWS provided three type of daily T a , known as maximum, minimum, and average T a . The maximum T a is the maximum temperature and highly associated with the highest solar radiation during the day approximately at 11:00 a.m. until 2:00 p.m. [36] while the minimum T a is the minimum temperature during the night of a particular day and the average T a is the average temperature for the whole particular day. Nonetheless, for this study, only maximum T a were utilized for further analysis due to its compatibility with day T s which are acquired at approximately 10:30 a.m. and 1:30 p.m. local time.
The AWSs are measured 2 m above open cleared ground and commonly located at the government facilities such as hospitals, government offices, and schools. Additionally, the AWSs location are located away from human activities on the land, i.e., agricultural activity, that may change the earth's surface in order to ensure there are no drastic changes in the T a readings. Figure 1 illustrate the location of AWSs in Peninsular Malaysia, while Table 1 tabulates the information in greater details. The AWSs are mainly distributed around the east and west coastal region of the peninsular, while the digital elevation model (DEM) highlights that the AWS are located at lowland region with elevation below than 60 m.

MODIS LST Product
Four MODIS LST products were used in this study: MOD11A1, MOD11A2, MYD11A1, and MYD11A2 with 1 km spatial resolution in a 1200 by 1200 km grid that were obtained from https://search.earthdata.nasa.gov/search (accessed on 2 January 2020). A total of 2 tiles data was used to cover Peninsular Malaysia with horizontal line 27 to 28 and vertical line 08 (H27V08 and H28V08). The data were downloaded from January 2003 until December 2016. The MOD11A1 and MOD11A2 were retrieved from the MODIS sensors onboard Terra satellite with the acquisition day-time at approximately 10:30 a.m. and night-time 10:30 p.m. for daily and 8-days composite data, respectively. Additionally, the MYD11A and MYD11A2 were retrieved from the MODIS sensors onboard Aqua satellite with the acquisition day-time at approximately 1:30 a.m. and 1:30 p.m. for daily and 8-days composite data, respectively. The AWSs data were first pre-processed for outlier identification and removal using the boxplot outlier analysis [37], following Equations (1) and (2).
where Q1 and Q3 is the first and third quartile of the data. The extreme values exceeding these quartiles were masked out. Then, the 8-days average T a data were calculated by averaging each 8-days of T a data from 2003 until 2016. A total of 5113 daily T a and 630 8-days T a rows data were prepared for pairing with daily T s and 8-days T s , respectively. Figure 2 illustrates the mean T a distribution from 2003 to 2016. The east coast region possessed lower mean T a values than the west coast region, of which the center to northern part showed higher mean T a than its southern counterpart.

Pre-Processing of T s
The downloaded tiles were projected to WGS 1984 UTM Zone 47N and 48N. The original temperature in Kelvin (K) were then converted to degree Celsius ( • C) according to Equation (3) [38].
where 0.02 is the scale factor, T is the absolute temperature in Kelvin and T s is the land surface temperature in • C. Cloud-contaminated pixels were identified using boxplot outlier analysis as in Equations (1) and (2). On the average, high numbers of the cloudcontaminated pixels showed unreasonably low temperature value that was less than 15 and 10 • C for lowland and highland region, respectively. The pixels with extreme values were removed and consequently, resulted in the missing data values. Further, the T s at the same location with the AWSs was extracted and the available pixels of T s and T a of the same locality were paired. Figure 3 illustrates the map of T s distribution in Peninsular Malaysia for clear days on 30 March 2016, 9 March 2014, and 14 June 2012 whereby the T s were able to depict the heterogeneity of surface temperature for the Peninsular Malaysia, despite presence of white patches indicating the missing data; mostly because of cloud cover. The images reflect that T s are highly associated with elevation: towards the center region of Peninsular Malaysia where the mountain ranges are located at, and the temperatures are lower compared to the coastal region. Besides, the east coast region reflects lower temperatures range than the west coast, because of the Northeast monsoon that outsets from November until March each year bringing heavy rains on the east coast region and results in lower surface temperature. In general, the MOD11A1 product depicts a lower temperatures range than the MYD11A1 due to the acquisition time of MODIS sensors onboard of Terra and Aqua.  Figure 4. It can be noticed, MYD possess a higher surface temperature range (from 30.86 to 41.65 • C) than MOD (from 29.13 to 38.10 • C). Besides, the west coastal area demonstrates a higher temperature range than the east coastal.

Region Delineation
To portray the heterogenous nature of T a , the estimation of T a from T s at a regional scale is necessary. Besides the analysis for each station, the analysis is also carried out for specific climatic regions outlined following Suhaila and Yusop [39]. The determination of the boundary for each region is based from the temperature and rainfall trends and the geographic condition of Peninsular Malaysia. The west and east regions depict different temperature trends as the east regions revealed lower mean temperature values than the west region for both T a and T s (Figures 2-4). The Titiwangsa Range that runs through the central Peninsula acts as a barrier to hinder the heavy rainfall and cooler weather from reaching the west region during the Northeast Monsoon [39]. In addition, for this study, the west region [39] is further divided into central and northcentral. This is because the Klang Valley areas (Kuala Lumpur, Putrajaya, and Central of Selangor) revealed higher mean temperature values than the northcentral region (North of Selangor) [40]. Besides, the central region is associated with urban and industrial area and exhibits the effect of urban heat island (UHI) [1,39,40]. Figure 5 illustrates the lowland areas along east and west coast of Peninsular Malaysia (DEM = 1-60 m) that are applicable for T s to T a estimation and the outline for the five regions, i.e., northwest, northcentral, central, southwest, and east regions.

Preparing Regional Data
The pre-processed T a and T s from 36 stations were stack as regional data according to the information in Figure 5 and Table 2. There were 5 sets of regional datasets, and the T a and T s were labelled as T a_r and T s_r , respectively. Four sets of T s /T s_r data from MOD/MYD11A1 and MOD/MYD11A2 with the same locality as the AWSs were paired with their daily and 8-days T a /T s_r . Then, the evaluation was carried out in leave-two-years-out cross validation manner for both by station and by region analysis [41]. The leave-two-years-out validation indicates that for each station/region, 2 years data from the collocated data were excluded from the training set and was used as the validation set. The evaluation was repeated for 13 times with different sets of leave-two-years-out data, and the results were averaged to gain the final values for the evaluation metrics.
The assessment of the T s /T s_r against T a /T a_r was carried out for both pre-and post-bias corrected data using four evaluation metrics: root mean square error (RMSE), percentage bias (PBIAS), mean absolute error (MAE), and correlation coefficient (r) [10,26]. The equation for RMSE, PBIAS, MAE, and r as shown as Equations (4)-(7), respectively.
where T si = land surface temperature data, T ai = air temperature data, and N = sample size. Root mean square error (RMSE) is a measure to determine error between two variables comparing between simulated and observed value. A low RMSE indicates a good fit between two variables.
Percentage bias (PBIAS) is a measurement of the average tendency of the simulated values to be larger or smaller than the observed one. PBIAS of 0.0 indicates optimal value while positive values indicate overestimation and negative values indicate underestimation bias.
Mean absolute error (MAE) determines the average absolute difference between two variables. Likewise, RMSE, a low MAE indicates a good fit between two variables.
where T s is the average land surface temperature data and T a is the average air temperature data. Correlation coefficient, r, measures strength and direction of relationship between both variables. The r value of 1 indicates a very strong positive correlation, 0 indicates no correlation, and −1 indicates strong negative correlation.

Mean Bias Correction (MBC)
There are many types of MBC method that are commonly used in climate study to reduce the error value in estimated data; they vary from the simple methods such as linear scaling, daily translation, and variance scaling to the more complicated approaches such as distribution mapping and empirical quantile mapping [32]. In this study, linear scaling (LS-MBC) and quantile mapping (QM-MBC) techniques are used for the estimation of T a from T s . Both techniques can yield high bias correction and suitable for climate study in region with less extreme temperature range [32]. Moreover, in our study, the LS technique is modified by manipulating the CF values according to a number of set conditions in the equation.
The LS-MBC technique utilizes a single correction factor retrieved from the subtraction result between the mean observation and simulated data over a period of time [32]. For temperature, the simulated data are corrected by an additive function of CF. The LS-MBC can perform the correction either using monthly or daily mean value throughout the year [34]. Nonetheless, the former is more suitable for extreme weather region with extreme daily temperature throughout the year such as in region with four seasons, while the latter is appropriate for region with moderate temperature range such as Peninsular Malaysia. Despite that, this study will examine the correction factor value from both monthly and constant daily average. The LS-MBC are shown in Equation (8): where T sc,d is known as the corrected T s,d for the particular d day that has undergone additional term from subtraction of µT a = mean T a to µT s = mean T s . The mean values are the average obtained from the whole tested days and the average monthly data throughout 2003 until 2016. For the correction, each AWS station will have one constant daily CF and 12 monthly CF values as in Equation (9) and the calculation can be simplified to Equation (10), T sc,d = T s,d + CF.
Since each station will undergo an addition of CF values as written in the LS-MBC formulation, Equation (10) will result in an over correction of bias values especially for the T s,d that already approximate the T a range. Therefore, a modification of the LS-MBC formula is needed to change the bias correction method into a more robust formula. This problem could be overcome by identifying the product of T s,d − µT a and further manipulation of the CF range value depending on the product of T s,d − µT a .
The CF values are manipulated into 6 conditions, based on the positive or negative value of the particular CF. A larger µT a than the µT s will result in positive CF. When the CF is positive, the calculation will be as follows (Equation (11)): The followings are the calculation for the modified MBC, when the µT a is smaller than the µT s and results in negative CF values (Equation (12)): Depending on the subtraction result from T s,d − µT a as the condition, an adjustment of CF values will take place. Each T s,d − µT a will only satisfy one condition.
Meanwhile, quantile mapping mean bias correction (QM-MBC) works by identifying the empirical cumulative density functions (ecdf) of T a as the absolute data and T s as the simulated data. The generated correction function depends on the quantile values between 0% and 100%. Simultaneously, the correction function is applied to unbiasedness the T s data according to the temperature pattern distribution [42]. The calculation can be expressed as Equation (13) [32]: where ecd f −1 denotes the inverse ecd f . Table 3 tabulates the maximum, minimum, and average value of RMSE, MAE, PBIAS, and r for pre-MBC (T s ), post-MBC linear scaling using constant daily CF (T scd ), post-MBC linear scaling using monthly CF (T scm ) and post-MBC quantile mapping (T scq ) against T a according to the leave-two-years-out evaluation for 36 AWS station. For the pre-MBC (T s ), the average RMSE values were much higher for the MYD products (MYD11A1/11A2), i.e., 5.13 and 5.09 • C compared to the MOD products (MOD11A1/MOD11A2), i.e., 3.42 and 3.48 • C. The error value for MYD products was approximately 1.5 times higher than the MOD products. Generally, the same pattern was observed for the maximum and minimum RMSE values.
A greater detail as illustrated in Figure 6 revealed that for pre-MBC (T s ), MOD11A1/11A2 had lower RMSE values ranging from 1.85 to 5.83 • C with only three of the stations showing RMSE > 5 • C: CA0002 for MOD11A1 and CA0019 and CA0058 for MOD11A2. Meanwhile, MYD11A1/11A2 depicted higher RMSE values ranging from 2.37 to 8.28 • C with 54.67% of the stations showing RMSE > 5 • C. MYD11A1/MYD11A2 data exhibited 21 and 18 stations with RMSE > 5 • C. Generally, the distribution of stations with higher RMSE values was mainly in the west coast region rather than the east coast region. Table 3. RMSE, MAE, PBIAS, and r for pre-MBC (T s ), post-MBC linear scaling using constant daily CF (T scd ), post-MBC linear scaling using monthly CF (T scm ), and post-MBC quantile mapping (T scq ) against T a according to the leave-two-yearsout evaluation for 36 AWS station. Nonetheless, the RMSE distributions for post-MBC products (T scd , T scm , and T scq ) showed that the majority of the RMSE values were between 1 and 3 • C. Specifically, for T scd and T scm , only 9.72%/12.50% and 6.94%/9.72 stations showed RMSE value > 3 • C for MOD11A1/11A2 and MYD11A1/11A2, respectively. Meanwhile, for T scq , 16.67% stations showed RMSE values > 3 • C for both MOD11A1/2 and MYD11A1/2. Lower numbers of stations with RMSE values > 3 • C were recorded for the 8-days MOD/MYD11A2 (T scd = 1 and 3 stations, T scm = 1 and 2 stations, and T scq = 3 stations for both) than the daily MOD/MYD11A1 (T scd = 6 stations for both, T scm = 4 and 5 stations, and T scq = 9 stations for both). Table 3, the pre-MBC T s products tended to overestimate the T a . The average PBIAS of MOD11A1 and MOD11A2 depicted lower values (1.36% and 4.40%) than the MYD11A1 and MYD11A2 (8.88% and 10.74%). The same trend was observed for the maximum PBIAS values as MYD11A1 and MYD11A2 manifested high maximum PBIAS of 22.02% and 25.56%, respectively, while MOD11A1 and MOD11A2 showed low maximum PBIAS of 9.90% and 14.22%, respectively. In contrast, for the minimum PBIAS values, MOD11A1 and MOD11A2 recorded −14.27% and −11.70%, suggesting higher magnitudes of underestimation in comparison to the MYD11A1 and MYD11A2 that had −11.77% and −6.74%. It can be noted that MOD11A1 demonstrated the highest minimum and lowest maximum and average PBIAS, while MYD11A2 depicted the lowest minimum and highest maximum and average PBIAS. Meanwhile, the average PBIAS values for the post-MBC (T scd , T scm , and T scq ) exhibited reduced over-/underestimations with average PBIAS values between −1.85% and 2.46%, −1.47% and 2.55%, and −3.26% and −0.37%, respectively. In particular, the daily MOD11A1 depicted negative average PBIAS values (T scd = −1.85%, T scm = −1.47%, and T scq = −3.20%). Moreover, all post-MBC T scq products demonstrated negative average PBIAS values with daily MOD/MYD11A1 having higher magnitude of negative bias values (−3.20% and −3.26%) than the 8-days MOD/MYD11A2 (−0.37% and −0.42%), respectively. Meanwhile, the maximum PBIAS values for T scd and T scm exhibited lower magnitudes of positive PBIAS values for the daily MOD/MYD11A1 (T scd = 2.43% and 3.35% and T scm = 2.10% and 3.15%) than the 8-days MOD/MYD11A2 (T scd = 3.76% and 4.78% and T scm = 3.82% and 4.73%). The maximum PBIAS values for T scq were low negative for daily MOD/MYD11A1 (−0.62% and −1.07%) and low positive for 8-days MOD/MYD11A2 (1.20% and 1.02%). Moreover, it is noted that the minimum PBIAS values for post-MBC (T scd , T scm , and T scq ) for the daily MOD/MYD11A1 depicted higher magnitude of negative PBIAS values (T scd = −6.33% and −6.34%, T scm = −5.63% and −5.51%, and T scq = −6.11% and −5.51%) than the 8-days MOD/MYD11A2 (T scd = −3.55% and −3.37%, T scm = −2.85% and −2.71%, and T scq = −2.11% and −1.65%)

Assessing the average PBIAS values in
As illustrated in Figure 7, the pre-MBC T s demonstrated that overestimation of T a was more frequent than the underestimation, with 69.     Table 4 tabulates the maximum, minimum, and average value of RMSE, MAE, PBIAS, and r from pre-MBC (T s_r ), post-MBC linear scaling using constant daily CF (T scd_r ), post-MBC linear scaling using monthly CF (T scm_r ), and post-MBC quantile mapping (T scq_r ) against T a_r according to the leave-two-years-out evaluation for 5 regions. The tabulated data indicated that the average RMSE values for the pre-MBC T s_r of the by region comparison reflect same pattern as the by station comparison. The MYD products (5.23 and 5.42 • C) indicate 1.5 times higher error than the MOD (3.33 and 3.55 • C). The same pattern is shown by maximum and minimum RMSE values. Meanwhile, the post-MBC products showed an overall low average RMSE < 3 • C. Specifically, the RMSE range values were recorded between 1.96 and 3.71 • C, 1.95 and 3.51 • C, and 1.88 and 3.61 • C for T scd_r , T scm_r , and T scq_r , respectively. Figure 10 illustrates the RMSE for pre-and post-MBC products. It can be observed that the error values were well reduced for all regions and products.  Table 4. RMSE, MAE, PBIAS, and r from by region comparison for pre-MBC (T s_r ), post-MBC linear scaling using constant daily CF (T scd_r ), post-MBC linear scaling using monthly CF (T scm_r ), and post-MBC quantile mapping (T scq_r ) against T a_r according to the leave-two-years-out evaluation.  10. Comparison of evaluation metrics score between T s_r , T scd_r , T scm_r , and T scq_r against T a_r .

MAE
A similar pattern was observed between pre-MBC T s_r MAE and RMSE values as the MYD products result in 1.5 times higher error values than the MOD ( Table 4). The average values were recorded at 2.79/3.09 • C and 4.68/4.89 • C for MOD11A1/11A2 and MYD11A1/11A2 products, respectively. The maximum and minimum MAE values revealed two different range of error values between the MOD and MYD products, i.e., 2.36 to 5.75 • C and 4.36 to 5.75 • C, respectively. The post-MBC products depicted an overall low average MAE value < 2.5 • C (T scd_r = 2.12 to 2.27 • C, T scm_r = 2.10 to 2.27 • C, and T scq_r = 1.78 to 2.08 • C). Figure 10 showed an apparent reduction in MAE values between the pre-and post-MBC products except for MOD11A1.

Correlation Coefficient (r)
The r values observed from Table 4 illustrated low r values for all pre-MBC products with average r between 0.16 and 0.29. In detail, the daily products tabulated slightly higher r range between 0.23 and 0.29 than the 8-days products with r range between 0.16 and 0.17. The post-MBC products showed a slight improvement in the r values with T scq_r revealing the highest r range for average and maximum values (average = 0.25 to 0.29 and maximum = 0.40 to 0.52), followed by T scm_r (average = 0.24 to 0.29 and maximum = 0.27 to 0.39) and T scd_r (average = 0.19 to 0.26 and maximum = 0.25 to 0.31). Figure 10 shows that the performance of post-MBC varies according to the region. Notably, the T scq_r performed best for central and northcentral region. Other than that, an overall slightly improved r values can be observed throughout the MBC products.

By Station Performance
The assessment of MODIS T s against T a (pre-MBC) indicated high variabilities in the relationship between both data. For instance, MOD11A1 was found to have the closest relationship to T a data with the lowest average RMSE of 3.42 • C and PBIAS of 1.36%, followed by MOD11A2 with RMSE of 3.48 • C and PBIAS of 4.40%, MYD11A1 with RMSE of 5.13 • C and PBIAS of 8.88% and finally, MYD11A2 with average RMSE of 5.09 • C and highest PBIAS of 10.74%. These metrics suggested that MOD11A1/11A2 was better at estimating the T a than MYD11A1/11A2. This could be attributed to the time of data acquisition. The former is acquired at 10:30 a.m. local time, while the latter is acquired at 1:30 p.m. local time whereby later in the noon, solar activity for the tropical region such as Peninsular Malaysia is the highest. Consequently, the MOD data are a better proxy for maximum daily T a given that the MYD T s tended to overestimate the T a .
Additionally, the daily T s , i.e., MOD/MYD11A1 and 8-days T s , i.e., MOD/MYD11A2 depicted different characteristics of performance as the former illustrated lower errors and lower under-/overestimations than the latter. Daily data are known to be better at estimating T a than 8-days data due to the inconsistency of the numbers of days used to count the 8-days T s [38]. Although the 8-days T a are averaged from the daily T a for a complete duration of 8-days, 8-days T s are subject to missing data resulted from cloud coverage and hence, could be averaged from daily data as low as from 2 to 8 days. This inconsistency is believed to contribute to low accuracies in estimation between both datasets.
Furthermore, although researchers such as El Kenawy et al. [15] recommended that corrections of T s against T a were needed if only the score of PBIAS > ± 5%, the usage of MOD11A1 (average PBIAS = 1.36%) for lowland area in Peninsular Malaysia region without any correction to represent the T a can still be considered inappropriate given its average error (3.42 • C) is too large to be considered for many applications [43].
Post-MBC (T scd , T scm , and T scq ) data assessment against T  [43], estimated products at RMSE approximately < 3 • C are considered acceptable, and most of the recent studies for T a estimation reported RMSE values < 3 • C.
Although the modified MBC demonstrated a good performance at reducing PBIAS over-/underestimation range, the underestimation in the MOD11A1 data was significantly high in post-MBC, i.e., T scd = 86.11%, T scm = 77.78%, and T scq = 100% versus pre-MBC, i.e., 61.11%. A possible explanation for this is because for pre-MBC, the MOD11A1 T s and T a temperature range were already close to each other and hence, the application of the MBC tended to overcorrect the MOD11A1 T s and result in higher magnitude of underestimation. Contrarily, the MBC only contributed a slight change for the post-MBC r values. This is because the application of the MBC technique could only adjust the range of the T s values by shifting the T s values towards the T a according to the CF values that satisfy the condition given (for LS-MBC) or according to the respective quantile calculation (for QM-MBC) without changing the original temporal pattern of the T s data. While these techniques have resulted in a significant improvement of error and bias terms for the T scd, T scm , and T scq against T a , the relationship between both datasets showed a moderate to no improvement.
Additionally, the assessment on post-MBC (T scd , T scm , and T scq ) showed that ranges for error and bias values were not as high as in the pre-MBC data. All post-MBC T scd , T scm , and T scq products depicted low and similar ranges of average bias and errors (average RMSE and MAE < 3 • C and average PBIAS between −5% and 5%), while the pre-MBC demonstrated that the MYD11A1/11A2 T s had higher error (average RMSE and MAE = 4.56 to 5.13 • C) and bias (PBIAS = 8.88 to 10.74%) range than the MOD11A1/11A2 (average RMSE and MAE = 2.88 to 3.48 • C, PBIAS = 1.36 to 4.40%). This suggested that both MBC approaches (LS-MBC and QM-MBC) were suitable for all types of MODIS T s tested in this study.
For comparison, Figure 11 displays the boxplot presentation for the average temperature values between pre (T s ) and post-MBC (T scd , T scm , and T scq ) against T a from 36 AWS. Significant improvements could be observed in the post-MBC values as the range of T s has decreased tremendously approaching the range of T a with T scq depicting the most approximate temperature range to the T a. This indicates that the QM-MBC method can capture well the distribution of T a by removing the biases by quantile [44] rather than the LS-MBC technique that removes the biases using CF values from overall constant daily and monthly values.  Figure 12 illustrates the boxplot presentation of evaluations metrics score between T s , T scd , T scm , and T scq against T a for 36 AWSs. The bias correction method consistently reduced the RMSE and MAE to 2-5 • C and 1-4 • C, respectively. The PBIAS also showed a reasonable reduction in the over-/underestimation to range of bias between −6.34% and 4.78%. In comparison, T scd , T scm , and T scq demonstrated a small-range of error and bias, but not for r. T scm showed highest improvement in r values than the T scd and T scq . This indicates that although these post-MBC products produced an acceptable result for T a approximation, the application of LS-MBC technique using the monthly CF was better than the constant daily CF (T scd ) and quantile mapping technique (T scq ). This is because application of monthly CF is calculated from a more detail calculation (by month) than the daily CF that is calculated using single CF for the whole 13 years. Application of monthly CF enable changes in temperature according to the changes in the monsoon season, if present, to be applied. Figure 13 displays the comparison between evaluation metrics between T s_r , T scd_r , T scm_r , and T scq_r against T a_r for 5 regions. The LS and QM-MBC has tremendously reduced the bias and error values as the boxplot for post-MBC (T scd_r , T scm_r , and T scq_r ) were located nearer to zero value than the pre-MBC (T s_r ). Meanwhile, the r values depicted modest improvement with the r maximum values pre-MBC (T s_r = 0.24-0.32) slightly increased (T scd_r = 0.25-0.31, T scm_r = 0.27-0.39, and T scq_r = 0.40-0.52).

By Region Performance
Moreover, the T scq_r depicted the most improved values of RMSE, MAE, and PBIAS. This indicates that QM-MBC can perform better when more data are stack together (regional data). With more collocated pixels, the temperature distribution can be better defined, resulting in better quantile calculation [32,33]. The same pattern was observed for r as T scq_r showed the highest improvement of r for MOD11A1, MOD11A2, and MYD11A2.
The low correlation coefficient (r) for both pre-and post-MBC data compared to T a , despite the low RMSE obtained for post-MBC are possible and acceptable. In this study, the reason for low r values could be because of the low temperature range for between 10 and 15 • C and less variation in the temperature pattern [16]. In the tropic region, there are only two distinct seasons, the hot and dry and cold and wet season driven by the two monsoonal wind seasons, i.e., the Southwest and Northeast monsoons. Besides, El Kenawy [15] expresses low correlation values for summer compared to other seasons. This indicates that solar radiation duration may affect the correlation between T a and T s . Figure 12. Comparison of evaluation metrics score between T s , T scd , T scm , and T scq against T a .
The information obtained from each station (total of 36 stations) and region (total of 5 regions) including (1) the validated range of T s /T s_r and T a /T a_r , (2) constant daily and monthly CF values (LS-MBC), (3) constant daily and monthly mean T s /T s_r and T a /T a_r (LS-MBC), and (4) quantile values for T s /T s_r and T a /T a_r (QM-MBC) can be used for future T a estimation from MODIS T s (MOD11A1, MOD11A2, MYD11A1, and MYD11A2) by following the MBC equation. The similar pattern of result between the regional (T s_r and T a_r ) and the by station (T s and T a ) analysis, proved that the information obtained from this study are applicable for the lowland areas according to their identified regions. Overall, the developed technique by modifying the CF application for LS-MBC is rather a simple one but managed to gain a good accuracy of T a retrievals, whereas the QM-MBC worked best for regional study.

Limitations and Future Works
First, we acknowledge that the findings of this study are only limited to lowland areas in Peninsular Malaysia at altitude less than 60 m above sea level. In the future, a study should be conducted for highland area given that the T s is known to be affected by the increase in elevation [22]. Second, the missing data resulted from cloud cover (T s ), partial cloud cover (T s ), and outlier's analysis (T s and T a ) in this study were treated as not available (NA). To overcome the gaps in T s data, interpolation [45], similar pixel method [46], data fusion [47], or machine learning [48] have been applied in other works but not in our study that only used original datasets. Finally, we acknowledge the limitation in delivering the correction factor values given unknown gaps in the constant daily and monthly T s , as opposed to Yoo et al. [49] who performed linear interpolation. Addressing further, the second and third gaps can possibly avoid distorted CF values and thus give better results. Figure 13. Comparison of evaluation metrics score between T s_r , T scd_r , T scm_r , and T scq_r against T a_r .

Conclusions
This study found that (1) the error values for pre-MBC T s in comparison to the T a were high especially for the MYD11A1/11A2 than the MOD11A1/11A2; chiefly due to different acquisition times that is influenced by solar activity. (2) The high bias and error values for the pre-MBC T s suggested that the application of correction is necessary. (3) The MBC for station analysis (T scd , T scm , and T scq ) has successfully reduced the errors and biases in pre-MBC (T s ), i.e., average RMSE values from 3.42-5. 13 (5). Both MBC methods (LS-MBC and QM MBC) resulted in an acceptable range of correction but by a comparison, the MBC using monthly CF depicted a slightly better correction than the MBC using daily constant CF, while QM-MBC worked best for regional study. This study also illustrated that although the proposed LS-MBC technique in this study is rather a simple mathematical equation than the more sophisticated QM-MBC technique, it was adequate for improving the estimation T a from MODIS T s for a tropical region.

Data Availability Statement:
The data presented in this study are available and freely downloadable from https://search.earthdata.nasa.gov/search (accessed on 2 January 2020). The other data are available on request from the corresponding author.