Optimized AVHRR land surface temperature downscaling method for local scale observations: case study for the coastal area of the Gulf of Gdańsk

2.1 Downscaling theory

Downscaling is the process of converting a set of spatial information (e.g. 2- or 3-dimensional) from a low to a high resolution. In the context of LST estimation based on satellite imagery, downscaling refers to estimation of LST values within the sub-pixel area. In practice, in order to downscale a satellite image, thus improving its spatial resolution, the approach of merging at least two datasets is used; a low resolution image is merged to a high resolution image (from the same or different sensor), or a low resolution image obtains spatial details of a high resolution image [24, 25].

According to the downscaling principles, the proposed approach assumes to retain radiometric characteristics of input images and aims to minimize the error. Let us assume two input data sources - such that both sources represent acquisition over the same area, and that a physical quantity they measure is comparable, i.e. they measure the same quantity in the same units (energy, radiance, temperature, etc.). Another point in the methodology is the definition of relation between two data sources that enables transformation from one source domain to another. Both sources have different but spatially constant resolution (Fig. 1).

Figure 1

Diagram representing the relation between the spatial distribution of measured quantity in high resolution and low resolution

Let us assume that for a given flat area, corresponding to a pixel in coarser resolution, the value of measured quantity is the weighted average of this quantity for all subareas corresponding to pixels in finer resolution contained by the coarse resolution pixel area. It is mathematically expressed as:

where

where

T_low - is the value of the pixel in coarse resolution,

Thigh(n) - is the value of the n-th pixel in high resolution image belonging to the mentioned high pixel resolution area,

s⁽ⁿ⁾ - is the area of the given pixel in high resolution,

S - is the area of a single pixel in coarse resolution,

N - is the number of pixels in high resolution imagery that are contained in the corresponding pixel of low resolution image.

In general, N depends on the ratio of data sources’ spatial resolution taken under consideration and may also depend on the number of high resolution pixels being excluded from the analysis for a given low resolution pixel due to cloud presence and other limitations.

The dependency mentioned above enables transformation between the high and low scale. However, as the transformation of pixel values from the high to low scale results in the loss of information, the reverse transformation is obviously ambiguous. The estimation of pixel values of high resolution having given pixel values of low resolution requires some additional, local, high scaled information about characteristics of the specific area that influences the pixel values T_low and T_high, i.e. some independently measured quantity Q usually derived from satellite imagery or auxiliary database as well. The influence of Q on T ought to be stable, which means that the dependence between Q and PV should be the same for the entire investigated area. Also, Q values should be stable in time (to a higher degree than the quantity expressed by PV), which makes it easier to collect them (i.e. lower satellite revisit time is needed). In the case of LST as T, the surface thermal properties expressed by its vegetation index or effective emissivity may be used as Q [26]. The latter was utilized by the authors in the current study. One of the main assumptions applied is that the generalization of the rule connecting the local emissivity and LST value must be allowed with respect to at least a single scene of low resolution imagery.

The way of utilizing the high resolution information on thermal properties of the land for downscaling the LST image is described in the next subsections.

2.2 PBIM method

Having defined the concept of downscaling, let us focus on the particular method of downscaling to be optimized. Since the objective of the study is to present the optimized downscaling AVHRR LST method in view of its applicability in coastal zone LST estimation, we have taken under closer analysis the PBIM (Pixel Block Intensity Modulation) downscaling process, described in the next subsection, and proposed a way of its optimization.

PBIM is proposed for downscaling thermal products to higher resolution using the high resolution transformation kernel. The general form of the PBIM method is expressed as:

where T_low - is the pixel value of the image in low resolution, ε_high - is the high effective emissivity pixel value, ε_high→mean - is the high effective emissivity value averaged for a current low resolution pixel area, T_high - is the resulted current pixel value of high resolution downscaled image. It should be noted that according to formula (4), the strict linear relation between ε and T is assumed.

The appropriate definition of the used high resolution data source, referred to here as ε, is the most important issue in the utilized downscaling method. The better the description of local properties of the area the image of which is to be downscaled, the smaller the obtained error of transformation. The effective emissivity composite product values calculated based on several high resolution images were utilized as the HRK [27].

