Landsat 15-m Panchromatic-Assisted Downscaling (LPAD) of the 30-m Reflective Wavelength Bands to Sentinel-2 20-m Resolution

The Landsat 15-m Panchromatic-Assisted Downscaling (LPAD) method to downscale Landsat-8 Operational Land Imager (OLI) 30-m data to Sentinel-2 multi-spectral instrument (MSI) 20-m resolution is presented. The method first downscales the Landsat-8 30-m OLI bands to 15-m using the spatial detail provided by the Landsat-8 15-m panchromatic band and then reprojects and resamples the downscaled 15-m data into registration with Sentinel-2A 20-m data. The LPAD method is demonstrated using pairs of contemporaneous Landsat-8 OLI and Sentinel-2A MSI images sensed less than 19 min apart over diverse geographic environments. The LPAD method is shown to introduce less spectral and spatial distortion and to provide visually more coherent data than conventional bilinear and cubic convolution resampled 20-m Landsat OLI data. In addition, results for a pair of Landsat-8 and Sentinel-2A images sensed one day apart suggest that image fusion should be undertaken with caution when the images are acquired under different atmospheric conditions. The LPAD source code is available at GitHub for public use.


Introduction
The Landsat-8 polar-orbiting satellite (launched 2013) carries the Operational Land Imager (OLI) that has nine reflective wavelength (435 nm to 1384 nm) spectral bands; eight 30-m visible, near-infrared, short wave infrared, and cirrus bands, and also one 15-m panchromatic band [1].The multi-spectral instrument (MSI) onboard the Sentinel-2A (launched 2015) and Sentinel-2B (launched 2017) polar-orbiting satellites acquires multi-spectral reflective wavelength medium resolution data with a design heritage from the Landsat and Satellite Pour l'Observation de la Terre (SPOT) sensors [2].The MSI has 13 reflective wavelength (443 nm to 2190 nm) spectral bands; four 10-m visible and near-infrared bands, six 20-m red edge, near-infrared and short wave infrared bands, and three 60-m bands.Combination of near-contemporaneous Landsat-8 OLI and Sentinel-2 MSI data will provide increased opportunities for more frequent cloud-free surface observations [3,4].For example, recent research suggests that the integration of the Landsat-8 and Sentinel-2 data would be beneficial for applications including agricultural mapping and monitoring [5][6][7], phenological studies [8], and glacier extent mapping [9].There are a number of pre-processing issues that need to be addressed before the well calibrated Landsat-8 [10,11] and Sentinel-2A [2,12] data can be used together or treated as effectively being sensed from the same sensor.These include handling the different sensor spectral response functions and correction for atmospheric effects [13,14], correction of surface reflectance anisotropy [15,16], and handling image tiling and geolocation differences [17][18][19].
This paper is concerned with the Landsat-8 and Sentinel-2 spatial resolution differences.The different Landsat-8 and Sentinel-2 spatial resolutions can be overcome by reprojecting their data into the same coordinate system and then upscaling or downscaling the data from one sensor to the other.The Sentinel-2 10-m and 20-m data could be upscaled to match the Landsat-8 30-m data.For example, by resampling the 10-m and 20-m Sentinel-2 bands to Landsat 30-m resolution for sensor comparison of spectral indices [20] or to derive between-sensor registration information [21].Conversely, the Landsat-8 30-m data could be downscaled to match the Sentinel-2 20-m or 10-m data.For example, by resampling the Landsat-8 30-m data to Sentinel-2 resolution for applications including agriculture monitoring [7], land cover classification [22], glacier extent mapping [9], or to derive between-sensor registration information [18].Downscaling of Landsat-8 30-m data to the Sentinel-2 20-m or 10-m resolution is attractive as the more spatially resolved Sentinel-2 information is not lost.
In this paper a new method to downscale the Landsat-8 30-m bands to Sentinel-2A 20-m resolution is presented and compared to conventional resampling approaches.The Landsat 15-m Panchromatic-Assisted Downscaling (LPAD) method has two principal steps: (i) downscaling the Landsat-8 30-m reflective band data to 15-m resolution using the detailed spatial information provided by the Landsat-8 15-m panchromatic band and then (ii) reprojecting and resampling the downscaled 15-m data into registration with the Sentinel-2A 20-m data.It is shown to introduce less spectral and spatial distortion and to provide visually more coherent data than conventionally resampled Landsat 20-m data for pairs of images selected on the same day over diverse environments.In addition, results for Landsat-8 and Sentinel-2A images sensed one day apart are included and suggest that fusion of these sensor data should be undertaken with caution when they are acquired under different atmospheric conditions.The LPAD source code is made available at GitHub for public use.

