Create your own Timelapse from Satellite Imagery

TL;DR This blog discusses the use of satellites in our life, explains considerations when working with satellite imagery, and presents a script to create a timelapse for any location on Earth. We utilize open data and provide code examples to help you get started to incorporate satellite imagery in your own projects.

Example Timelapse

How satellites make the world go round

Satellites are used for everything, everywhere, all the time. Whether finding the route to the restaurant where to eat, broadcasting the World Cup, or when navigating to your bike-rental. Their varied uses are a recent development. Ever since Laika made her rounds around the world, satellites have been sent to Earths’ orbit. Nowadays, satellites are used in agriculture, weather forecasting, spreading the internet, mining, defense, banking, and more. Measurements from satellites are often made within the electromagnetic spectrum, allowing us to make beautiful time lapses from their data.

What is possible

Open source data from satellites is something special. It is by definition data for ALL of the world, most of which have great resolution. It is widely available, thanks to the help of programs like the Copernicus Programme (ESA) and the LandSat program (NASA). The data comes as measurements in the electromagnetic spectrum (“photo’s”), positional information (“GNSS”, “GPS”) , as well as height measurements (“InSAR”), with each satellite having its own specs. Within this article, we focus on using electromagnetic measurements in the visible spectrum, to create a timelapse.

Resolution

Before starting to code, we need to manage expectations by discussing resolution. Each satellite has a different resolution, which limits what can be done with data from a given satellite. For satellite imagery measuring in the electromagnetic spectrum, there are four different types of resolution:

• Spatial
• Temporal
• Spectral
• Radiometric

Spatial resolution refers to the area contained in each measurement. In the case of our time series, it refers to the “pixel” size — the smallest discernible measurement within the dataset. Temporal resolution is how often the same position is measured by the satellite, how many images are made per unit of time. Spectral resolution deals with the narrowness of the different bands measured in the electromagnetic spectrum. Humans’ spectral resolutions are for example the red, green and blue wavelengths. Radiometric resolution has to do with the information stored within each pixel, which is the number of bits used to store information for each pixel.

It is important to keep in mind what limitations in resolution mean for your project.

Doing it yourself

The imagery we will use in this demonstration is imagery from the Sentinel 2 satellites. The Sentinel 2 mission has 2 satellites orbiting the Earth with a combined revisit time of 5 days. It has great Python API’s to extract data with, and a spatial resolution of 10 meters. The raw data in processed form is available in 10, 20, or 60 meters spatial resolution. We’ll create a timelapse from downloaded images. The following shows code to create a timelapse for a certain location.

To create a simple timelapse we execute the following steps:

1. Download imagery of a certain location
2. Convert imagery jp2 data into images
3. Convert images into video

Define parameters
Before we can start, we need to define what we want, over what timespan. For now, we use Amsterdam as a proxy for the Aurai HQ. To connect with the API, we need to have our account details made here.

Amsterdam is part of two tiles of Sentinel product, hence we select just one for now. We choose the satellites’ orbit where the whole tile contains data. Importantly for selecting the right imagery for a timeseries, we choose a cloud cover percentage to select images without much clouds.

Sentinel 2 comes in two product types: Level 1C and level 2A. Level 2A provides surface reflectances, compared to top-of-atmosphere reflectances for Level 1C. In laymen’s terms, this means that Level 2A contains processed measurements for the Earth’s surface, whereas Level 1C contains less processed data. We thus choose product type Level 2A for both orbiting satellites.

import os
from sentinelsat import SentinelAPI, placename_to_wkt

# Define folders to download products to and to save images in
folder_images_data = os.getcwd()+'/Amsterdam Image Data/'
folder_images_rgb = os.getcwd()+'/Amsterdam RGB Images/'

# Filename to save timeseries to
ts_out = "Timeseries"

# Connection with Sentinel API
username = '<username>'
password = '<password>'
api = SentinelAPI(username, password, 'https://scihub.copernicus.eu/dhus/')

# Define the time interval and location for which we want to get the images
location = 'Amsterdam'
tileID = '31UFU'    
rel_orbit = "R051"  
date_from = "20190101"
date_to = "20221231"

# Define the product options
product_type = ['S2MSI2A','S2MSI2B']       
res = '60m'                     # 10m/20m/60m
bands = ['B04', 'B03', 'B02']   # if you don't want a RGB image, you can play around with these     
place, _ = placename_to_wkt(location)
max_cloudcover = 15

# Define the image format in the download
image_format_raw = 'jp2'
image_format_ts = 'jpg'

Note that we select only online imagery currently. Offline products can be requested through the Sentinelsat API. We chose not to include these here, as this is but a demonstration. To find a suitable tile ID and relative orbit, you should check some images to find the best match to your liking.

Now we query our options and download the data. Note that this can take some time.

# Query the products
products = api.query(place,
                    date=(date_from, date_to),
                    platformname='Sentinel-2',
                    producttype=product_type,
                    cloudcoverpercentage=(0,max_cloudcover))

