Visualising the RGB channels of Satellite Images with Python | by Conor O’Sullivan | Apr, 2023
How to deal with multiple spectral bands, large pixel values and skewed RGB channels when visualising satellite images
Satellite images contain an abundance of information. The downside is that it is not straightforward to visualise them. The data differs from normal images as they can have:
- Over 12 channels
- Large pixel values
- Skewed pixel values
We will discuss these key considerations. We then factor them into a Python function that will give you flexibility when combining the RGB channels. Specifically, it allows you to adjust the brightness and hue of your images. The code is given and you can find the full project on GitHub.
Our first import is the Geospatial Data Abstraction Library (gdal). This can be useful when working with remote sensing data. We also have more standard Python packages (lines 4–5). Finally, glob is used to handle file paths (line 7).
# Imports
from osgeo import gdalimport numpy as np
import matplotlib.pyplot as plt
import glob
We will be visualising images from the test set Sentinel-2 waters edge dataset (SWED). This contains 98 images of coastlines from 49 locations. We load all the paths to these images below.
#Load paths
paths = glob.glob("../data/SWED/test/images/*")
Number of spectral bands
We load the first image (line 2) and output its shape (line 3). This gives us a value of (12,256,256). In other words, each image consists of 12 channels of 256×256 pixels. If this was a normal image, we would expect dimensions of (256,256,3). That is each pixel has 3 channels — RGB.
#Load first image
img = gdal.Open(paths[0]).ReadAsArray()
img.shape #(12,256,256)
The 12 channels are the sentinel-2 spectral bands. Each band is a measurement of a different wavelength of electromagnetic radiation reflected off the Earth’s surface. This includes infrared light that humans cannot see. When visualising the images, we are only interested in the visible light bands.
For our particular image, the visible light bands are found in positions 3 (red), 2 (green) and 1 (blue). We select these below. We also transpose the array so it has dimensions (256,256,3).
#Get RGB image
rgb = img[[3,2,1]].transpose(1,2,0)
Large pixel values
To capture as much detail as possible, satellite images have a large dynamic range. As a result, pixel values do not fall into the typical range of [0, 255]. For our image, printing out the min and max pixel values gives us 150 and 8,600.
# Pixel range
print(rgb.min(),rgb.max()) #150 8,600
The sentinel-2 images have a maximum reflectance value of 10000. Although pixels can occasionally have values higher than this, we can ignore these large values when visualising the RGB channels. So, we scale images by dividing by 10000 and clipping them between 0 and 1. This ensures all pixel values will be between 0 and 1.
#Scale image
rgb = np.clip(rgb/10000,0,1)
Skewed pixel values
We can now display our satellite image using matplotlib (line 2). You will notice the image has a low brightness/contrast. This is typical when displaying RGB channels.
#Display RGB image
plt.imshow(rgb)
As mentioned, the dynamic range of satellite images is large. This is so they can capture a wider range of pixels — from very bright to very dark. The result is that pixels tend to be skewed towards lower values. You can see this in the histograms below. These give the pixel frequencies for the RGB channels of the above image.
#Display histograms of pixel intesity for each band
fig, axs = plt.subplots(1,3,figsize=(18,5))
fig.patch.set_facecolor('xkcd:white')labels = ['Red','Green','Blue']
for i,ax in enumerate(axs):
ax.hist(rgb[:,:,i].flatten(),bins=100)
ax.set_title(labels[i],size=20,fontweight="bold")
ax.set_xlabel("Pixel Value",size=15)
ax.set_ylabel("Frequency",size =15)
A simple approach to getting better images is to clip the range of pixel values for each channel (line 2). We take only the pixel values from 0 to 0.3 and scale them back to between 0 and 1. In Figure 3, you can see the resulting image is brighter.
# Clip RGB image to 0.3
rgb = np.clip(rgb,0,0.3)/0.3plt.imshow(rgb)
For a quick visualisation, the above is all you will need. If you want more control over the process, you can use the function below. This will scale the image, select the RGB channels and clip each channel using a different cutoff. There is also a display option that will output the RGB histograms with the selected cutoffs.
def visualise_rgb(img,clip=[0.3,0.3,0.3],display=True):
"""Visulaise RGB image with given clip values and return image"""# Scale image
img = np.clip(img/10000,0,1)
# Get RGB channels
rgb = img[[3,2,1]]
#clip rgb values
rgb[0] = np.clip(rgb[0],0,clip[0])/clip[0]
rgb[1] = np.clip(rgb[1],0,clip[1])/clip[1]
rgb[2] = np.clip(rgb[2],0,clip[2])/clip[2]
rgb = rgb.transpose(1,2,0)
if display:
#Display histograms of pixel intesity with given clip values
fig, axs = plt.subplots(1,4,figsize=(22,5))
fig.patch.set_facecolor('xkcd:white')
labels = ['Red','Green','Blue']
for i,ax in enumerate(axs[0:3]):
ax.hist(img[3-i].flatten(),bins=100)
ax.set_title(labels[i],size=20,fontweight="bold")
ax.axvline(clip[i],color="red",linestyle="--")
ax.set_yticks([])
#Display RGB image
axs[3].imshow(rgb)
axs[3].set_title("RGB",size=20,fontweight="bold")
axs[3].set_xticks([])
axs[3].set_yticks([])
return rgb
You can see this function in action below. In this case, we clip each channel using a cutoff of 0.3. Figure 4 shows the histograms with these cutoffs and the resulting RGB image. This can be useful when experimenting with different cutoffs.
img = gdal.Open(paths[0]).ReadAsArray()
rgb = visualise_rgb(img,[0.3,0.3,0.3])
Adjusting brightness
Different images will have different optimal cutoffs. In fact, the cover image of this article was created using 3 different values. The function above makes it easy to adjust them. In Figure 5, you can see how adjusting the cutoffs will change the brightness.
rgb_1 = visualise_rgb(img,[0.15,0.15,0.15],display=False)
rgb_2 = visualise_rgb(img,[0.3,0.3,0.3],display=False)
rgb_3 = visualise_rgb(img,[0.45,0.45,0.45],display=False)
Adjusting colour (hue)
It also gives you more control over each colour channel. As seen in Figure 6, if we reduce the cutoff for a given channel that colour will become more prominent. This can help you fine-tune your visualisations.
rgb_1 = visualise_rgb(img,[0.2,0.3,0.3],display=False)
rgb_2 = visualise_rgb(img,[0.3,0.3,0.2],display=False)
rgb_3 = visualise_rgb(img,[0.3,0.2,0.3],display=False)
When working with remoting sensing data is often important to visualise your images. It will allow you to build intuition around your problem. Perhaps even more important, you can create beautiful images used to explain your work to others. The function above will allow you to create the most eye-catching images possible. Keep in mind, you may need to adapt it to your satellite image dataset.
How to deal with multiple spectral bands, large pixel values and skewed RGB channels when visualising satellite images
Satellite images contain an abundance of information. The downside is that it is not straightforward to visualise them. The data differs from normal images as they can have:
- Over 12 channels
- Large pixel values
- Skewed pixel values
We will discuss these key considerations. We then factor them into a Python function that will give you flexibility when combining the RGB channels. Specifically, it allows you to adjust the brightness and hue of your images. The code is given and you can find the full project on GitHub.
Our first import is the Geospatial Data Abstraction Library (gdal). This can be useful when working with remote sensing data. We also have more standard Python packages (lines 4–5). Finally, glob is used to handle file paths (line 7).
# Imports
from osgeo import gdalimport numpy as np
import matplotlib.pyplot as plt
import glob
We will be visualising images from the test set Sentinel-2 waters edge dataset (SWED). This contains 98 images of coastlines from 49 locations. We load all the paths to these images below.
#Load paths
paths = glob.glob("../data/SWED/test/images/*")
Number of spectral bands
We load the first image (line 2) and output its shape (line 3). This gives us a value of (12,256,256). In other words, each image consists of 12 channels of 256×256 pixels. If this was a normal image, we would expect dimensions of (256,256,3). That is each pixel has 3 channels — RGB.
#Load first image
img = gdal.Open(paths[0]).ReadAsArray()
img.shape #(12,256,256)
The 12 channels are the sentinel-2 spectral bands. Each band is a measurement of a different wavelength of electromagnetic radiation reflected off the Earth’s surface. This includes infrared light that humans cannot see. When visualising the images, we are only interested in the visible light bands.
For our particular image, the visible light bands are found in positions 3 (red), 2 (green) and 1 (blue). We select these below. We also transpose the array so it has dimensions (256,256,3).
#Get RGB image
rgb = img[[3,2,1]].transpose(1,2,0)
Large pixel values
To capture as much detail as possible, satellite images have a large dynamic range. As a result, pixel values do not fall into the typical range of [0, 255]. For our image, printing out the min and max pixel values gives us 150 and 8,600.
# Pixel range
print(rgb.min(),rgb.max()) #150 8,600
The sentinel-2 images have a maximum reflectance value of 10000. Although pixels can occasionally have values higher than this, we can ignore these large values when visualising the RGB channels. So, we scale images by dividing by 10000 and clipping them between 0 and 1. This ensures all pixel values will be between 0 and 1.
#Scale image
rgb = np.clip(rgb/10000,0,1)
Skewed pixel values
We can now display our satellite image using matplotlib (line 2). You will notice the image has a low brightness/contrast. This is typical when displaying RGB channels.
#Display RGB image
plt.imshow(rgb)
As mentioned, the dynamic range of satellite images is large. This is so they can capture a wider range of pixels — from very bright to very dark. The result is that pixels tend to be skewed towards lower values. You can see this in the histograms below. These give the pixel frequencies for the RGB channels of the above image.
#Display histograms of pixel intesity for each band
fig, axs = plt.subplots(1,3,figsize=(18,5))
fig.patch.set_facecolor('xkcd:white')labels = ['Red','Green','Blue']
for i,ax in enumerate(axs):
ax.hist(rgb[:,:,i].flatten(),bins=100)
ax.set_title(labels[i],size=20,fontweight="bold")
ax.set_xlabel("Pixel Value",size=15)
ax.set_ylabel("Frequency",size =15)
A simple approach to getting better images is to clip the range of pixel values for each channel (line 2). We take only the pixel values from 0 to 0.3 and scale them back to between 0 and 1. In Figure 3, you can see the resulting image is brighter.
# Clip RGB image to 0.3
rgb = np.clip(rgb,0,0.3)/0.3plt.imshow(rgb)
For a quick visualisation, the above is all you will need. If you want more control over the process, you can use the function below. This will scale the image, select the RGB channels and clip each channel using a different cutoff. There is also a display option that will output the RGB histograms with the selected cutoffs.
def visualise_rgb(img,clip=[0.3,0.3,0.3],display=True):
"""Visulaise RGB image with given clip values and return image"""# Scale image
img = np.clip(img/10000,0,1)
# Get RGB channels
rgb = img[[3,2,1]]
#clip rgb values
rgb[0] = np.clip(rgb[0],0,clip[0])/clip[0]
rgb[1] = np.clip(rgb[1],0,clip[1])/clip[1]
rgb[2] = np.clip(rgb[2],0,clip[2])/clip[2]
rgb = rgb.transpose(1,2,0)
if display:
#Display histograms of pixel intesity with given clip values
fig, axs = plt.subplots(1,4,figsize=(22,5))
fig.patch.set_facecolor('xkcd:white')
labels = ['Red','Green','Blue']
for i,ax in enumerate(axs[0:3]):
ax.hist(img[3-i].flatten(),bins=100)
ax.set_title(labels[i],size=20,fontweight="bold")
ax.axvline(clip[i],color="red",linestyle="--")
ax.set_yticks([])
#Display RGB image
axs[3].imshow(rgb)
axs[3].set_title("RGB",size=20,fontweight="bold")
axs[3].set_xticks([])
axs[3].set_yticks([])
return rgb
You can see this function in action below. In this case, we clip each channel using a cutoff of 0.3. Figure 4 shows the histograms with these cutoffs and the resulting RGB image. This can be useful when experimenting with different cutoffs.
img = gdal.Open(paths[0]).ReadAsArray()
rgb = visualise_rgb(img,[0.3,0.3,0.3])
Adjusting brightness
Different images will have different optimal cutoffs. In fact, the cover image of this article was created using 3 different values. The function above makes it easy to adjust them. In Figure 5, you can see how adjusting the cutoffs will change the brightness.
rgb_1 = visualise_rgb(img,[0.15,0.15,0.15],display=False)
rgb_2 = visualise_rgb(img,[0.3,0.3,0.3],display=False)
rgb_3 = visualise_rgb(img,[0.45,0.45,0.45],display=False)
Adjusting colour (hue)
It also gives you more control over each colour channel. As seen in Figure 6, if we reduce the cutoff for a given channel that colour will become more prominent. This can help you fine-tune your visualisations.
rgb_1 = visualise_rgb(img,[0.2,0.3,0.3],display=False)
rgb_2 = visualise_rgb(img,[0.3,0.3,0.2],display=False)
rgb_3 = visualise_rgb(img,[0.3,0.2,0.3],display=False)
When working with remoting sensing data is often important to visualise your images. It will allow you to build intuition around your problem. Perhaps even more important, you can create beautiful images used to explain your work to others. The function above will allow you to create the most eye-catching images possible. Keep in mind, you may need to adapt it to your satellite image dataset.
Denial of responsibility! Techno Blender is an automatic aggregator of the all world’s media. In each content, the hyperlink to the primary source is specified. All trademarks belong to their rightful owners, all materials to their authors. If you are the owner of the content and do not want us to publish your materials, please contact us by email – [email protected]. The content will be deleted within 24 hours.