The world’s leading publication for data science, AI, and ML professionals.

Downscaling a Satellite Thermal Image from 1000 m to 10 m (Python)

Thermal Sharpening of Sentinel-3 Images: From 1 Km to 10 m Using Python in Google Colab

Sentinel-3 thermal image downscaled from 1000 m to 10 m, visualized by the author.
Sentinel-3 thermal image downscaled from 1000 m to 10 m, visualized by the author.

Table of Contents

  1. πŸŒ… Introduction
  2. πŸ’Ύ Downloading Sentinel-3 (1000 m) and Sentinel-2 images (10 m)
  3. βš™οΈ Sentinel-3 Image Processing
  4. 🌑 ️ Temperature-NDVI Space
  5. πŸ“ Sharpening the Thermal Image (1000 m to 10 m)
  6. πŸ—Ί ️ Visualization of the Sharpened Thermal Image
  7. πŸ“„ Conclusion
  8. πŸ“š References

πŸŒ… Introduction

Downscaling the thermal imagery captured by satellites has been extensively studied due to the tradeoff between the spatial and temporal resolution of satellites that provide thermal images. For example, the revisit cycle of Landsat-8 is 16 days, with an original thermal resolution of 100 meters. In contrast, Sentinel-3 can provide daily thermal images, but at a spatial resolution of 1000 meters.

The trade-off between spatial and temporal resolution, Image by the author
The trade-off between spatial and temporal resolution, Image by the author

One approach to address the coarse resolution of thermal images could be launching more satellites equipped with thermal sensors, such as NASA’s Landsat-9, launched in September 2021. In this case, the temporal resolution for both Landsat-8 and Landsat-9 is 8 days (instead of 16 days with one satellite), assuming clear skies.

However, as you can guess, this approach requires a multimillion-dollar investment and several years of effort. Instead, researchers have focused on statistical methods, correlating the visible/near-infrared (VNIR) bands from satellites with higher spatial resolution (but lower temporal resolution) to thermal images from satellites with lower spatial resolution (but higher temporal resolution). For example, studies have shown that the Normalized Difference Vegetation Index (NDVI) calculated from VNIR bands of Sentinel-2 (10m, every 5 days) can be inversely correlated with thermal images from Sentinel-3 (1000m, daily).

But how can we use this relationship to sharpen the thermal images?

The short answer is that if the correlation between the NDVI of Sentinel-2 and the thermal image of Sentinel-3 is strong enough, we could apply that equation to the NDVI at a 10-meter resolution to generate a thermal image at 10-meter resolution as well.

If you want to learn more about satellite bands and the sensor spectrum, please refer to this:

Satellites Can See Invisible Lava Flows and Active Wildfires, But How? (Python)

In this post, we will download a low spatial resolution thermal image from Sentinel-3 and a high spatial resolution VNIR image from Sentinel-2. These two images are captured simultaneously by the respective satellites. Next, we will calculate the NDVI using VNIR bands ((NIR-Red)/(NIR+Red)), upscale it to 1000m, and explore the correlation between NDVI and the thermal band (both at 1000m resolution). Finally, we will use this correlation to generate the temperature map at 10m resolution.

If you enjoyed the beginning, there’s more to discover as you keep reading.


πŸ’Ύ Downloading Sentinel-3 (1000 m) and Sentinel-2 images (10 m)

I’ve already written three posts on downloading Sentinel-2 images in R and Python, as well as downloading Sentinel-3 images in Python. I don’t want to replicate those steps here; instead, I refer you to these posts:

Downloading Sentinel-2 images in R:

How to Download Sentinel-2 Images in R (Updated January 2024)

Downloading Sentinel-2 images in Python:

Downloading Sentinel-2 Imagery in Python with Google Colab (Updated Nov 2023)

Downloading Sentinel-3 images in Python:

Download and Visualize Land Surface Temperature and NDVI from Sentinel-3 Imagery using Python and…