The PBIM allows for merging multiresolution 2D datasets retrieved from sensors. Note that in this approach, the simulated high resolution image retains the radiometric characteristics and spatial distribution of the original low resolution image. In particular, the average of the simulated high-resolution sub-pixel values is equal to the pixel value of the original low resolution image.

The PBIM method can be applied to any satellite imagery derived products - from a relatively low level, where pixel values are spectral bands of the images, to higher levels like LST or SST products.

Note also that in this approach to re-scaling of thermal products it is important that spatial data considered for processing must represent relatively flat areas, thus the low scale can be expressed as a simple aerial average of the LSTs at the high scale [28]. Also, the re-scaling methods have limitations for areas with high surface type diversity, which particularly refers to wetlands, water containers situated among sub pixel areas. Water surfaces have different emissivity and heat capacity properties and thus they significantly influence the model [29].

2.3 Proposed method of PBIM optimization

In the presented approach, as in the PBIM method, the ε_high also plays a principal role, as it describes high resolution properties of an area corresponding to LST image to be downscaled. As mentioned, ε_high is computed based on high resolution quasi-static data that describes thermal properties of an area. The choice on what should be the basis for the construction of ε_high used for defining the rescaling transformation is somewhat arbitrary. The general rule is that it should possibly describe the thermal properties of the surface in the best possible way - particularly the surface’s susceptibility to absorb heat from shortwave radiation if the aim is to downscale day-time LST [30]. Also, what has been introduced in subsection 2.1, ε_high is considered to be quasi-static and much easier to be retrieved than data to be downscaled [31]. Note also that images used to construct composite data must be registered in similar vegetation and weather conditions [32]. In our approach, we have taken under consideration the summer period of the years of 2013 and 2014. The FVC and subsequently the effective emissivity ε have been calculated on the basis of the normalized differential vegetation index (NDVI) using the formulae (20) and (21) from subsection 3.2.

What is more, as mentioned in the previous section, the PBIM method relies on a relatively simple assumption that T_high depends linearly on the respective values (eq. 4). It can be observed that this is not always the case [33]. For instance, Fig. 2 presents the plot of effective emissivity from Landsat 8 OLI vs. corresponding values of LST from Landsat 8 TIRS for the investigated Pomerania region in Northern Poland (for details on data processing, and specifically on effective emissivity and LST calculating from Landsat 8 imagery – see section 3). The content of Fig. 2 is in general in line with the known observations of the “triangle” character of the VI – LST dependence. Fig. 2b shows the relation between effective emissivity values, namely, the central values of particular small ranges of effective emissivity, and LST values averaged for groups of pixels corresponding to those effective emissivity small ranges (red points). It is able to be seen from the figures that the relation between effective emissivity and LST tends not to be linear. The line of fit (corresponding to all data presented in Fig. 2a, not only to averaged data from Fig. 2b) is shown in Fig. 2b by black dots. It is visible that it cannot be used as a satisfactory approximation of ε_high - LST dependence.

Figure 2

LST/effective emissivity scatter plot for considered area (a). Right picture (b) represents effective emissivity and LST dependence for considered region: red plot vs. LST value represents averaged for groups of pixels corresponding to small ranges of effective emissivity. Black plot is a linear approximation of effective emissivity and LST dependence.

To take into account the non-linear dependence between the biophysical properties of the surface in the LST downsampling, the authors propose the following inverse downscaling approach (DS_opt). We assume that there exists a dependence between effective emissivity and LST expressed as

where f (·) need not to be linear (and does not have to be monotonic as well), and may be estimated for instance as the red plot in Fig. 2. Combining (2) and (5), and naming PV explicitly as LST, yields:

for a given low scale pixel. The form of f (·) is generally unknown for a given terrain, but assuming that it is stable for the entire scene and applying (6) to each low scale pixel, we may express the problem as a system of linear equations. For this purpose, let us calculate the histogram [h₁, h₂, ..., h_K] of ε_high (where h_k, k = 1, ..., K is the percentage of high scale pixels whose ε_eff value belongs to a given range) for each low scale pixel, using always the same ε_eff minimum/maximum values and the same number of bars K. Then, if we express f (·) as a vector of values (weights) [w₁, w₂,...,w_K] corresponding to the central values of ε_eff ranges used in calculating the histogram, and assuming that the single pixel area S_n is constant, we may write (6) as

