GEE Sentinel-2 Download

Author

Affiliation

Scott Forrest

Queensland University of Technology, CSIRO

Abstract

Download monthly composite Sentinel-2 imagery using Google Earth Engine (GEE). I’m not sure why there are so many maps rendered.

Load packages

import os
import ee
import geemap
import json
import requests
import folium

from IPython.display import display, HTML, clear_output

Initialise GEE

You will need to authenticate your Google account the first time you run this.

# Initialize Earth Engine
try:
    ee.Initialize()
    print("Earth Engine already initialized")
except Exception as e:
    ee.Authenticate()
    ee.Initialize()
    print("Earth Engine initialized")

To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions:

https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/cloud-platform%20https%3A//www.googleapis.com/auth/drive%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=8Hna6UzUyx7DzA-HmrkOImrFIJUws5RNTL1OMOUSoSI&tc=TnwUh7VPfdW0R1TCqAoa2clRJEZLEbXx9Bmk1GWfPTE&cc=h7xAclOjEUfu_HviuXD29tT42bo0XGKuCIEb0VUWgjU

The authorization workflow will generate a code, which you should paste in the box below.

Successfully saved authorization token.
Earth Engine initialized

Plot the basemap

The first time you run the code below, follow the link to authenticate your Google account. After following the instructions online, you will receive an authorization code, which you can paste back into the input box in your notebook. If you’re using VSCode, this will appear at the top of your window where the search bar is.

Map = geemap.Map(lite_mode=True, zoom=2)
Map.add_basemap("SATELLITE")
# Map

# Clean up existing HTML files
html_files = ['satellite_map.html', 'satellite_map2.html', 'satellite_map_S2.html']
for file in html_files:
    if os.path.exists(file):
        os.remove(file)

# Clear previous outputs before displaying new map
clear_output(wait=True)
Map.to_html(filename='satellite_map.html')
HTML('satellite_map.html')

Import AOI

# Option 1: Read JSON from file
def load_aoi_from_file(json_file_path):
    with open(json_file_path, 'r') as f:
        geojson = json.load(f)
    
    # Create an ee.Geometry from the GeoJSON
    return ee.Geometry(geojson)

# Alternative: Define AOI from coordinates
# aoi = ee.Geometry.Rectangle([-122.5, 37.5, -122.0, 38.0])  # San Francisco area

# Load AOI from file
aoi_name = 'mimal_test'
aoi = load_aoi_from_file(fr'AOIs/{aoi_name}.geojson')

# New basemap for AOI
Map2 = geemap.Map(lite_mode=True, zoom=2)
Map2.add_basemap("SATELLITE")
Map2

# Show the AOI on map
Map2.addLayer(aoi, {}, 'Area of Interest')
# Center the map on the AOI and zoom to it
Map2.centerObject(aoi, zoom=10)  # Adjust zoom level (1-20) as needed
# Map  # Display the map with the AOI

# Clear previous outputs before displaying new map
clear_output(wait=True)
Map2.to_html(filename='satellite_map2.html')
HTML('satellite_map2.html')

Get Sentinel-2 Collection with Cloud Masking

This function creates a cloud mask for Sentinel-2 imagery.

There are other approaches to filter clouds, such as the approach listed here: https://developers.google.com/earth-engine/tutorials/community/sentinel-2-s2cloudless.

