Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Images and Image Collections

import ee
import geemap
# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize(project='your-gee-project-id-here')
Output

Accessing a single Image

# Load the global Digital Elevation Model (DEM) from the Earth Engine data catalog
srtm = ee.Image('USGS/SRTMGL1_003')

Clipping an image to a bounding box

# Define the bounds of the area of interest
bounds = ee.Geometry.BBox(120.7, 14.5, 121.1, 15.2)

# Clip the SRTM DEM to the defined bounds
srtm_clipped = srtm.clip(bounds)
# Show the data
m = geemap.Map()

srtm_viz = {
    'min': 0,
    'max': 2000,
    'palette': ['black', 'white']
}

srtm_ph = {
    'min': 0,
    'max': 200,
    'palette': ['black', 'white']
}

m.add_layer(srtm, srtm_viz, name="SRTM")
m.add_layer(srtm_clipped, srtm_ph, name="SRTM Manila")
m

Accessing an ImageCollection

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

Filter

Filtering is used to to get a smaller ImageCollection from a larger one. In Earth Engine, a filter can limit a large ImageCollection using a combination of criterion such as date, location, or other image characteristics, and return a smaller ImageCollection.

📌 REMINDER

  1. You can always chain and combine filters. For example, you can combine a filterDate with a filterBounds to return an ImageCollection with images that both satisfy a given date criteria AND intersect a given geometry.

  2. Chained filter operations are executed from LEFT to RIGHT.

  3. The order of the filters affect the performance/speed of the filtering especially for large and complex processing. In Earth Engine, it is good practice to order the filters with the filterBounds first, followed by metadata filters in order of decreasing specificity.

# Set visualization parameters for RGB display
vis_rgb = {
    'min': 0.0,
    'max': 3000,
    'bands': ['B4', 'B3', 'B2']  # Red, Green, Blue bands for S2
}

# Define the point geometry for pseudo Puerto Princesa center
pp_point = ee.Geometry.Point([118.75, 9.78])

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

# Apply filters step-by-step
loc_filter = s2.filterBounds(pp_point)
# print(loc_filter.size(), 'Location Filter')

date_loc_filter = loc_filter.filterDate('2024-01-01', '2024-12-31')
# print(date_loc_filter.size(), 'Date Location Filter')

cloud_date_loc_filter = date_loc_filter.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
# print(cloud_date_loc_filter.size(), 'Cloud Date Location Filter')

# Chained filter
chained_filter = s2.filterBounds(pp_point) \
                    .filterDate('2024-01-01', '2024-12-31') \
                    .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)

m1 = geemap.Map()
m1.centerObject(pp_point, 10);
m1.add_layer(cloud_date_loc_filter.first(), vis_rgb, 'Step-by-Step Filtered Image')
m1.add_layer(chained_filter.first(), vis_rgb, 'Chain Filtered Image')

m1

Map

When we map a function to an ImageCollection, it applies the operations to each image resulting in a processed ImageCollection as output. For example, in this exercise, we will add an EVI band to the images we filtered above.

# Function to compute EVI and add it as a band to the input image
# EVI = 2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))
def compute_evi(image):
    nir = image.select('B8')   # NIR for Sentinel-2
    red = image.select('B4')   # Red for Sentinel-2
    blue = image.select('B2')  # Blue for Sentinel-2

    evi = image.expression(
        '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
            'NIR': nir,
            'RED': red,
            'BLUE': blue
        }
    )

    return image.addBands(evi.rename('EVI'))

# Function to compute NDVI and add it as a band to the input image
def compute_ndvi(image):
    ndvi = image.normalizedDifference(['B8', 'B4'])
    return image.addBands(ndvi.rename('NDVI'))

# Map the ndvi function over the filtered collection
filtered_with_ndvi = chained_filter.map(compute_ndvi)

print(filtered_with_ndvi.first().bandNames().getInfo())

Reduce