# Filter more specifically on relative orbit, tile and whether product is online
from collections import OrderedDict

products = OrderedDict({key: values 
                        for key, values in products.items()
                        if tileID in values['title'] 
                        and rel_orbit in values['title']
                        and api.is_online(values['uuid'])})

# Download the products to our created directory
api.download_all(products, directory_path=folder_images_data)

We define a function for later, to select files ending on a certain string from a given directory.

def find_filepath(extension: str, path: str) -> str:
    """ Function to get the correct filepath within a unzipped products 
    :par extension  file ending for which to look
    :par path       starting point of search
    returns         filepath string
    """
    for filename in os.listdir(path):
        if filename.endswith(extension):
            return path+'/'+filename
        elif os.path.isdir(path+'/'+filename):
            return find_filepath(extension, path+'/'+filename)
        else:
            continue

Sentinel products are downloaded as .zip files, each of which has a file size of ~1Gb. We only want to use three bands for our image, so we unzip only these.

We loop over each product in chronological order to create images from said product.

# Create images chronologically
imagenumber = 0
products = [value for value in products.values()][::-1]
for product in products: 
    imagenumber += 1
    f_png = folder_images_rgb+f"image_{imagenumber:003}.{image_format_ts}"
    
    product_title = product['title']
    product_path = folder_images_data+product_title

We extract bands 2, 3, and 4 from our product, the visible spectrum.

from zipfile import ZipFile
from contextlib import closing

with closing(ZipFile(product_path+'.zip', 'r')) as zObject:
  # Get the date of the image for naming later
  _, _, date, _, _, _, _ = product_path.split("_")
  yyyymmdd = date.split('T')[0]
  yyyy, mm, dd = yyyymmdd[:4], yyyymmdd[4:6], yyyymmdd[6:]
  # Extracting specific file bands from in the zip
  for file in zObject.infolist():
      if any([band in file.filename for band in bands]) \
          and file.filename.endswith(image_format_raw)\
          and date in file.filename\
          and res in file.filename:
          zObject.extract(file, path=folder_images_data)

The three bands are combined into a single .tiff file.

import rasterio

# Get the different bands jp2 images
band2 = rasterio.open(find_filepath(f'_B02_{res}.{image_format_raw}', product_path+'.SAFE'))
band3 = rasterio.open(find_filepath(f'_B03_{res}.{image_format_raw}', product_path+'.SAFE'))
band4 = rasterio.open(find_filepath(f'_B04_{res}.{image_format_raw}', product_path+'.SAFE'))

# Create a profile to merge the bands in - based on one of the native profiles
new_profile = band2.profile
new_profile.update({"count": 3})

# Convert bands to .tiff
f_tiff = folder_images_rgb+date+'.tiff'

with rasterio.open(f_tiff, 'w', **new_profile) as rgb_image:
    rgb_image.write(band2.read(1),1)
    rgb_image.write(band3.read(1),2)
    rgb_image.write(band4.read(1),3)

Values need to be rescaled to create RGB images from the created .tiff files. Note that there is more or less sun, dependent on the season or weather conditions. Play around with rescaling until you are comfortable with the results.

import numpy as np
from skimage import exposure

# Open .tiff
img = rasterio.open(f_tiff)

# Rescale the image reflectance using skimage exposure
image = np.array([img.read(3), img.read(2), img.read(1)]).transpose(1,2,0)
p10, p90 = np.percentile(image, (10,90)) # ignore outliers
image = exposure.rescale_intensity(image, in_range=(p10, p90)) / 100000

Create an image using matplotlib.

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

fig = plt.figure()
ax = show(image.transpose(2,0,1), transform=img.transform)
fig.add_axes(ax)
plt.axis('off')
plt.title(f"Aurai HQ at date: {dd}-{mm}-{yyyy}")
plt.savefig(f_png, bbox_inches='tight', dpi=1200)
plt.close(fig)    

Let’s check how that looks.

Example RGB image from Sentinel 2 data

Success! From the images, we create a video using the ffmpeg command. Note that the example .GIFs are shortened for presentation purposes.

os.chdir(folder_images_rgb)
os.system(f'ffmpeg -framerate 1 -pattern_type glob -i \'*.{image_format_ts}\' -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" -pix_fmt yuv420p -r .7 {ts_out}.gif')

Et voila, you have made your very own timelapse from satellite imagery!

The resulting GIF

Conclusion

Satellites are used in many applications, encompassing most of everyday technology and more. A lot of data from satellites is freely available. When working with this data, it is important to take into consideration the limitations for your specific project. To demonstrate a simple use-case for satellite imagery, we’ve presented a simple script to create timelapses from Sentinel 2 data, useable for your project, at any time, anywhere.

Aurai provides custom data solutions that help companies gain insights into their data. We engineer your company’s future through simplifying, organizing and automating data. Your time is maximized by receiving the automated knowledge effortlessly and enacting better processes on a foundation of relevant, reliable, and durable information. Interested in what Aurai can mean for your organisation? Don’t hesitate to contact us!