Data and Study Areas
The Landsat-8 satellite was launched February 2013 and carries the OLI and Thermal Infrared Sensor (TIRS) that sense the Earth's land surface with a 16 day repeat cycle over a 185-km-wide swath acquired from a sun-synchronous polar orbit [1,23].The TIRS has two 100-m thermal bands but unfortunately a stray-light issue means that they are sub-optimal and land surface temperature retrieval using split-window techniques is not recommended [24].The Sentinel-2A MSI was launched June 2015 and the follow on Sentinel-2B MSI was launched March 2017 and each has a 10-day repeat cycle and together a five-day repeat cycle and provide 290-km-wide swaths from sun-synchronous orbits [2].Both the OLI and MSI sensors have 12-bit radiometric resolution.
Table 1 shows the OLI and spectrally equivalent MSI reflective wavelength bands used in this study.The shortest wavelength sensor blue bands were not used because they are particularly sensitive to atmospheric effects [13,25].The 30-m OLI cirrus (1363-1384 nm) band and the three 60-m MSI bands (433-453 nm, 935-955 nm, and 1360-1390 nm) were not used because they are primarily designed for atmospheric remote sensing [1,2].Five pairs of contemporaneous Landsat-8 and Sentinel-2A images over sites with different land cover and containing anthropogenic and natural features with straight and irregular boundaries were selected (Table 2, images illustrated in Section 4).The five sites include urban areas, crop fields, burned areas, landslides, and tropical forest, and were sensed under cloud-free conditions.Four of the image pairs were selected so that the Landsat-8 and Sentinel-2A images were sensed on the same day (from only about 5 min to 19 min apart) in order to minimize the impact of land surface and atmospheric changes.Consequently, none of the image data were atmospherically corrected.The Congo image pair was selected because the sensor data were acquired one day apart and the Sentinel-2A image contained evident smoke aerosols due to biomass burning.The impact of land surface and atmospheric changes on the Congo image pair is discussed in Section 6.The recently available Landsat-8 Collection 1 Tier-1 OLI data, that are geolocated with a radial root mean square error <12 m and have a band-to-band registration accuracy less than 4.5 m [1,26], were used.The data are available in approximately 185 km × 180 km scenes defined in the Universal Transverse Mercator (UTM) projection and in the heritage Worldwide Reference System (WRS-2) of path (ground track parallel) and row (latitude parallel) coordinates [27].The stored Landsat-8 digital numbers were converted to top of atmosphere (TOA) spectral reflectance (ρ λ , unitless) using metadata conversion coefficients and then dividing by the sine of the scene center solar elevation angle (θ s ).Sentinel-2A L1C data that are available in 109 km × 109 km tiles with the UTM projection and US Military Grid Reference System (MGRS) [17,28] were used.The Sentinel-2A L1C data have geo-location accuracy specifications of less than 12.5 m and band-to-band registration accuracies of less than 3 m, 6 m, and 18 m for the 10-m, 20-m, and 60-m bands, respectively [29].