If you don’t want to download images with writing codes, check this post:

How to Download Satellite Images Without Writing a Single Line of Code

The crucial step in applying the statistical method to downscale thermal images is to find clear images captured by satellites simultaneously. You can filter the metadata based on dates, cloud cover, and your area of interest (AOI) before downloading images. If you want to learn more about filtering and working with Sentinel metadata (AOI, cloud cover, etc.), read this:

How Many Cloud-Free Sentinel-2 Images Have Been Captured Near Your Location?

In this post, my AOI is located in California and I found clear images captured by Sentinel-2 and Sentinel-3 on June 19, 2023. Feel free to search for other locations or dates. But if you want to work with the images I downloaded, here is the information:

Sentinel-2: _S2B_MSIL2A_20230620T183919_N0509_R070_T10SFG20230620T224951

satellite = "SENTINEL-2"

level = "S2MSI2A"

AOI = "POLYGON ((-121.0616 37.6391, -120.966 37.6391, -120.966 37.6987, -121.0616 37.6987, -121.0616 37.6391))"

_start_date = "2023–06–19" ; enddate = "2023–06–21"

Sentinel-3: _S3A_SL_2_LST____20230620T181455_20230620T181755_20230620T201625_0179_100_141_2340_PS1_O_NR004

satellite = "SENTINEL-3"

level= "LST"

AOI = "POLYGON ((-121.0616 37.6391, -120.966 37.6391, -120.966 37.6987, -121.0616 37.6987, -121.0616 37.6391))"

_start_date = "2023–06–19" ; enddate = "2023–06–21"

After downloading the NIR and red bands from Sentinel-2 at 10m (the bands that we need to calculate NDVI) and the thermal image from Sentinel-3, you should have these three files in your directory:

Snapshot of the directory, Image by the author
Snapshot of the directory, Image by the author

βš™οΈ Sentinel-3 Image Processing

Sentinel-3 images cover a much larger scene than Sentinel-2. Since we need average NDVI values for each pixel of Sentinel-3, we need to clip the Sentinel-3 image based on the extent of the Sentinel-2 image. To do this, the first step is to reproject the Sentinel-3 image to the same projection as the Sentinel-2 images. This can be accomplished by:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from pyproj import Transformer

input_raster = 'Sentinel-3_L2_LST_reproj.tif'
output_raster = 'Sentinel-3_L2_LST_reproj_32610.tif'
dst_crs = 'EPSG:32610'

# Read the input raster file
with rasterio.open(input_raster) as src:
    # Get the transform, width, and height for the destination CRS
    transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)

    # Set up the destination
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height,
        'dtype': np.float32,
    })

    # Create the destination and write the reprojected data
    with rasterio.open(output_raster, 'w', **kwargs) as dst:
        # Perform the reprojection
        for i in range(1, src.count + 1):
            reproject(
                source=src.read(1).astype(np.float32)* src.scales[0] + src.offsets[0],
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.bilinear)

In this script, we read the Sentinel-3 raster image that we downloaded from the previous step, calculated the transformation parameters for reprojection to the specified CRS, and then exported the new raster file with the reprojected data. After these steps, you should have these four files in your directory:

Snapshot of the directory, Image by the author
Snapshot of the directory, Image by the author

Now that the Sentinel-3 thermal image has been re-projected into the Sentinel-2 coordinate system, we can proceed to clip the thermal image based on the extent of the Sentinel-2 data using the following code:

import rasterio
import numpy as np

# Open the two raster files
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as small_raster:
    with rasterio.open('Sentinel-3_L2_LST_reproj_32610.tif') as big_raster:

        # Get the extent of the smaller raster
        min_x, min_y, max_x, max_y = small_raster.bounds

        # Read the data from the bigger raster within the extent of the smaller raster
        window = rasterio.windows.from_bounds(min_x, min_y, max_x, max_y, big_raster.transform)
        data = big_raster.read(window=window)

        # Update the metadata of the bigger raster to match the extent of the smaller raster
        clipped_meta = big_raster.meta.copy()
        clipped_meta.update({
            'height': window.height,
            'width': window.width,
            'transform': rasterio.windows.transform(window, big_raster.transform),
            'dtype': data.dtype
        })

        # Write the clipped data 
        with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif', 'w', **clipped_meta) as clipped_raster:
            clipped_raster.write(data)