where m is the number of a “coarse” pixel in a scene (m = 1,..., M). The system of M equations (7) may be also written in a matrix form:

and may be shown as a diagram – Fig. ??.

Referring to (7) and (8), each row in H matrix represents information on thermal properties of terrain within a certain pixel of low resolution image (Fig. 3). Row values of H represent histogram values calculated using a corresponding set of N pixel values of ε_eff from high resolution image (in reality, usually less than N pixel values is used for a given coarse pixel, due to cloud presence or other limitations in high scale data from which ε_high has been derived). LST_low is a column vector containing low resolution image LST values and w is the vector of weights to be found. Having two data sources (LST in a coarse scale and ε_high in a fine scale) which fulfil the conditions stated above, for the whole scene a set of M equations with K variables is obtained.

Figure 3

Diagram of matrix equation representing dependence between high and low resolution datasets

Finally, as w_k values express the LST themselves, after finding them with use of the method described in the next subsection, the LST_high may be estimated directly as w_k corresponding to a given ε_high value, or with the application of some kind of interpolation of w values.

2.4 Finding an optimal solution

Now our task is to estimate w having given LST_low and H. Let us denote the observation vector LST_low as y and the searched vector of parameters w as x, and include the vector v of bias which is always present in a real measurement situation:

Estimation of x having known H and measured y biased by unknown v is a typical linear inverse problem [34]. Unfortunately, in reality we usually have to take into account that the postulated model does not describe the relation between x and y perfectly, e.g. it is linear only in approximation, or the assumed values in H are not fully correct, or the assumed number of independent variables (the number of x elements) is not correct. Moreover, the K and M values are not equal and thus H matrix is not invertible. Thus, one of the possible ways to obtain the solution in such a situation is to estimate x in a sense of least squares, i.e. by minimizing the ∥Hx − b∥² where ∥ · ∥² is the Euclidean norm. Estimate x^ may be then calculated as:

if the H^TH matrix is not singular. In such a case, a regularization term is introduced to provide the existence of the solution, changing (10) to:

where Γ is a suitably chosen Tikhonov matrix [35-38] and x₀ is the expected solution calculated from the PBIM method. We used the simple form of Γ = λI, where I is the identity matrix and the λ value was found using the nearoptimal parameter calculation of singular value decomposition method proposed by O-Leary [34].

The solution x^ is the optimal vector that minimizes the error of reconstruction of simulated high resolution image based on low resolution AVHRR imagery and utilized HRK that is stored in H matrix.

At the last step of the proposed approach, the image is reconstructed to high resolution using optimal weights stored in x^. Reconstruction is performed from x^a> vector using local low resolution LST correction, denoted as DS_opt:

where: - (i¨˙,j¨) denotes low resolution pixel coordinates that correspond to high resolution (i, j) coordinates, - ε_high - is the effective emissivity.

3.1 The experiment scheme and data flow

In order to compare the results obtained using the proposed approach with other existing solutions, we used two basic types of datasets: AVHRR and Landsat 8, registered at possibly the same time. Then, we assumed that AVHRR LST is a low spatial resolution (1.1 km) input data that is downscaled to higher resolution (100 m) and we treat Landsat 8 LST product as a kind of proxy validation dataset.

In this case, the effective emissivity, as the static high resolution data to be merged with low resolution image, was retrieved from a composition of two Landsat 8 images registered on 2013-08-05 at 09:45:37 UTC (scene center time) and on 2014-07-23 at 09:43:29 UTC (scene center time). For both of the scenes, fractional conditions for the analyzed area are similar (middle of the summer) and repeatable. The Metop-B/AVHRR3 imagery was registered on 2013.08.05 on 9:10:26 UTC.

The overview diagram of processing and validation strategy of the proposed approach is presented in Fig. 4. After creating the composite FVC high resolution image from two Landsat 8 scenes, it was used along with the AVHRR low resolution LST image as an input to our downscaling procedure. Also, for comparison of the results, the PBIM downscaling was applied for the same data. At the end, a comparison of the down-scaling results, obtained both by DS_opt and PBIM with the high resolution Landsat 8 LST data, was performed.