Landsat 15-m Panchromatic-Assisted Downscaling (LPAD)
The LPAD method consists of two principal steps: (i) downscaling the Landsat-8 30-m reflective band data to 15-m resolution using the detailed spatial information provided by the Landsat-8 15-m panchromatic band and then (ii) reprojecting and resampling the downscaled 15-m data into registration with the Sentinel-2A 20-m data.These steps are described below.Downscaling is a process to enhance the spatial resolution of sensor data usually by using information provided by other higher resolution ancillary data [30][31][32][33][34].It is sometimes termed data fusion, or more specifically, pansharpening, when the higher resolution ancillary data is the panchromatic band [35][36][37][38][39]. Hereafter, for consistency, the term downscaling is used throughout this paper.Downscaling using panchromatic bands has been studied for several decades and the detailed critical surveys can be found in [40][41][42].Such methods are often grouped into component substitution (CS) and multiresolution analysis (MRA) methods [36,41,43,44].
In this study an implementation of the CS Brovey method [45] that was designed for Landsat-8 data and is computationally efficient and so appropriate for global application is implemented [46].The method includes three steps: (1) resampling all the 30-m Landsat-8 bands to 15-m resolution; (2) derivation of a 15-m Landsat-8 intensity image; (3) downscaling of all the resampled 15-m band values in (1) by modulating and adding the difference (the spatial details) between the derived intensity image in (2) and the 15-m panchromatic image.These steps are described below.
Step 1: All the 30-m Landsat-8 bands used in this study (Table 1) are resampled to 15 m resolution to align precisely with the Landsat-8 15-m panchromatic band.Unlike many commercial sensors, the 15-m band is not perfectly nested within the lower resolution 30-m bands [46], rather, the panchromatic pixel grid is shifted by 7.5 m in the row and column directions relative to the 30-m pixel grid (Figure 1).values in (1) by modulating and adding the difference (the spatial details) between the derived intensity image in (2) and the 15-m panchromatic image.These steps are described below.
Step 1: All the 30-m Landsat-8 bands used in this study (Table 1) are resampled to 15 m resolution to align precisely with the Landsat-8 15-m panchromatic band.Unlike many commercial sensors, the 15-m band is not perfectly nested within the lower resolution 30-m bands [46], rather, the panchromatic pixel grid is shifted by 7.5 m in the row and column directions relative to the 30-m pixel grid (Figure 1).The 30-m bands are resampled to the 15-m panchromatic grid using a method to compensate for this kind of pixel shift [47].First, the number of columns and rows in the Landsat-8 30-m image are doubled by inserting a zero-value pixel every second row and column.Then a 7 × 7 low-pass filter (Table 2 in [47]), equivalent to cubic convolution resampling, is applied.
Step 2: A 15-m intensity image is derived as: * , = 0.4030 * , + 0.5177 * , + 0.0802 * , , where * , is the intensity image reflectance at 15-m pixel location , ; * , , * , , and * , are the Landsat-8 15-m red (band 4), green (band 3) and blue (band 2) reflectance values resampled from the 30-m data (Step 1).The spectral weights in Equation ( 1) were derived considering a large number of laboratory spectra and the Landsat-8 OLI spectral response functions [46].The intensity image is derived using the red, green, and blue bands as they overlap spectrally with the panchromatic band (Table 1).The * symbol denotes that the 15-m data that were resampled from the 30-m data.
Step 3: The downscaling is then undertaken as: where , is the downscaled 15-m reflectance value for a given band (i.e., in this study for any of Landsat-8 bands 2-7 in Table 1) for 15-m pixel location ( , ), * , is the 15-m reflectance resampled from the 30-m data (Step 1), , is the 15-m panchromatic reflectance value, and * , is the 15-m intensity image value (Step 2).
Figure 2 illustrates the LPAD downscaling correctly handling the 7.5 m grid shift between the 15-m panchromatic and 30-m reflective bands (middle) and incorrectly handling the grid shift by simply assuming that the 15-m and 30-m bands nest perfectly (right).Qualitatively the LPAD 15-m data (middle and right) provide clearer field boundaries and more apparent spatial detail than the 30-m data (left).Also evident is that the correct handling of the 7.5 m grid shift between the 15-m panchromatic and 30-m Landsat-8 bands is needed with less apparent distortion around the circular The 30-m bands are resampled to the 15-m panchromatic grid using a method to compensate for this kind of pixel shift [47].First, the number of columns and rows in the Landsat-8 30-m image are doubled by inserting a zero-value pixel every second row and column.Then a 7 × 7 low-pass filter (Table 2 in [47]), equivalent to cubic convolution resampling, is applied.
Step 2: A 15-m intensity image is derived as: where I * (i, j) is the intensity image reflectance at 15-m pixel location (i, j); ρ * red (i, j), ρ * green (i, j), and ρ * blue (i, j) are the Landsat-8 15-m red (band 4), green (band 3) and blue (band 2) reflectance values resampled from the 30-m data (Step 1).The spectral weights in Equation ( 1) were derived considering a large number of laboratory spectra and the Landsat-8 OLI spectral response functions [46].The intensity image is derived using the red, green, and blue bands as they overlap spectrally with the panchromatic band (Table 1).The * symbol denotes that the 15-m data that were resampled from the 30-m data.
Step 3: The downscaling is then undertaken as: Remote Sens. 2017, 9, 755 5 of 18 where ρband (i, j) is the downscaled 15-m reflectance value for a given band (i.e., in this study for any of Landsat-8 bands 2-7 in Table 1) for 15-m pixel location (i, j), ρ * band (i, j) is the 15-m reflectance resampled from the 30-m data (Step 1), ρ pan (i, j) is the 15-m panchromatic reflectance value, and I * (i, j) is the 15-m intensity image value (Step 2).
Figure 2 illustrates the LPAD downscaling correctly handling the 7.5 m grid shift between the 15-m panchromatic and 30-m reflective bands (middle) and incorrectly handling the grid shift by simply assuming that the 15-m and 30-m bands nest perfectly (right).Qualitatively the LPAD 15-m data (middle and right) provide clearer field boundaries and more apparent spatial detail than the 30-m data (left).Also evident is that the correct handling of the 7.5 m grid shift between the 15-m panchromatic and 30-m Landsat-8 bands is needed with less apparent distortion around the circular field edges.