In this script, we read two raster files (Sentinle-3 and Sentinel-2 images), extract the extent of the smaller raster (Sentinel-2), read the corresponding data from the larger raster (Sentinle-3 thermal image), update the metadata of the clipped Sentinle-3 to match the Sentinel-2 image, and then write the clipped thermal image to a new TIFF file.

Snapshot of the directory, Image by the author
Snapshot of the directory, Image by the author

Let’s plot NDVI and the clipped temperature map, side by side:

import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show

# File paths
red_path = '/content/T10SFG_20230620T183919_B04_10m.jp2'
nir_path = '/content/T10SFG_20230620T183919_B08_10m.jp2'

clipped_temperature_path = '/content/Sentinel-3_L2_LST_reproj_32610_clipped.tif'

# Read raster data
with rasterio.open(red_path) as red_src:
  red = red_src.read(1)

with rasterio.open(nir_path) as nir_src:
  nir = nir_src.read(1)

with rasterio.open(clipped_temperature_path) as clipped_temp_ds:
  clipped_temperature = clipped_temp_ds.read(1)

# Calculate NDVI
ndvi = (nir - red) / (nir + red)

# Plot NDVI and temperature side by side
fig, (ax1,ax2) = plt.subplots(ncols=2, figsize=(12, 6))

# Plot NDVI
im1 = ax1.imshow(ndvi, cmap=ndvi_cmap, vmin=0, vmax=0.6)
ax1.set_title('NDVI', fontweight='bold', fontsize=14)
fig.colorbar(im1, ax=ax1, shrink=0.5)

# Plot clipped temperature
im2= ax2.imshow(clipped_temperature, cmap=ndvi_cmap.reversed(), vmin=300, vmax=315)
ax2.set_title('Clipped Temperature', fontweight='bold', fontsize=14)
fig.colorbar(im2, ax=ax2, shrink=0.5)
plt.show()
Sentinel-2 NDVI at 10-meter resolution (left) and clipped Sentinel-3 thermal image at 1000 m resolution (right), visualized by the author
Sentinel-2 NDVI at 10-meter resolution (left) and clipped Sentinel-3 thermal image at 1000 m resolution (right), visualized by the author

As you can see, we were able to successfully clip the Sentinel-3 thermal image based on the extent of the Sentinel-2 map (NDVI in this case). However, the Sentinel-2 image is also truncated, and if we run the zonal statistics to obtain the average NDVI values for thermal pixels, we will get many NaN values:

Sentinel-2 NDVI in 10m (background), Sentinel-3 pixels (grids, 1000m), Image by the author
Sentinel-2 NDVI in 10m (background), Sentinel-3 pixels (grids, 1000m), Image by the author

We will fix this in the next step by excluding the NaN values from the dataframe.


🌑 ️ Temperature-NDVI Space

With two clear images, one thermal image from Sentinel-3 and one VNIR from Sentinel-2, at the same projection and extent, we can execute zonal statistics to obtain an average of NDVI values for each temperature pixel. Basically, through the zonal statistics approach, we aggregate the NDVI map from 10 m to 1000 m, enabling the plotting of NDVI values against thermal values in a temperature-NDVI space.

As mentioned in the introduction, thermal values should be inversely correlated with NDVI values. Higher NDVI indicates a greater fractional cover of vegetation, corresponding to colder pixels, while lower NDVI indicates less vegetation or bare soil, corresponding to warmer pixels. Let’s do the zonal statistics to explore the space between thermal and NDVI values for our AOI:

import rasterio
import rasterio.features
import rasterio.mask
import pandas as pd
import geopandas as gpd
import rasterstats

import rasterio
from rasterio.features import shapes
mask = None

# Open the input raster
with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif') as src:
    # Read the raster band
    image = src.read(1).astype(np.float32)* src.scales[0] + src.offsets[0]
    results = (
        {'properties': {'Temperature': v}, 'geometry': s}
        for i, (s, v)
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
    geoms = list(results)
    gpd_polygonized_raster  = gpd.GeoDataFrame.from_features(geoms)

# Open the rasters
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as nir_src:
  with rasterio.open('T10SFG_20230620T183919_B04_10m.jp2') as red_src:

        # Read the data as float32
        nir = nir_src.read(1).astype(np.float32)* nir_src.scales[0] + nir_src.offsets[0]
        red = red_src.read(1).astype(np.float32)* red_src.scales[0] + red_src.offsets[0]

        # Calculate NDVI
        ndvi = (nir - red) / (nir + red)

        # Calculate the zonal statistics for each polygon
        stats = rasterstats.zonal_stats(gpd_polygonized_raster.geometry, ndvi, affine=nir_src.transform, stats='mean')

        # Add the mean value of NDVI to the dataframe
        gpd_polygonized_raster['NDVI'] = [s['mean'] for s in stats]

# Save the polygon layer to a shapefile
gpd_polygonized_raster.to_file('output_polygons.shp')

# Create a pandas dataframe from the geodataframe
stats_df = pd.DataFrame(gpd_polygonized_raster.drop(columns='geometry'))

# Print the dataframe
print(stats_df)

In this script, we polygonize a thermal raster, calculate NDVI from Sentinel-2 imagery, and compute zonal statistics (mean NDVI) for each temperature pixel.

Sentinel-3 thermal values versus Sentinel-2 NDVI values
Sentinel-3 thermal values versus Sentinel-2 NDVI values

As shown in the table and discussed previously, there are NaN values for some pixels in which we have thermal data from Sentinel-3, due to the truncation of Sentinel-2 images. We will omit the rows that have NaN values from Sentinel-2:

# Drop rows with NaN values 
df_clean = stats_df.dropna(subset=['NDVI'])
df_clean

the output will be:

Sentinel-3 thermal values versus Sentinel-2 NDVI values after removing nan values
Sentinel-3 thermal values versus Sentinel-2 NDVI values after removing nan values

With this clean dataframe, we can plot the temperature-NDVI space by:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import linregress

# Create a scatter plot 
plt.figure(figsize=(8, 6))
sns.scatterplot(x='NDVI', y='Temperature', data=df_clean, palette='coolwarm')

# Set plot title and axis labels
plt.title('1-1 Plot of NDVI vs Temperature', fontsize=16, fontweight='bold')
plt.xlabel('NDVI', fontsize=14)
plt.ylabel('Temperature', fontsize=14)

# Show the plot
plt.show()

The output will be:

The scatter plot of Sentinel-3 temperature vs. Sentinel-2 NDVI, image by the author
The scatter plot of Sentinel-3 temperature vs. Sentinel-2 NDVI, image by the author

We can observe the inverse correlation between NDVI and temperature in this plot, but some points appear to be outliers. This could be attributed to pixels in the images that represent features other than vegetation and bare soil, such as open water, urban areas, or simply pixels with low quality. To remove those pixels, we will only keep NDVI values between 0.1 and 0.6 and temperature values between 300 Kelvin and 330 Kelvin, then update the plot accordingly:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import linregress

# Filter the DataFrame based on the specified conditions
filtered_df = df_clean[(df_clean['NDVI'] >= 0.1) &amp; (df_clean['NDVI'] <= 0.6) &amp;
                 (df_clean['Temperature'] >= 300) &amp; (df_clean['Temperature'] <= 330)]

# Create a scatter plot of Temperature versus NDVI
plt.figure(figsize=(8, 6))
sns.scatterplot(x='NDVI', y='Temperature', data=filtered_df, palette='coolwarm')

# Fit a linear regression model
slope, intercept = np.polyfit(filtered_df['NDVI'], filtered_df['Temperature'], 1)

# Plot the fitted line
x_line = np.linspace(min(filtered_df['NDVI']), max(filtered_df['NDVI']), 100)
y_line = slope * x_line + intercept
plt.plot(x_line, y_line, 'r', label='Fitted line')

# Create the equation string
equation_str = f'y = {slope:.2f}x + {intercept:.2f}'

# Display the equation on the plot
plt.text(min(filtered_df['NDVI']), max(filtered_df['Temperature']), equation_str, fontsize=12, color='red')

# Set plot title and axis labels
plt.title('1-1 Plot of Temperature vs NDVI', fontsize=16, fontweight='bold')
plt.xlabel('Temperature', fontsize=14)
plt.ylabel('NDVI', fontsize=14)

# Show the plot
plt.show()

The output will be:

The scatter plot of Sentinel-3 temperature vs. Sentinel-2 NDVI, image by the author
The scatter plot of Sentinel-3 temperature vs. Sentinel-2 NDVI, image by the author

Now, the inverse correlation is more visible, along with the equation that illustrates the relationship between NDVI values and temperature.


πŸ“ Sharpening the Thermal Image (1000 m to 10 m)

In this step, we will assume that the equation we have found between aggregated NDVI and temperature at 1000 m can also be valid for a 10-meter resolution. By applying the equation to the original NDVI map, we can estimate the temperature at a 10-meter resolution.

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# Open the NIR and RED files
with rasterio.open('T10SFG_20230620T183919_B08_10m.jp2') as src:
    nir = src.read(1)
    meta = src.meta

with rasterio.open('T10SFG_20230620T183919_B04_10m.jp2') as src:
    red = src.read(1)

# Calculate NDVI
ndvi = (nir - red) / (nir + red)

# Estimate temperature using NDVI
temp = -21.85 * ndvi + 314.9

# Create a color ramp for NDVI
ndvi_cmap = plt.cm.RdYlGn

# Plot NDVI and temperature side by side
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6))