Figure 4

Overall diagram of data processing flow and verification

It should be noted that for Landsat 8, in the climatic zone of the considered region several requirements need to be fulfilled to provide the appropriate analysis. The acquisition should be made under clean sky conditions in order not to interfere with the values registered in thermal and visible channels. The presence of clouds, aerosols, and diverse spatial distribution of water vapor content highly disturbs the process of acquisition. Also, note that Landsat 8 was launched in 2010 and together with a long revisit time (16 days), it highly reduces the number of available acquired imageries where all conditions stated above can be met.

We also performed accurate geometric calibration of both data sources to provide a unified reference system and cover the same areas of analysis. Data from both sources were reprojected to UTM zone 34 projection where the imageries cover area from 271872 to 434655 easting and 6004014 to 6092171 northing based on WGS 84 ellipsoid. The whole area covers 14,256 km², in which land consists of about 51% of the data. The height of the terrain is not moderately diverse as most of it is covered by flat areas or small hills. Diverse vegetation conditions of the urban areas of Gdańsk, Sopot and Gdynia municipality (Tricity), and green areas (forests, farmlands) located nearby, cause the NDVI to vary from approx. 0.12 to 0.67.

Fig. 5 presents the map of the location of the area of interest: the land at the Gulf of Gdańsk, Southern Baltic Sea in Poland, along with LST maps calculated using imageries from both sensors – high resolution Landsat 8 and low resolution MetopB/AVHRR.

Figure 5

On the left: the map depicting the location of the area of interest – the Gulf of Gdansk, South Baltic, Poland (source: Google Earth), on the right: LST maps calculated from both sensors – Landsat 8/TIRS (top-right) and MetopB/AVHRR3 (bottom-right). Blue pixels in the upper and the lower right images represent sea or cloud areas that were masked out of the analysis.

3.2 LST estimation

In the study, we selected MetOp-B/AVHRR3 imageries as the source of low resolution data to be downscaled and Landsat 8/TIRS LST as a validation dataset.

The set of NOAA and MetOp satellites utilizes the AVHRR radiation-detection moderate resolution remote imager dedicated to land and atmosphere. Due to quite a large number of currently operating NOAA and MetOp satellites, the AVHRR sensor is characterized by a relatively high time resolution (the revisit time for the respective region is approx. 2-3 hours). Tab. 1 presents the characteristics of AVHRR channels according to NOAA KLM User’s Guide.

As it is known, LST estimation requires emissivity correction [39-42]. Therefore, during this process we calculate the average emissivity for 4^th and 5^th AVHRR channel, denoted as ε_AVHRR and the difference of these channels’ emissivity Δε_AVHRR. The calculation process was different for several NDVI ranges. Namely, for NDVI < 0.2 the pixel is considered to be bare soil, then:

where R1 is AVHRR channel 1 reflectance.

Pixels with 0.2 < NDVI < 0.5 are considered to be partly vegetated and the ε_AVHRR and Δε_AVHRR are calculated using the following formulae:

where

For NDVI > 0.5 pixel is considered to be full vegetation and ε_AVHRR and Δε_AVHRR is fixed as:

The validation dataset was retrieved from the Landsat 8 satellite, which images the entire Earth every 16 days in an 8-day offset from Landsat 7. Landsat 8 carries two instruments: the Operational Land Imager (OLI) sensor that includes refined heritage bands, along with three new bands with respect to Landsat 7 [43]a deep blue band for coastal/aerosol studies, a shortwave infrared band for cirrus detection, and a Quality Assessment band. The Thermal Infrared Sensor (TIRS) provides two thermal bands. Both these sensors provide improved signal-to-noise radiometric performance quantized over a 12-bit dynamic range. The improved signal to noise performance of Landsat 8 enables better characterization of the land cover state and conditions. Products are delivered as 16-bit images (scaled to 55,000 grey levels). Table 2 presents the characteristics of Landsat 8 channels.