def maskS2clouds_SCL(image):
    """
    Function to mask clouds using the SCL band from Sentinel-2 L2A data.
    SCL values:
    1: Saturated or defective
    2: Dark Area Pixels 
    3: Cloud Shadows
    4: Vegetation
    5: Bare Soils
    6: Water
    7: Clouds Low Probability / Unclassified
    8: Clouds Medium Probability
    9: Clouds High Probability
    10: Cirrus
    11: Snow / Ice
    """
    # Get the Scene Classification Layer
    scl = image.select('SCL')

    # Create a water mask to keep water bodies (value 6 is water)
    waterMask = scl.eq(6)
    
    # Keep pixels that are not clouds (values 7, 8, 9, 10) and not shadows (3)
    # Also exclude saturated/defective pixels (1)
    # Values to mask: 1, 3, 7, 8, 9, 10
    cloudShadowMask  = scl.neq(1).And(scl.neq(3)).And(
        scl.neq(8)).And(scl.neq(9)).And(scl.neq(10)) #.And(scl.neq(7)) # removed waterbodies
    
    # Combine masks: keep pixels that are either water OR pass the cloud/shadow mask
    finalMask = waterMask.Or(cloudShadowMask)
    
    # Scale the image to 0-1 range (default is 0-10000 for SR data)
    scaled = image.divide(10000)
    
    # Apply the mask
    return scaled.updateMask(finalMask).copyProperties(image, ['system:time_start'])

Get Monthly Composites of Sentinel-2

Define a function to get monthly composites

We want to create monthly composites of Sentinel-2 imagery. This function will filter the Sentinel-2 image collection by date and area of interest (AOI), apply the cloud mask, and then compute the median of each ‘cloud-free’ pixel for the month.

def get_monthly_composites(start_date, end_date, aoi):
    """Generate monthly S2 composites for the given date range and AOI."""
    # Convert date strings to ee.Date objects
    start = ee.Date(start_date)
    end = ee.Date(end_date)
    
    # Get number of months
    months = end.difference(start, 'month').round().int()

    # First, check if we have any S2 data for this month (without filtering)
    raw_collection = ee.ImageCollection('COPERNICUS/S2_SR') \
        .filterBounds(aoi) \
        .filterDate(start_date, end_date)
    
    raw_count = raw_collection.size().getInfo()
    print(f"  Found {raw_count} raw S2 images for {start_date} to {end_date}")
    
    # Function to get image for a specific month
    def get_monthly_image(month_index):
        # Calculate current month start and end
        current_month_start = start.advance(month_index, 'month')
        current_month_end = current_month_start.advance(1, 'month')
        
        # Format dates for naming
        date_format = current_month_start.format('YYYY-MM')
        
        # Get S2 collection for this month
        s2_collection = ee.ImageCollection('COPERNICUS/S2_SR') \
            .filterBounds(aoi) \
            .filterDate(current_month_start, current_month_end) \
            .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 90)) \
            .map(maskS2clouds_SCL)
        
        # # Check if we have any images left after cloud filtering
        # filtered_count = s2_collection.size().getInfo()
        # print(f"  After cloud filtering: {filtered_count} images remain")
        
        # Create a median composite and clip to AOI
        composite = s2_collection.median().clip(aoi)
        
        # Select RGB bands and scale for visualization
        # Bands: B2 (blue), B3 (green), B4 (red)
        # These get re-ordered to RGB by the image_cutting_support.create_padded_png_S2 function
        rgb = composite.select(['B2', 'B3', 'B4'])
        
        return rgb.set({
            'system:index': date_format,
            'system:time_start': current_month_start.millis()
        })
    
    # Create a list of months
    month_indices = ee.List.sequence(0, months)
    
    # Map the function over the months
    composites = ee.ImageCollection.fromImages(
        month_indices.map(get_monthly_image)
    )
    
    return composites

Define a function to export and download images

If the images are small enough, you can download them directly to your local machine. If they are too large, you can export them to your Google Drive.