The last step of the Filter-Map-Reduce process is Reducing. Reducers in GEE are just functions that aggregate data.

In this example, we will get the mean value of the NDVI band for each pixel in each Image in the Image Collection and display the results.

# Visualization parameters for NDVI
vis_ndvi = {
    'min': 0,
    'max': 0.8,
    'palette': ['black', 'blue', 'orange', 'green']
}

# Reduce the collection to mean, min, and max NDVI
mean_ndvi = filtered_with_ndvi.select(['NDVI']).reduce(ee.Reducer.mean())
min_ndvi = filtered_with_ndvi.select(['NDVI']).reduce(ee.Reducer.min())
max_ndvi = filtered_with_ndvi.select(['NDVI']).reduce(ee.Reducer.max())

# Create a geemap Map and center on the point
m2 = geemap.Map()
m2.centerObject(pp_point, 10)
m2.addLayer(mean_ndvi, vis_ndvi, 'Mean NDVI')
m2.addLayer(min_ndvi, vis_ndvi, 'Min NDVI')
m2.addLayer(max_ndvi, vis_ndvi, 'Max NDVI')

m2

Time-series data

Temporal filter

import ee
import geemap
import pandas as pd
import matplotlib.pyplot as plt

# Define the point of interest in Puerto Princesa (e.g. airport)
pp_point = ee.Geometry.Point(118.75767702, 9.7421659)

chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/PENTAD")
start_date = '2025-01-01'
end_date = '2025-12-31'
chirps_2025 = chirps.filter(ee.Filter.date(start_date, end_date))

# Count the number of images
print(chirps_2025.size().getInfo())

Aggregation

year = 2025

# Create a list for months of 2024 (Jan-Sep)
months_2025 = ee.List.sequence(1, 12)

# Write a function that takes a month number and returns a monthly image
def create_monthly_image(start_month):
    start = ee.Date.fromYMD(year, start_month, 1)
    end = start.advance(1, 'month')
    month_filtered = chirps_2025.filter(ee.Filter.date(start, end))

    # Calculate total precipitation
    total = month_filtered.reduce(ee.Reducer.sum())

    # Set metadata for the monthly image
    return total.set({
        'system:time_start': start.millis(),
        'system:time_end': end.millis(),
        'year': year,
        'month': start_month
    })

# Map the function on the list of months
monthly_images = months_2025.map(create_monthly_image)

# Create an ee.ImageCollection from the monthly images
monthly_collection = ee.ImageCollection.fromImages(monthly_images)
print(monthly_collection.getInfo()['features'])
# View the monthly images
m3 = geemap.Map()

month_names = {
    1: 'January', 2: 'February', 3: 'March', 4: 'April',
    5: 'May', 6: 'June', 7: 'July', 8: 'August', 9: 'September',
    10: 'October', 11: 'November', 12: 'December'
}

vis_params = {                                            
  'min': 0,
  'max': 500,
  'palette': ['white', 'green', 'blue']
}

# Bounding box of the Philippines
ph_bounds = ee.Geometry.Rectangle([116.928, 4.587, 126.537, 21.321])
m3.centerObject(ph_bounds, 6)

# for m in range(1, 13):
#       img = monthly_collection.filter(ee.Filter.eq('month', m)).first()
#       m3.addLayer(img, vis_params, f'Month {m}')

# m3

# Add all monthly images, clipped to the Philippines bounding box
for m in range(1, 13):
    img = monthly_collection.filter(ee.Filter.eq('month', m)).first()
    clipped = img.clip(ph_bounds)
    m3.addLayer(clipped, vis_params, f'{month_names[m]} 2025')

m3
# Add plot

# Extract precipitation values at the point and build a DataFrame for plotting
# ui.Chart is a Code Editor-only UI component — not available in Python
# Use getRegion() to extract values, then plot with matplotlib
region_data = monthly_collection.getRegion(pp_point, scale=5566).getInfo()
headers = region_data[0]
rows = region_data[1:]
df = pd.DataFrame(rows, columns=headers)
df['date'] = pd.to_datetime(df['time'], unit='ms')
df = df[['date', 'precipitation_sum']].dropna()
df = df.sort_values('date')