The LST from Landsat 8 imagery was calculated by applying the split window technique algorithm that uses two thermal bands located in the atmospheric window between 10 and 12 µm. Then, we applied radiometric calibration procedure for Landsat 8 data using metadata of the provided imagery and used expansion of the Planck equation according to (17) for TIRS data. Spectral radiances delivered in Landsat datasets were converted into satellite brightness temperature using the following relationship [44, 45]:

where:

K₁ - band specific thermal conversion coefficient delivered with Landsat 8 imagery metadata files, namely K₁ = 774.89 for channel 10 and 440.89 for channel 11 [Wm2sterμm],

K₂ - band specific thermal conversion coefficient delivered with Landsat 8 imagery metadata files, namely K₂ = 1321.08 for channel 10 and 1201.14 for channel 11 [Wm2sterμm]˙

TOA_i - top of the atmosphere spectral radiance [Wm2sterμm] registered for i-th thermal channel of Landsat 8 TIRS.

Landsat 8 LST emissivity correction, denoted as ε_TIRS, was performed according to a procedure proposed by Skokovic et al. [46], Shaouhua Zhao et al. [47], according to the values presented in the Tab. 3, with use of the Corine Land Cover database and FVC data, in order to distinguish urban and non-urban area.

The split window technique (SWT), utilizing the thermal channels of AVHRR3 and TIRS, has been used here for LST calculation. This is a relatively simple and robust method of retrieving LST estimation for large or regional scale applications. Its main property is that it removes atmospheric fluctuations, taking advantage of the fact that radiation absorption differs in particular wavelengths of electromagnetic spectra. Consequently, the TIRS and AVHRR LST images were corrected pixel by pixel for emissivity and water vapor content, according to Jimenez-Munoz [48–51] with values of c₀-c₆ (Tab. 4) according to the following formula:

where W is a water vapor content retrieved from an external database [50–54], T_i and T_j are thermal channels of AVHRR (i = 4 and j = 5) and TIRS (i = 10 and j = 11) respectively, ϵ^ is the average emissivity for thermal channels of AVHRR3 and TIRS, respectively (ε_AVHRR and ε_TIRS) and Δε is the emissivity difference between the thermal channels of AVHRR3 and TIRS.

With the above consideration in mind, in our methodology we proposed the following approach:

1. we applied the split window technique for MetopB/AVHRR3 and for Landsat 8 TIRS,

2. in both algorithms we used the surface emissivity and water vapor content correction approach.

Effective emissivity ε computation was based on high resolution NDVI data over the respective region. The maximum NDVI (NDVI_max) and the minimum NDVI (NDVI_min) values for the study area were determined, which were further used for computing the fractional vegetation cover (FVC) [57, 58]:

and the effective emissivity ε as:

4.1 Comparison of LST downscaling results for the entire scene

In this subsection we present the results obtained using the proposed optimized downscaling process of LST. We compared the results obtained by DS_opt with the PBIM method of downscaling in global (i.e. for the entire scene Fig. 6, and using a more detailed analysis in Fig. 7) and local scale (i.e. for the locations along the two selected transect lines within the investigated land area, Fig. 8, 9 and 10) using objective criteria, namely, the root mean squared error (RMSE) and the Pearson correlation coefficient (R). Pixels representing either clouds or sea on any of the images were masked during the analysis so their values did not interfere with the overall process of error minimization.

Figure 6

Comparison of downscaling results for the investigated area: top-left: LST calculated from Landsat 8 image, top-right: LST downscaled from AVHRR using PBIM method, bottom-right: LST downscaled from AVHRR using DS_opt method

Figure 7

The obtained downscaled to high resolution AVHRR LST image and scatter plot (every 100th point) of Landsat 8 LST and downscaled AVHRR LST (right) using DS_opt. Red line in right image represents “1:1” line. Correlation coeflcient and root mean squared error (RMSE) were: R=0.828 and RMSE=2.255°C, respectively.

Figure 8

Landsat LST image with marked two transect lines defining the small subareas selected for detailed analysis of the downsampling results. The first line (AB) starts from suburban areas of the Tricity municipality (A), passes through a vegetated area of the Tricity Landscape Park and ends in the Tricity municipal area near the border between Gdańsk and Sopot (B). The second transect starts from the suburban area of Gdynia (C), passes through forest areas and ends in the center of Gdynia city (D).