Reprojection and Resampling of the Downscaled 15-m Data into Registration with the Sentinel-2A 20-m Data
Landsat-8 and Sentinel-2A data are misaligned because their tiling schemes are different and because the Landsat-8 geolocation framework contains residual errors that vary among Landsat path/rows [19].The pairs of Sentinel-2A and Landsat-8 images (Table 2) were aligned using an automated feature and area-based least squares matching approach that provides sub-pixel precision matching [18].First the Landsat-8 30-m near infrared (NIR) band was resampled to 20-m resolution and then matched with the Sentinel-2A NIR 20-m band.Then following the recommendations of Yan et al. [18] an affine transformation was used to characterize the spatial relation between the matched sensor data, defined as: where , are the UTM coordinates in the Landsat-8 image, , are the UTM coordinates in the Sentinel-2A 20-m data, and , are the affine transformation coefficients.Each band of the Landsat-8 15-m data, derived as Equation ( 2), was reprojected into registration with the Sentinel-2A 20-m grid using the affine transformation coefficients.In this way the Landsat-8 15-m band data were only resampled once.The indirect approach [48] was used whereby the location of each Sentinel-2A 20-m image pixel was reprojected as Equation (3) into the Landsat-8 image and then the Landsat-8 15-m band values were resampled to 20-m resolution.Bilinear and cubic convolution resamplers were considered and for brevity the results are termed as bilinear-based LPAD and cubic convolution-based LPAD results respectively.

Evaluation Methodology
To provide a control data set, Landsat-8 downscaling from 30 m to 20 m was undertaken using conventional resampling, i.e., without the LPAD approach.Nearest neighbor, bilinear and cubic convolution resamplers were considered.The nearest neighbor resampling method selects the closest pixel and so preserves the input data values but is known to create pixel "blockiness" and locational shifts up to √2 2 ⁄ pixels [49].Bilinear and cubic convolution resamplers smooth the data by interpolating from the 4 and 16 surrounding pixels, respectively [50,51].Under certain assumptions  Landsat-8 and Sentinel-2A data are misaligned because their tiling schemes are different and because the Landsat-8 geolocation framework contains residual errors that vary among Landsat path/rows [19].The pairs of Sentinel-2A and Landsat-8 images (Table 2) were aligned using an automated feature and area-based least squares matching approach that provides sub-pixel precision matching [18].First the Landsat-8 30-m near infrared (NIR) band was resampled to 20-m resolution and then matched with the Sentinel-2A NIR 20-m band.Then following the recommendations of Yan et al. [18] an affine transformation was used to characterize the spatial relation between the matched sensor data, defined as: where x L , y L are the UTM coordinates in the Landsat-8 image, x S , y S are the UTM coordinates in the Sentinel-2A 20-m data, and a i , b i are the affine transformation coefficients.Each band of the Landsat-8 15-m data, derived as Equation ( 2), was reprojected into registration with the Sentinel-2A 20-m grid using the affine transformation coefficients.In this way the Landsat-8 15-m band data were only resampled once.The indirect approach [48] was used whereby the location of each Sentinel-2A 20-m image pixel was reprojected as Equation (3) into the Landsat-8 image and then the Landsat-8 15-m band values were resampled to 20-m resolution.Bilinear and cubic convolution resamplers were considered and for brevity the results are termed as bilinear-based LPAD and cubic convolution-based LPAD results respectively.