# Plot NDVI
im1 = ax1.imshow(ndvi, cmap=ndvi_cmap, vmin=0, vmax=0.6)
ax1.set_title('NDVI', fontweight='bold', fontsize=14)
fig.colorbar(im1, ax=ax1, shrink=0.7)

# Plot temperature
im2 = ax2.imshow(temp, cmap=ndvi_cmap.reversed(), vmin=300, vmax=315)
ax2.set_title('Temperature', fontweight='bold', fontsize=14)
fig.colorbar(im2, ax=ax2, shrink=0.7)

plt.show()

Here we read NIR and red bands, calculate (NDVI) and estimate temperature values based on the derived equation. Next, we visualize the results side-by-side plots for NDVI and temperature. The maps are:

Sentinel-2 NDVI (10 m) is on the left, and downscaled Sentinel-3 temperature (10 m) is on the right. Visualized by the author.
Sentinel-2 NDVI (10 m) is on the left, and downscaled Sentinel-3 temperature (10 m) is on the right. Visualized by the author.

πŸ—Ί ️ Visualization of Sharpened Thermal Image

In this section, we will focus further on the visualization aspect. We plan to zoom in on the central regions of the images, clip the images around the center, and present NDVI maps from Sentinel-2, the original Sentinel-3 thermal image (at 1000 m resolution), and the sharpened Sentinel-3 image (at 10 m) based on our analysis. These visualizations will be presented side by side for comparison and a more detailed examination:

# Plot NDVI and temperature side by side
fig, axs = plt.subplots(ncols=3, figsize=(15, 5))