def export_monthly_images(composites, folder_path, scale=10):
    """Export monthly images to local files."""
    # Create folder if it doesn't exist
    os.makedirs(folder_path, exist_ok=True)
    
    # Get the list of images
    image_list = composites.toList(composites.size())
    num_images = image_list.size().getInfo()
    
    # Loop through each image and export
    for i in range(num_images):
        image = ee.Image(image_list.get(i))
        date_str = image.get('system:index').getInfo()
        
        print(f"Processing image {i+1}/{num_images} for date {date_str}...")
        
        # Define output path
        geotiff_path = os.path.join(folder_path, f'{aoi_name}_{date_str}.tif')
        
        # Get download URL for the GeoTIFF
        download_url = image.getDownloadURL({
            'scale': scale,
            'crs': 'EPSG:4326',
            'fileFormat': 'GeoTIFF',
            'region': image.geometry()
        })
        
        # Download the file
        response = requests.get(download_url, stream=True)
        response.raise_for_status()  # Raise an exception for HTTP errors
        
        # Save the file
        with open(geotiff_path, 'wb') as f:
            for chunk in response.iter_content(chunk_size=1024*1024):
                f.write(chunk)
        
        print(f"Saved GeoTIFF: {geotiff_path}")

# For larger or higher resolution images, use Google Drive export:
def export_to_drive(composites, drive_folder='Earth_Engine_Exports', aoi_name='S2_', scale=10):
    """Export images to Google Drive for larger areas or higher resolution."""
    # Get the list of images
    image_list = composites.toList(composites.size())
    num_images = image_list.size().getInfo()
    
    # Loop through each image and export
    for i in range(num_images):
        image = ee.Image(image_list.get(i))
        date_str = image.get('system:index').getInfo()
        
        # Create export task
        task = ee.batch.Export.image.toDrive(
            image=image.visualize(min=0, max=0.3, bands=['B4', 'B3', 'B2']),
            description=f'{date_str}_{aoi_name}_S2',
            folder=drive_folder,
            scale=scale,
            region=aoi.getInfo()['coordinates'],
            maxPixels=1e9
        )
        
        # Start the task
        task.start()
        print(f"Started export task for {date_str}")

Download images for a specified date range

We will get an image for June 2024 as an example.

# Define parameters
start_date = '2024-06-01'
end_date = '2024-06-30' # the function will 
output_folder = 'image/sentinel2_monthly_images'

Run the function

The function will print how many suitable images were found for the specified date range and AOI. If no images are found, you may need to adjust your date range or AOI.

# Get monthly composites
composites = get_monthly_composites(start_date, end_date, aoi)

  Found 12 raw S2 images for 2024-06-01 to 2024-06-30

# Visualize a preview
map2 = geemap.Map()
map2.centerObject(aoi, zoom=10)
map2.add_basemap("SATELLITE")
map2.addLayer(aoi, {}, 'Area of Interest')

# Add the first image to the map
first_image = ee.Image(composites.first())
viz_params = {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.5}
# map2.addLayer(first_image, viz_params, 'First Monthly Composite')
# map2

# Get the number of images in the collection
count = composites.size().getInfo()
print(f"Found {count} monthly composites")

# Add each monthly composite to the map with its date as the layer name
if count > 0:
    # Convert to list for iterating
    image_list = composites.toList(count)
    
    # Add each image as a separate layer
    for i in range(count):
        image = ee.Image(image_list.get(i))
        date_str = image.get('system:index').getInfo()  # or 'date' if you set that property
        
        # Add to map - only the first layer will be visible by default
        map2.addLayer(image, viz_params, f'S2 RGB {date_str}', shown=(i==0))
        
    print(f"Added {count} layers to the map. Use the layer control to toggle visibility.")
else:
    print("No composites available to display.")

# Add the layer control (if not already visible)
map2.add_layer_control()

# Display the map
map2

# Clear previous outputs before displaying new map
clear_output(wait=True)
map2.to_html(filename='satellite_map_S2.html')
HTML('satellite_map_S2.html')

print(f'AOI name:   {aoi_name}')

# Export images locally (for small areas)
# export_monthly_images(composites, output_folder) 

# Or export to Google Drive (for larger areas)
export_to_drive(composites, drive_folder='sentinel2_images', aoi_name=aoi_name, scale=10)

AOI name:   mimal_test
Started export task for 2024-06
Started export task for 2024-07