Evaluation Methodology
To provide a control data set, Landsat-8 downscaling from 30 m to 20 m was undertaken using conventional resampling, i.e., without the LPAD approach.Nearest neighbor, bilinear and cubic convolution resamplers were considered.The nearest neighbor resampling method selects the closest pixel and so preserves the input data values but is known to create pixel "blockiness" and locational shifts up to √ 2/2 pixels [49].Bilinear and cubic convolution resamplers smooth the data by interpolating from the 4 and 16 surrounding pixels, respectively [50,51].Under certain assumptions cubic convolution is theoretically optimal but usually these assumptions are not met and consequently cubic convolution resampled data may have edge-over and under-shoot issues along high contrast edges [52].
The Landsat-8 data downscaled to 20-m resolution using the LPAD method and the control, i.e., conventionally resampled to 20-m Landsat-8 data, were compared with the corresponding Sentinel-2A 20-m data for each site (Table 2).As the Sentinel-2A red, green, and blue bands are defined at 10 m (Table 1) these bands were degraded to 20 m by simple averaging without considering the Sentinel-2A point spread function as it is currently unknown.Note that, unlike for Landsat-8, four adjacent Sentinel-2A 10-m pixels nest within each 20-m pixel, and so averaging is straightforward.
The LPAD and control Landsat-8 20-m data were compared visually and quantitatively with the Sentinel-2A 20-m data.This was undertaken considering the three true color (red, green, and blue) bands and considering all six spectrally similar bands (Table 1).Our expectation is that because the Landsat-8 15-m panchromatic band has a bandwidth that includes only the Landsat-8 red, green, and blue bands the LPAD approach will work better for the true color than for the other bands.
Quantitative comparison of the data was undertaken using the Q2 n metric that quantifies the average spatial and spectral distortion between any number of bands between two images [53].The Q2 n metric is a single value bounded between 0 and 1 (where higher values indicate greater quality and less distortion) and is an extension of the single band universal image quality index (UQI) [54].The UQI is widely used as, similar to human visual perception, it is sensitive to the structural information differences, which is not the case for conventional image difference statistics.In this study the Q2 n values were calculated with respect to non-overlapping adjacent n × n pixel windows and then averaged to yield a single image value.The same n = 32 pixel window dimensions used to assess downscaling Landsat-8 30 m to 15 m [46] was implemented.

Conventional Resampling (Nearest Neighbor, Bilinear, and Cubic Convolution) Based Downscaling
Figure 3 illustrates the results of Landsat-8 downscaling from 30 m to 20 m using three conventional resampling approaches.The New York true color image subset is illustrated because it contains high spatial frequency features associated with closely packed city blocks, roads, rivers and parks.Compared to the Sentinel-2A 20-m data (a) the conventionally resampled Landsat-8 data (b), (c), and (d), do not contain the same level of spatial detail.In particular, the nearest-neighbor resampled data (b) are inferior with the aforementioned "blockiness" and locational shifts apparent with much of the detail in the areas of road and buildings quite obviously lost.The bilinear (c) and cubic convolution (d) resampled Landsat-8 data appear very similar and both smooth the spatial detail.
The Q2 n values derived by comparing the illustrated 20-m Sentinel-2A image with the nearest neighbor, bilinear, and cubic convolution resampled Landsat-8 20-m images are 0.7142, 0.7520, and 0.7760, respectively.Thus, quantitatively the nearest neighbor and cubic convolution resampled data were the most and least distorted (spectrally and spatially) respectively compared to the Sentinel-2A data.Given the better performance of the bilinear and cubic convolution resamplers, for the remainder of this study they are both considered in the LPAD downscaling, i.e., bilinear-based LPAD and cubic convolution-based LPAD results are presented.The 2 values derived by comparing the illustrated 20-m Sentinel-2A image with the nearest neighbor, bilinear, and cubic convolution resampled Landsat-8 20-m images are 0.7142, 0.7520, and 0.7760, respectively.Thus, quantitatively the nearest neighbor and cubic convolution resampled data were the most and least distorted (spectrally and spatially) respectively compared to the Sentinel-2A data.Given the better performance of the bilinear and cubic convolution resamplers, for the remainder of this study they are both considered in the LPAD downscaling, i.e., bilinear-based LPAD and cubic convolution-based LPAD results are presented.It is evident that the LPAD 20-m results provide more clear object boundaries and finer spatial detail than the conventionally resampled 20-m data.Visually, the LPAD generated 20-m data appear spatially and spectrally almost the same as the Sentinel-2A 20-m reference data.The texture of the New York city blocks and even other fine objects can be unambiguously identified (Figure 4d,f).In contrast, the object boundaries are blurry in the 20-m urban areas resampled from the original 30-m data using bilinear and cubic convolution resamplers (Figure 4c,e) and the fine spatial and textural details of the city blocks are degraded.The LPAD generated Landsat-8 20-m crop fields (Figure 5d,f) have sharper boundaries than the Landsat-8 20-m conventionally resampled band data (Figure 5c,e).Similarly, the LPAD generated Landsat-8 20-m burned areas (Figure 6d,f) boundaries appear more distinct than the bilinear and cubic convolution resampling generated 20-m results (Figure 6c,e).The LPAD generated Landsat-8 20-m landslides (Figure 7d,f) appear identical to the Sentinel-2A 20-m reference data (Figure 7b).The small and elongated landslides can be visually identified in LPAD generated 20-m data but they are blurry and cannot be accurately identified in the bilinear and cubic convolution resampled 20-m data (Figure 7c,e).The LPAD generated 20-m deforested patches (Figure 8d,f) appear almost consistent with the Sentinel-2A 20-m reference data (Figure 8b) and their boundaries are more blurred in the bilinear and cubic convolution resampled data (Figure 8c,e).It is evident that the LPAD 20-m results provide more clear object boundaries and finer spatial detail than the conventionally resampled 20-m data.Visually, the LPAD generated 20-m data appear spatially and spectrally almost the same as the Sentinel-2A 20-m reference data.The texture of the New York city blocks and even other fine objects can be unambiguously identified (Figure 4d,f).In contrast, the object boundaries are blurry in the 20-m urban areas resampled from the original 30-m data using bilinear and cubic convolution resamplers (Figure 4c,e) and the fine spatial and textural details of the city blocks are degraded.The LPAD generated Landsat-8 20-m crop fields (Figure 5d,f) have sharper boundaries than the Landsat-8 20-m conventionally resampled band data (Figure 5c,e).Similarly, the LPAD generated Landsat-8 20-m burned areas (Figure 6d,f) boundaries appear more distinct than the bilinear and cubic convolution resampling generated 20-m results (Figure 6c,e).The LPAD generated Landsat-8 20-m landslides (Figure 7d,f) appear identical to the Sentinel-2A 20-m reference data (Figure 7b).The small and elongated landslides can be visually identified in LPAD generated 20-m data but they are blurry and cannot be accurately identified in the bilinear and cubic convolution resampled 20-m data (Figure 7c,e).The LPAD generated 20-m deforested patches (Figure 8d,f) appear almost consistent with the Sentinel-2A 20-m reference data (Figure 8b) and their boundaries are more blurred in the bilinear and cubic convolution resampled data (Figure 8c,e).In general, objects with pixel-scale small axis dimensions, such as narrow roads, are visually discernable in the LPAD generated Landsat-8 20-m data but not in the bilinear and cubic convolution resampled 20-m data.