# Plot NDVI
vmin, vmax = 0, 0.6
ndvi_subset = ndvi[int(0.75 * ndvi.shape[0]):, int(0.75 * ndvi.shape[1]):]
im1 = axs[0].imshow(ndvi_subset, cmap=ndvi_cmap, vmin=vmin, vmax=vmax)
axs[0].set_title('NDVI', fontweight='bold', fontsize=14)
axs[0].set_xticks([])
axs[0].set_yticks([])
fig.colorbar(im1, ax=axs[0], shrink=0.7)

# Plot temperature
with rasterio.open('Sentinel-3_L2_LST_reproj_32610_clipped.tif') as src:
    original_temp = src.read(1)

vmin, vmax = 300, 315
temp_subset = original_temp[int(0.75 * original_temp.shape[0]):, int(0.75 * original_temp.shape[1]):]
im3 = axs[1].imshow(temp_subset, cmap=ndvi_cmap.reversed(), vmin=vmin, vmax=vmax)
axs[1].set_title('Temperature', fontweight='bold', fontsize=14)
axs[1].set_xticks([])
axs[1].set_yticks([])
fig.colorbar(im3, ax=axs[1], shrink=0.7)

# Plot temperature
vmin, vmax = 300, 315
temp_subset = temp[int(0.75 * temp.shape[0]):, int(0.75 * temp.shape[1]):]
im2 = axs[2].imshow(temp_subset, cmap=ndvi_cmap.reversed(), vmin=vmin, vmax=vmax)
axs[2].set_title('Temperature', fontweight='bold', fontsize=14)
axs[2].set_xticks([])
axs[2].set_yticks([])
fig.colorbar(im2, ax=axs[2], shrink=0.7)

# Adjust the spacing between the subplots
fig.subplots_adjust(wspace=0.2)

plt.show()

The maps are:

Sentinel-2 NDVI (10 m) is on the left, Sentinel-3 temperature (1000 m) is in the middle, and downscaled Sentinel-3 temperature (10 m) is on the right, visualized by the author.
Sentinel-2 NDVI (10 m) is on the left, Sentinel-3 temperature (1000 m) is in the middle, and downscaled Sentinel-3 temperature (10 m) is on the right, visualized by the author.

πŸ“„ Conclusion

The straightforward correlation between Visible and Near-Infrared (VNIR) bands and thermal imagery has proven to be a useful method for enhancing the resolution of thermal images. This technique has practical applications in estimating temperature at a higher spatial resolution, especially in the absence of satellites with proper spatial resolution. This downscaling method serves as a handy tool when a high-resolution thermal map is needed and provides details on small-scale temperature variations. Looking forward, the launch of more satellites with advanced thermal sensors will give us more frequent thermal images at higher resolution. Until then, this method remains a cost-effective option for achieving higher-resolution thermal imagery.


πŸ“š References

Copernicus Sentinel data [2024] for Sentinel data

Copernicus Service information [2024] for Copernicus Service Information.

Agam, N., Kustas, W. P., Anderson, M. C., Li, F., Neale, C. M. U. (2007). A vegetation index-based technique for spatial sharpening of thermal imagery. Remote Sensing of Environment, 107(4), 545–558. ISSN 0034–4257.

Gao, F., Kustas, W. P., Anderson, M. C. (2012). A Data Mining Approach for Sharpening Thermal Satellite Imagery over Land. Remote Sensing, 4, 3287–3319.

Huryna, H., Cohen, Y., Karnieli, A., Panov, N., Kustas, W. P., Agam, N. (2019). Evaluation of TsHARP Utility for Thermal Sharpening of Sentinel-3 Satellite Images Using Sentinel-2 Visual Imagery. Remote Sensing, 11, 2304.


πŸ“± Connect with me on other platforms for more engaging content! LinkedIn, ResearchGate, Github, and Twitter.

Here are the relevant posts available through these links:

Tracking The Great Salt Lake’s Shrinkage Using Satellite Images (Python)

Satellites Can See Invisible Lava Flows and Active Wildfires, But How? (Python)

Watching Storms from Space: A Python Script for Creating an Amazing View


Related Articles