Creating Satellite Image Timelapses
A while ago, I wrapped up the know-how of collecting and preparing satellite imagery Data from the European Space Agency's Sentinel satellites in my article titled Deep Dive into ESA's Sentinel API. Since then, ESA rolled out a major update unseen in years on the Sentinel Hub with updates API access methods. Hence, I briefly review how to get data from the current API. As an additional use-case, I also show how to merge the downloaded satellite image snapshots into animated gifs using pure Python. This article aims to get you started and be onboard, with the possibility to further explore the API depending on specific use cases and data sets.
All images created by the author.
1. Getting started
First, you will need to install the sentinelhub python library, which you can do in a Jupyter notebook by running the following cell:
import sys
!{sys.executable} -m pip install sentinelhub --upgrade Then, you will also need to sign up and create your access tokens. You may read about the how-to here and get your account up and running on your Sentinel Hub Dashboard here. The Dashboard also informs you about your data usage, including the amount of free credits you have and the possibility to upgrade to a premium plan. To that end, you may also be interested in reading about the different types of data layers the Sentinel Hub provides here under the Data tab. Long story short – sign up for the API, log into the Dashboard, hit User Settings, and g for the OAuth client button.
Let's config our notebook:
from sentinelhub import SHConfig
client_secret = 'GEZixbNoqGcFYiFMNSLz74HrzN03rfvi'
client_id = '10e7a03e-5783-413f-a5ba-3e1a452ce742'
config = SHConfig(sh_client_id = client_id, sh_client_secret = client_secret)2. Download satellite images
Here, I rely on sample codes built by the Sentinel team to show you how to get imagery data for various areas. I try to make it interesting by selecting a few more interesting locations, as you will see below.
First, we define the target area using WGS-84 (long-lat) coordinate representation, select the resolution of the image we wish to download (we can select from 10m, 20m, and 60m pixel size), and define the enclosing bounding box. Additionally, we compute the image size of the final query, which is measured in pixels.
Note: the image's dimensions in each direction are limited to 2500 pixels in the free trial, so make sure to adjust your bounding box and resolution accordingly.
# define the bounding box coordinates using long-lat terminology
target_area_coords = (46.16, -16.15, 46.51, -15.58)
# define the image resolution (10m, 20m, 60m)
resolution = 60
# define the bounding box object
from sentinelhub import CRS, BBox, bbox_to_dimensions
target_bbox = BBox(bbox=target_area_coords, crs=CRS.WGS84)
# a quick one-liner to compute the pixel-size of the final image:
target_size = bbox_to_dimensions(target_bbox, resolution=resolution)
print('The image's size at ' + str(resolution) + ' resolution: ', target_size)The output of this cell:
Next, using the API query samples from the Sentinel tutorials, we can call the API and query data fitting into the bounding box given our resolution setting as follows, also specifying the satellite bands we are interested in. Here, we use bands 2, 3, and 4, corresponding to the red, green, and blue color channels of the visible spectral range.
Additionally, by setting mosaicking_order to LEAST_CC, we will query the satellite image specific to the time interval with the lowest cloud coverage value.
# here we define the input query
from sentinelhub import DataCollection, MimeType, MosaickingOrder, SentinelHubRequest
evalscript_true_color = """
//VERSION=3
function setup() {
return {
input: [{
bands: ["B02", "B03", "B04"]
}],
output: {
bands: 3
}
};
}
function evaluatePixel(sample) {
return [sample.B04, sample.B03, sample.B02];
}
"""
request_true_color = SentinelHubRequest(
evalscript=evalscript_true_color,
input_data=[
SentinelHubRequest.input_data(
data_collection=DataCollection.SENTINEL2_L1C,
time_interval=("2023-01-01", "2024-01-01"),
mosaicking_order=MosaickingOrder.LEAST_CC,
)
],
responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
bbox=target_bbox,
size=target_size,
config=config,
)
true_color_imgs = request_true_color.get_data()
image = true_color_imgs[0]
imageThe output of this cell, the numpy array storing the rasterized satellite image data:
Now, plot the result by scaling the numpy array values in a way that gives the image nice contrast and good visibility in the true-color RGB range. One may play around with these parameters to enhance image quality and contrast further both in Python and even on the final exported images.
import matplotlib.pyplot as plt
import numpy as np
f, ax = plt.subplots(1,1,figsize=(10,10))
factor=2.05 / 255
clip_range=(0, 200)
ax.imshow(np.clip(image * factor, *clip_range))The output of this cell:
Now, let's wrap this up in a function and get visuals of a few interesting spots across the globe.
def download_image(target_area_coords, resolution):
# define the bounding box object
target_bbox = BBox(bbox=target_area_coords, crs=CRS.WGS84)
target_size = bbox_to_dimensions(target_bbox, resolution=resolution)
print('The image's size at ' + str(resolution) + ' resolution: ', target_size)
# init the query
evalscript_true_color = """
//VERSION=3
function setup() {
return {
input: [{
bands: ["B02", "B03", "B04"]
}],
output: {
bands: 3
}
};
}
function evaluatePixel(sample) {
return [sample.B04, sample.B03, sample.B02];
}
"""
# get the data
request_true_color = SentinelHubRequest(
evalscript=evalscript_true_color,
input_data=[
SentinelHubRequest.input_data(
data_collection=DataCollection.SENTINEL2_L1C,
time_interval=("2023-01-01", "2024-01-01"),
mosaicking_order=MosaickingOrder.LEAST_CC,
)
],
responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
bbox=target_bbox,
size=target_size,
config=config,
)
true_color_imgs = request_true_color.get_data()
image = true_color_imgs[0]
# plot the data
f, ax = plt.subplots(1,1,figsize=(10,10))
factor=2.0 / 255
clip_range=(0, 200)
ax.imshow(np.clip(image * factor, *clip_range))
ax.axis('off')A Sneak peak to lake Balaton in Hungary:
target_area_coords = (17.1055894908987192,47.15783093995728, 18.252391267583046, 46.60415867989036)
resolution = 60
download_image(target_area_coords, resolution)
The Perito Moreno Glacier in Argentina:
target_area_coords = (-73.15480000852375,-50.376712087300675,-72.92235409436927, -50.52489277557761)
resolution = 10
download_image(target_area_coords, resolution)
A beautiful atoll in the Maldives:
target_area_coords = (72.70900633500187, 3.032187530924483 , 73.12923489370033, 2.6351608765294827)
resolution = 20
download_image(target_area_coords, resolution)
The Bumbunga Lake in South Australia:
target_area_coords = (138.1467833742433, -33.827908049913276,138.23094555171292, -33.95376441094313)
resolution = 10
download_image(target_area_coords, resolution)
The Palm Jumeirah in Dubai:
target_area_coords = (55.076501959959515, 25.176112472405652, 55.1850748041356, 25.0523812683095)
resolution = 10
download_image(target_area_coords, resolution)
The Elephant foot glacier in Greenland:
target_area_coords = (-19.775848942751963,80.935666488578635, -19.142019641734414, 80.86476334882002)
resolution = 10
download_image(target_area_coords, resolution)
3. Create a timelaps animation
In the last part of this tutorial, I first create an output folder and then update the previous data downloader and visualizer function to save the output image.
Also, I set the values of the time_interval interval in a way that I would have one snapshot – the least cloudy, arriving at random times within each month but allowing the best visibility – from each month going back to 2016, meaning that the timelapse will have more than 50 frames total.
import os
foldout = 'frames'
if not os.path.exists(foldout):
os.makedirs(foldout)Then, I update the downloader function with a few smaller things – to take the year and month as inputs and save the downloaded image instead of showing it in the notebook. Additionally, I annotate each frame by year and month.
import calendar
def save_image(target_area_coords, resolution, year, month, foldout):
# get the naming terminology
month_names = [calendar.month_name[i] for i in range(1, 13)]
timestamp = str(year) + ' - ' + month_names[month-1]
title = foldout + '/' + str(year) + ' - ' + str(month) + ' - ' + month_names[month-1]
if not os.path.exists(title + '.png'):
# define the bounding box object
target_bbox = BBox(bbox=target_area_coords, crs=CRS.WGS84)
target_size = bbox_to_dimensions(target_bbox, resolution=resolution)
# init the query
evalscript_true_color = """
//VERSION=3
function setup() {
return {
input: [{
bands: ["B02", "B03", "B04"]
}],
output: {
bands: 3
}
};
}
function evaluatePixel(sample) {
return [sample.B04, sample.B03, sample.B02];
}
"""
# get the data
request_true_color = SentinelHubRequest(
evalscript=evalscript_true_color,
input_data=[
SentinelHubRequest.input_data(
data_collection=DataCollection.SENTINEL2_L1C,
time_interval=(str(year) + "-" + str(month) + "-01", str(year) + "-" + str(month) + "-28"),
mosaicking_order=MosaickingOrder.LEAST_CC,
)
],
responses=[SentinelHubRequest.output_response("default", MimeType.PNG)],
bbox=target_bbox,
size=target_size,
config=config,
)
true_color_imgs = request_true_color.get_data()
image = true_color_imgs[0].astype(np.uint8)
# plot the data
f, ax = plt.subplots(1,1,figsize=(10,10))
factor=2.0 / 255
clip_range=(0, 200)
ax.imshow(np.clip(image * factor, *clip_range))
ax.axis('off')
ax.annotate(timestamp, xy=(0.0125, 0.967), xycoords='axes fraction', ha='left', va='center',
bbox=dict(boxstyle='round', alpha=0.0), fontsize=26, color='w')
plt.savefig(title + '.png', dpi = 200, bbox_inches = 'tight')
plt.tight_layout()
plt.close()Now run this function and get a series of images about the Perito Moreno Glacier:
target_area_coords = (-73.15480000852375,-50.376712087300675,-72.92235409436927, -50.52489277557761)
resolution = 10
for year in range(2016, 2025):
for month in range(1,11):
#for month in range(1,3):
if year == 2024:
if month<4:
save_image(target_area_coords, resolution, year, month, foldout)
else:
save_image(target_area_coords, resolution, year, month, foldout)Finally, let's assemble the frames into an animation:
from PIL import Image
png_files = [f for f in os.listdir(foldout) if f.endswith('.png')]
png_files.sort()
frames = []
for png_file in png_files:
file_path = os.path.join(foldout, png_file)
img = Image.open(file_path)
frames.append(img)
output_gif_path = 'footage_complete.gif'
frames[0].save(
output_gif_path,
save_all=True,
append_images=frames[1:],
duration=200, # Set the duration between frames in milliseconds
loop=0 # Set loop to 0 for an infinite loop, or any positive integer for a finite loop
)
In this footage, we may see that there are a couple of completely cloudy frames – so in practice, it may be worth dropping those frames to have a smoother animation. One may tweak the frame merger algorithm for the same reason.
Finally, let's have a look and do a year-by-year comparative time lap by only taking each year's January. This footage will better illustrate the really saddening magnitude of how the glacier is decreasing year by year in a clearly detectable and visible way. This example also shows a potent example of using satellite imagery for environmental monitoring.
from PIL import Image
png_files = [f for f in os.listdir(foldout) if 'January' in f]
png_files.sort()
frames = []
for png_file in png_files:
file_path = os.path.join(foldout, png_file)
img = Image.open(file_path)
frames.append(img)
output_gif_path = 'footage_januaries.gif'
frames[0].save(
output_gif_path,
save_all=True,
append_images=frames[1:],
duration=200,
loop=0
)
UPDATE
Using the following code block you can also add crossfading between the frames:
from PIL import Image, ImageDraw
png_files = [f for f in os.listdir(foldout) if f.endswith('.png')]
png_files.sort()
frames = []
for png_file in png_files:
file_path = os.path.join(foldout, png_file)
img = Image.open(file_path)
frames.append(img)
# Define the duration for each frame in milliseconds
frame_duration = 100
# Define the number of frames for the crossfade transition
crossfade_frames = 6
# Create a list to store the crossfaded frames
crossfaded_frames = []
# Iterate through pairs of consecutive frames
for i in range(len(frames) - 1):
# Extract current and next frames
current_frame = frames[i]
next_frame = frames[i + 1]
# Create a sequence of crossfaded frames between the current and next frames
for j in range(crossfade_frames + 1):
# Calculate the alpha value for blending
alpha = j / crossfade_frames
# Blend the current and next frames using alpha blending
blended_frame = Image.blend(current_frame, next_frame, alpha)
# Append the blended frame to the list of crossfaded frames
crossfaded_frames.append(blended_frame)
# Add the last frame without crossfading
crossfaded_frames.append(frames[-1])
output_gif_path = 'footage_complete_2.gif'
# Save the GIF with crossfaded frames
crossfaded_frames[0].save(
output_gif_path,
save_all=True,
append_images=crossfaded_frames[1:],
duration=frame_duration, # Set the duration between frames in milliseconds
loop=0 # Set loop to 0 for an infinite loop, or any positive integer for a finite loop
)