LPAD Downscaling
Remote Sens. 2017, 9, 755 8 of 18 In general, objects with pixel-scale small axis dimensions, such as narrow roads, are visually discernable in the LPAD generated Landsat-8 20-m data but not in the bilinear and cubic convolution resampled 20-m data.Table 3 summarizes the results considering just the red, green, and blue bands.For all the sites, except the Congo forest site, the conventional bilinear resampling approach has the most distortion, followed in descending order of distortion by the conventional cubic convolution resampling, cubic convolution-based LPAD, and bilinear-based LPAD approaches.The Q2 n differences between the conventional bilinear and cubic-convolution resampling site results are expected for the reasons discussed above (Section 5.1).However, less immediately expected is the higher performance of the bilinear-based LPAD approach compared to the cubic convolution-based LPAD approach.
The Congo forest site results are anomalous compared to the other sites (Table 3).There is negligible difference between the conventional bilinear and bilinear-based LPAD approaches (both have Q2 n = 0.733 to 3 d.p), the least distortion is for the conventional cubic convolution resampling approach (0.743), and the cubic convolution-based LPAD approach has the most distortion (0.708).Unlike for the other sites, the Congo Sentinel-2A and Landsat-8 data were acquired one day apart (Table 2) and the atmosphere changed during this period, with aerosol from biomass burning apparent in the Sentinel-2A data with much less apparent aerosols the day before in the Landsat-8 data (Figure 8).  4 summarizes the same results as for Table 3 but considering all six corresponding bands and not just the red, green and blue bands.The same pattern of Q2 n values in Table 4 as Table 3 is found.For all the sites, except the Congo forest site, the conventional bilinear resampling approach has the most distortion, followed in descending order of distortion by the conventional cubic convolution resampling, cubic convolution-based LPAD, and then the bilinear-based LPAD approaches.