Figure 9

Comparison of LST high resolution data with respect to AB transect line: Landsat 8 LST, AVHRR LST downscaled by DS_opt method, AVHRR LST downscaled by PBIM method and AVHRR LST downscaled by nearest neighbor technique

Figure 10

Comparison of LST high resolution data with respect to CD transect line: Landsat 8, AVHRR DS_opt downscaled by the proposed method, AVHRR downscaled by the PBIM method (AVHRR PBIM) and AVHRR resampled to LANDSAT 8 resolution using nearest neighbor technique

We compared the obtained downscaled high resolution image to LST calculated from Landsat 8 image acquired on 13 August 2013. The RMSE between Landsat 8 LST and downscaled image using DS_opt was 2.255°C with correlation coefficient R equal to 0.828 and Bias = 0.557°C. It shows a significant improvement in comparison to the PBIM method, where we obtained RMSE = 2.832°C, R = 0.775 and Bias = 0.997°C.

The analysis of spatial distribution images shows that the results obtained by DS_opt algorithm better represent thermal characteristics for the analyzed area and are visually closer to the original Landsat 8 LST image than the image obtained using the PBIM method (Fig. 6). Specifically, it can be observed that the influence of local surface high resolution characteristics is reflected in DS_opt rather than local LST retrieved from AVHRR observations. What is more, as in the PBIM, the LST_low term is directly used in the formula and strongly influences the obtained LST_high value, it results in the appearance of coarse resolution artifacts in the downscaled product (Fig. 6b). This effect does not occur in the DS_opt result case. In Fig. 7 we presented the scatter plot of this case.

4.2 Detailed analysis of the retrieved results

In many cases, for instance for small areas where we deal with high diversity of surface conditions, obtaining the local accurateness of the applied downscaling is even more important. For this purpose, we took under consideration two relatively small sub-areas defined as narrow strips of land situated along selected “cross section” lines, in order to verify how the algorithm manages with diverse heterogeneous conditions (Fig. 8).

The first line starts from the West, at point A located at 53.9356°N, 17.6832°E, which is situated in the suburbs of the Tricity area. Then, the line transects through the Tricity Landscape Park and ends in the urban area of the Tricity municipality (Gdańsk-Oliva district) - point B (53.9309°N 17.8077°E). The length of the line is 6 km and it transects through low as well as very high vegetated areas, so the thermal conditions along this transect are significantly diversified.

The second line starts from suburban area of Gdynia (point C – 54.045 °N, 17.6513 °E), passes through forest areas and ends in the center of Gdynia city (point D – 53.9544 °N, 17.6601 °E). The diversity of thermal conditions for the area along the second transect is quite similar to that in the case of the first line.

Afterwards we performed analysis of the obtained results along selected transects from four datasets, namely:

1. LST calculated from Landsat 8 TIRS dataset (validation dataset),

2. AVHRR LST image downscaled using the proposed DS_opt method,

3. AVHRR LST image downscaled using the PBIM method,

4. AVHRR LST image downscaled transformed (resampled) to high resolution by simple nearest neighbor method.

Plotted LST high resolution values taken from particular datasets along the AB and CD transect, respectively are presented in Fig. 9 and 10. Tables 5 and 6 contain the comparison of the correlation coefficient and RMSE between Landsat 8 LST and the results retrieved by different methods of AVHRR LST downscaling: DS_opt proposed by the authors, PBIM, and nearest neighbor for AB and CD transect, respectively.

For this relatively small area, the above results show that the optimization method presented in DS_opt algorithm yields significantly better results than other approaches. Note that the aim of using the inverse problem algorithm is to minimize the downscaling error in a global, i.e. the entire imagery scale, so results corresponding to local sub-datasets do not always have to be better than when using other methods. However, this case shows that data obtained by the DS_opt fits the original Landsat 8 LST the best, which serves here as a kind of proxy verification tool.