# Plot monthly rainfall
plt.figure(figsize=(10, 5))
plt.plot(df['date'], df['precipitation_sum'], marker='o', linewidth=1)
plt.title('Monthly Rainfall at Puerto Princesa Airport')
plt.xlabel('Month')
plt.ylabel('Rainfall (mm)')
plt.tight_layout()
plt.show()

Combining images

"""
Mosaics and Composites are two simple ways to combine Images in an Image Collection (IC)
that returns 1 Image that you can process.

Mosaic (.mosaic()) just shows an IC with the latest pixels on top. This is the default
behaviour when loading ICs.

Composites allow us to define how the pixels of the Images in an IC will be combined.
For example, we can use the median() function so that the resulting pixel value is the
median of all pixels from the IC.
"""

# Let's compare a mosaic and a composite using a filtered ImageCollection

import ee
import geemap

ee.Initialize()

# Set visualization parameters for RGB display
vis_rgb = {
    'min': 0.0,
    'max': 3000,
    'bands': ['B4', 'B3', 'B2'],  # Red, Green, Blue bands for S2
}

# Define the point geometry for pseudo Puerto Princesa center
pp_point = ee.Geometry.Point([118.75, 9.78])

Map = geemap.Map()
Map.centerObject(pp_point, 10)

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

# Chained filter
chained_filter = (s2.filterBounds(pp_point)
                    .filterDate('2025-01-01', '2025-12-31')
                    .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20))

mosaic = chained_filter.mosaic()

median_composite = chained_filter.median()

Map.addLayer(chained_filter.first(), vis_rgb, 'First')
Map.addLayer(mosaic, vis_rgb, 'Mosaic')
Map.addLayer(median_composite, vis_rgb, 'Median Composite')

Map

Bands and indices

"""
Spectral Indices are crucial in remote sensing with some of the most basic forms
being Normalized Differences (i.e. [Band 1 - Band 2] / [Band 1 + Band 2])

We can manually compute for these Normalized Difference indices but we can also use
GEE's built-in normalizedDifference function.

In this example, let's compute for the NDVI (Normalized Difference Vegetation Index),
NDWI (Normalized Difference Water Index), and NDBI (Normalized Difference Built-up Index).
"""

import ee
import geemap

ee.Initialize()

# Visualization parameters
vis_rgb = {
    'min': 0.0,
    'max': 3000,
    'bands': ['B4', 'B3', 'B2'],  # Red, Green, Blue bands for S2
}

vis_ndvi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'green']
}

vis_ndwi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'blue']
}

vis_nbdi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'red']
}

# Define the point geometry for pseudo Puerto Princesa center
pp_point = ee.Geometry.Point([118.75, 9.78])

Map = geemap.Map()
Map.centerObject(pp_point, 10)

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

# Chained filter sort by CLOUDY_PIXEL_PERCENTAGE (increasing)
filtered_collection = (s2.filterBounds(pp_point)
                         .filterDate('2024-01-01', '2024-12-31')
                         .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)
                         .sort('CLOUDY_PIXEL_PERCENTAGE'))

# print(filtered_collection.getInfo())

# Get the first image (lowest CLOUDY_PIXEL_PERCENTAGE)
image = filtered_collection.first()

# Compute for the normalized difference indices
# (NIR - Red) / (NIR + Red)
ndvi = image.normalizedDifference(['B8', 'B4']).rename(['NDVI'])

# (Green - NIR) / (Green + NIR)
ndwi = image.normalizedDifference(['B3', 'B8']).rename(['NDWI'])

# (SWIR - NIR) / (SWIR + NIR)
nbdi = image.normalizedDifference(['B11', 'B8']).rename(['NBDI'])