Discussion
The integration of Landsat-8 and Sentinel-2 data for both research and applications is expected to receive growing attention in the global change research community [3,4].One challenge is the spatial resolution difference between the sensors.Rather than downscaling by conventional resampling, that we showed (Figure 3) will cause either smoothing of the downscaled results (e.g., with convention bilinear or cubic convolution resampling) or introduce locational shifts (e.g., with nearest neighbor resampling or pixel duplication approaches), the Landsat-8 15-m panchromatic band was used to assist the downscaling of the Landsat-8 OLI 30-m bands to 20 m.We note that it is not meaningful to downscale Landsat-8 data to 10 m using the LPAD method as the Landsat-8 panchromatic band is defined at 15 m.
Visual comparisons demonstrated that the LPAD method introduced less spectral and spatial distortion than conventional bilinear and cubic convolution resampling.The bilinear-based LPAD approach performed systematically better than the cubic convolution-based LPAD approach.This was because after downscaling from 30 m to 15 m the Landsat-8 15-m data were resampled to 20 m.This resampling from 15 m to 20 m is more like a pixel aggregation problem than an interpolation problem.In such cases, previous studies have noted the advantage of bilinear over cubic convolution resampling [55,56].
Our expectation was that because the Landsat-8 15-m panchromatic band has a bandwidth that includes only the Landsat-8 red, green and blue bands [46] the LPAD approach would work better for these three OLI bands than for the other bands (Table 1).However, the same pattern of the Q2 n values considering just the red, green, and blue bands (Table 3) and considering all six spectral bands (Table 4) was found.Perhaps this reflects the reduced atmospheric scattering that occurs at longer wavelengths [57,58].The results of this study indicate that the LPAD approach can be applied to downscale all of the Landsat-8 30-m OLI bands.
The Democratic Republic of the Congo results were anomalous, with negligible difference between the Q2 n values for the conventional bilinear and bilinear-based LPAD approaches.They were included, not because there is anything particularly different about the Congo scene content compared to the other four sites, but because their results illustrate an important issue.Specifically, atmospheric and perhaps surface differences (due to biomass burning) between the Landsat-8 and Sentinel-2A images meant that differences between the downscaled Landsat-8 20 m and Sentinel-2A 20 m were not due only to the downscaling.The Sentinel-2A image was acquired a day after the Landsat-8 image and contained spatially extensive aerosols that are known to scatter and absorb radiation and that effectively reduce the image contrast [57,58].Atmospheric correction methodologies that use radiative transfer algorithms and atmospheric characterization data are most suitable for large area and/or repeat atmospheric correction [25].Unfortunately, to date, there is no publically available Landsat-8 and Sentinel-2A atmospheric correction approach that uses the same radiative transfer algorithm and atmospheric characterization approach.Evidently, fusion of Landsat-8 and Sentinel-2A data should be undertaken with caution when images are acquired at different times [59] unless the data are acquired with similar surface and atmospheric conditions or are pre-processed very robustly.
In this study we did not consider sensor differences due to different viewing geometry conditions combined with surface reflectance anisotropy [15,16], or sensor spectral band pass differences [60,61], or differences in the interaction of radiation with the atmosphere across high contrast edges ("adjacency effects") that are dependent on the atmospheric contents and the sensor point spread function [62,63].These sensor differences are scene dependent and are likely to cause systematic distortions between the Sentinel-2A and Landsat-8 data.Despite this, with the exception of the Congo results, the bilinear-based LPAD approach provided the best downscaling results with less spatial and spectral distortion than conventional resampling approaches.
In the LPAD method the Landsat-8 30-m band data are first downscaled to 15 m using the simple Brovey method [45] with fixed spectral weights derived by considering a large number of laboratory spectra and the Landsat-8 OLI spectral response functions [46].This is computationally efficient and so is appropriate for global application.Computation efficiency is important for large-area or long-term Landsat-8 and Sentinel-2 studies due to the high satellite data volumes.However, for small data volumes, and/or when computational constraints are not an issue [64], advanced downscaling methods [41] could be used to provide less spatial and spectral distortion, although such improvements are usually subtle and less than the noise introduced by atmospheric and sun-sensor-target geometry effects [46].

Conclusions
A new way to downscale Landsat-8 30-m OLI data to Sentinel-2 20-m resolution was presented.Rather than downscale by conventional resampling, that will cause either smoothing of the downscaled results or introduce pixel-level locational shifts, the Landsat-8 15-m panchromatic band is used to assist the Landsat-8 30-m downscaling.The Landsat panchromatic-assisted downscaling (LPAD) approach was tested over five study areas with different land covers and containing anthropogenic and natural features with straight and irregular boundaries.Near-contemporaneous Sentinel-2A 20-m data were used as reference data for qualitative and quantitative comparison.Results demonstrated that (1) the LPAD method introduced less spectral and spatial distortion and provided visually more coherent data than traditional bilinear and cubic convolution resampling approaches; (2) the LPAD approach can be effectively applied to downscaling the Landsat-8 30-m OLI bands; (3) the bilinear-based LPAD approach provided the least spectral and spatial distortion; and (4) Landsat-8 and Sentinel-2A data fusion should be undertaken with caution when images are acquired at different times under different atmospheric conditions.The LPAD source code is made available at GitHub for public use and evaluation at https://github.com/jasonleepolyu/Landsat-Panchromatic-Assisted-Downscaling-LPAD-/tree/master.