As shown above, the obtained downscaled products retain radiometric characteristics of original images for the whole imagery, however, for subdatests the BIAS results suggest that some deviation from this rule occurs. Moreover, the correlation coefficients and RMSE indicate that DS_opt technique yields results closer to the original highresolution spatial pattern than the PBIM method. A more detailed investigation of the selected transect cross sections of the analyzed area indicates that also for small sub-datasets the usage of DS_opt is justified as the RMSE error is reduced from 2.6851°C for the PBIM case to 2.0042°C for AB transect and from 2.5262°C to 1.5401°C for CD transect. The improvement of the proposed methodology is further verified by correlation between downscaled products and time-coincident Landsat 8 imagery. Specifically, the correlation coefficient increases from 0.8207 to 0.9212 in the case of AB, and from 0.6038 to 0.8953 in the case of CD transect.

4.3 Optimized downscaling for local scale

In the previous subsection, we have examined and evaluated the results of downscaling using the inverse technique for the whole scene imagery and for particular transect sub-datasets. We have shown that the proposed optimization technique yields better results in comparison with the direct approach. In this section, we have taken under consideration three sub-datasets of the whole scene, which differ between each other in terms of vegetation and surface conditions, and we have verified how our approach functions in such cases. The cases taken under consideration are as follows:

Case 1 - area covering nearly 300 km² of Gdynia municipality urban area and Tricity National Park bounded by 53.7451^°N, 18.4904^°E (north-west) and 53.5237^°N, 18.6786°E (south-east) limits. This area is rather unique and is to large extent characterized by surface type heterogeneity typical for urban area. According to the Corine Land Cover database, 51% of the mentioned area is urban, while 49% is non-urban, high vegetation area of the Tricity National Park. The results obtained by PBIM and DS_opt methods for test case 1 are shown visually in Fig. 11. The quantitative results on the obtained correlation coefficient, the RMS error and the bias for this case are shown in Table 7.

Figure 11

Scatterplot and visual presentation of the results retrieved by PBIM and DS_opt methods for local sub-dataset test case 1

Case 2 – area of the Stęzyca Lake which is a part of the Pomeranian National Park that mainly comprise vegetated areas (forests and grasslands), which constitute 99% of the surface. The area is bounded by coordinates: 53.5222°N, 17.9436°E (north-east) and 53.3807°N, 18.3026°E (south-west). The results obtained by PBIM and DS_opt methods for test case 2 are shown visually in Fig. 12. The quantitative results on the obtained correlation coefficient, the RMS error and the bias for this case are shown in Table 8.

Figure 12

Scatterplot and visual presentation of the results retrieved by PBIM and DS_opt methods for local sub-dataset test case 2

Case 3 - urban area of Gdańsk municipality, bounded by 53.5281°N, 18.5208°E and 53.3891°N, 19.1395°E, that mainly comprise urban areas (over 80% according to the Corine Land Cover database). The results obtained by PBIM and DS_opt methods for test case 3 are shown visually in Fig. 13. The quantitative results on the obtained correlation coefficient, the RMS error and the bias for this case are shown in Table 9.

Figure 13

Scatterplot and visual presentation of the results retrieved by PBIM and DS_opt methods for local sub-dataset test case 3.

Tables 7-9 present analytical results obtained by the proposed optimization approach, PBIM method and original AVHRR3 data resampled using the nearest neighbor technique (presented in the third row) for cases 1-3. The second column of the tables represents the Pearson correlation coefficient, while the third represents the RMSE with respect to reference LST Landsat 8 OLI dataset. The last column shows BIAS between reference dataset and a result of downscaling products, in order to verify retention of radiometric properties of the results.

The results presented in these cases show that DS_opt method yields results closest to Landsat 8/TIRS validation dataset having overall lowest RMSE and highest correlation coefficient among the methods presented. This can be particularly observed for test case 1 (Tab. 7) when algorithm finds minimized error solution for a heterogeneous area, including high vegetation (low effective emissivity) and urban areas (high effective emissivity), yielding RMSE reduction up to 32% (for DS_opt) and correlation increase from 0.7659 (for PBIM) to 0.8866 (DS_opt). As this trend is also observed in cases 2 and 3 (Tab. 8 and Tab. 9, respectively), it may be concluded that the increase of the proposed downscaling method performance is visible for local sub-datasets.