# Add the indices to the map
Map.addLayer(image, vis_rgb, 'Sentinel-2')
Map.addLayer(ndvi, vis_ndvi, 'NDVI')
Map.addLayer(ndwi, vis_ndwi, 'NDWI')
Map.addLayer(nbdi, vis_nbdi, 'NBDI')

Map

Thresholds

"""
Applying a threshold is an easy way to reclassify an image. For example, if we just want to show which areas
are possible vegetation in Puerto Princesa, we can use a threshold to generalize the NDVI value as either
vegetation or non-vegetation. We can do this by deciding on a minimum value to use for vegetation, for example: 0.6
"""

import ee
import geemap

ee.Initialize()

# Visualization parameters
vis_veg = {
    'min': 0,
    'max': 1,
    'palette': ['white', 'green']
}

vis_ndvi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'green']
}

# Define the point geometry for pseudo Puerto Princesa center
pp_point = ee.Geometry.Point([118.75, 9.78])

Map = geemap.Map()
Map.centerObject(pp_point, 10)

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

# Chained filter sort by CLOUDY_PIXEL_PERCENTAGE (increasing)
filtered_collection = (s2.filterBounds(pp_point)
                         .filterDate('2024-01-01', '2024-12-31')
                         .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)
                         .sort('CLOUDY_PIXEL_PERCENTAGE'))

# Get the image with the lowest CLOUDY_PIXEL_PERCENTAGE
image = filtered_collection.first()

# Compute for the NDVI
# (NIR - Red) / (NIR + Red)
ndvi = image.normalizedDifference(['B8', 'B4']).rename(['NDVI'])

# Apply a threshold of 0.6, those higher = 1, those lower = 0
veg_ndvi = ndvi.gt(0.6)

# Add the indices to the map
Map.addLayer(ndvi, vis_ndvi, 'NDVI')
Map.addLayer(veg_ndvi, vis_veg, 'Vegetation vs Non-vegetation')

Map

Masks

"""
Masking involves the removal of specific areas in an image from being displayed or processed by using a mask.
The .mask() function in Earth Engine allows you to view and update a mask.
"""

import ee
import geemap

ee.Initialize()

# Visualization parameters
vis_ndvi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'green']
}

vis_ndwi = {
    'min': -0.2,
    'max': 0.8,
    'palette': ['white', 'blue']
}

# Define the point geometry for pseudo Puerto Princesa center
pp_point = ee.Geometry.Point([118.75, 9.78])

Map = geemap.Map()
Map.centerObject(pp_point, 10)

# Load the Sentinel-2 harmonized image collection
s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')

# Chained filter sort by CLOUDY_PIXEL_PERCENTAGE (increasing)
filtered_collection = (s2.filterBounds(pp_point)
                         .filterDate('2024-01-01', '2024-12-31')
                         .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)
                         .sort('CLOUDY_PIXEL_PERCENTAGE'))

# print(filtered_collection.getInfo())

# Get the first image (lowest CLOUDY_PIXEL_PERCENTAGE)
image = filtered_collection.first()

# Compute for the normalized difference indices
# (NIR - Red) / (NIR + Red)
ndvi = image.normalizedDifference(['B8', 'B4']).rename(['NDVI'])

# (Green - NIR) / (Green + NIR)
ndwi = image.normalizedDifference(['B3', 'B8']).rename(['NDWI'])

# Create a mask using the NDWI layer. .lt() means less than. This masks the water pixels (> 0).
water_mask = ndwi.lt(0)

# Mask the NDVI using the water_mask
masked_ndvi = ndvi.updateMask(water_mask)

# Add the layers to the map
Map.addLayer(ndvi, vis_ndvi, 'NDVI')
Map.addLayer(water_mask, {}, 'Water Mask')
Map.addLayer(masked_ndvi, vis_ndvi, 'Masked NDVI')

Map
GEE
Feature data (geometries and attributes)
GEE
Comparing Machine Learning Classifiers for Crop Mapping in the Philippines