Figure 1 .
Figure 1.Illustration of the 7.5 m row and column shifts that occur between the 15-m panchromatic (gray) and the 30-m reflective band (red) Landsat-8 pixels.

Figure 1 .
Figure 1.Illustration of the 7.5 m row and column shifts that occur between the 15-m panchromatic (gray) and the 30-m reflective band (red) Landsat-8 pixels.

3. 2 .
Reprojection and Resampling of the Downscaled 15-m Data into Registration with the Sentinel-2A 20-m Data

Figures 4 -
Figures 4-8 illustrate, for the five sites (Table 2), the Sentinel-2A 20-m true color images (a) and 256 × 256 20-m pixel Sentinel-2A sub-sets selected over regions of high image spatial contrast (b).In addition, for each sub-set, the conventional Landsat-8 bilinear (c) and cubic convolution (e) resampled 20-m data, and the bilinear-based LPAD (d) and cubic convolution-based LPAD (f) 20-m results are illustrated.It is evident that the LPAD 20-m results provide more clear object boundaries and finer spatial detail than the conventionally resampled 20-m data.Visually, the LPAD generated 20-m data appear spatially and spectrally almost the same as the Sentinel-2A 20-m reference data.The texture of the New York city blocks and even other fine objects can be unambiguously identified (Figure4d,f).In contrast, the object boundaries are blurry in the 20-m urban areas resampled from the original 30-m data using bilinear and cubic convolution resamplers (Figure4c,e) and the fine spatial and textural details of the city blocks are degraded.The LPAD generated Landsat-8 20-m crop fields (Figure5d,f) have sharper boundaries than the Landsat-8 20-m conventionally resampled band data (Figure5c,e).Similarly, the LPAD generated Landsat-8 20-m burned areas (Figure6d,f) boundaries appear more distinct than the bilinear and cubic convolution resampling generated 20-m results (Figure6c,e).The LPAD generated Landsat-8 20-m landslides (Figure7d,f) appear identical to the Sentinel-2A 20-m reference data (Figure7b).The small and elongated landslides can be visually identified in LPAD generated 20-m data but they are blurry and cannot be accurately identified in the bilinear and cubic convolution resampled 20-m data (Figure7c,e).The LPAD generated 20-m deforested patches (Figure8d,f) appear almost consistent with the Sentinel-2A 20-m reference data (Figure8b) and their boundaries are more blurred in the bilinear and cubic convolution resampled data (Figure8c,e).

Figures 4 -
Figures 4-8 illustrate, for the five sites (Table 2), the Sentinel-2A 20-m true color images (a) and 256 × 256 20-m pixel Sentinel-2A sub-sets selected over regions of high image spatial contrast (b).In addition, for each sub-set, the conventional Landsat-8 bilinear (c) and cubic convolution (e) resampled 20-m data, and the bilinear-based LPAD (d) and cubic convolution-based LPAD (f) 20-m results are illustrated.It is evident that the LPAD 20-m results provide more clear object boundaries and finer spatial detail than the conventionally resampled 20-m data.Visually, the LPAD generated 20-m data appear spatially and spectrally almost the same as the Sentinel-2A 20-m reference data.The texture of the New York city blocks and even other fine objects can be unambiguously identified (Figure4d,f).In contrast, the object boundaries are blurry in the 20-m urban areas resampled from the original 30-m data using bilinear and cubic convolution resamplers (Figure4c,e) and the fine spatial and textural details of the city blocks are degraded.The LPAD generated Landsat-8 20-m crop fields (Figure5d,f) have sharper boundaries than the Landsat-8 20-m conventionally resampled band data (Figure5c,e).Similarly, the LPAD generated Landsat-8 20-m burned areas (Figure6d,f) boundaries appear more distinct than the bilinear and cubic convolution resampling generated 20-m results (Figure6c,e).The LPAD generated Landsat-8 20-m landslides (Figure7d,f) appear identical to the Sentinel-2A 20-m reference