Time Series Analysis

Learn to analyze temporal patterns and changes using Earth Engine image collections.

What You’ll Learn

  • Loading and filtering image collections

  • Creating time series datasets

  • Temporal analysis and trend detection

  • Seasonal pattern analysis

  • Change detection techniques

  • Time series visualization

Prerequisites

  • Understanding of image collections

  • Knowledge of filtering techniques

  • Basic statistical concepts

  • Familiarity with temporal data

Loading Time Series Data

Basic time series setup:

import ee

# Initialize Earth Engine
ee.Initialize(project='your-project-id')

# Define study area and time range
geometry = ee.Geometry.Point([-122.4, 37.8])  # San Francisco
start_date = '2020-01-01'
end_date = '2023-12-31'

# Load Landsat collection
collection = (ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
             .filterDate(start_date, end_date)
             .filterBounds(geometry)
             .filter(ee.Filter.lt('CLOUD_COVER', 20)))

print(f"Found {collection.size().getInfo()} images")

Adding temporal properties:

def add_date_properties(image):
    """Add useful date properties to each image."""
    date = ee.Date(image.get('system:time_start'))
    return image.set({
        'year': date.get('year'),
        'month': date.get('month'),
        'day_of_year': date.getRelative('day', 'year'),
        'decimal_year': date.difference(ee.Date('1970-01-01'), 'year')
    })

# Apply to collection
collection = collection.map(add_date_properties)

Creating Time Series

Extract time series values:

def extract_time_series(image_collection, geometry, scale=30):
    """Extract time series data from image collection."""

    def calculate_indices(image):
        """Calculate spectral indices for time series."""
        # Apply scale factors
        scaled = image.select('SR_B.').multiply(0.0000275).add(-0.2)

        # Calculate indices
        ndvi = scaled.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI')
        ndwi = scaled.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI')
        evi = scaled.expression(
            '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))',
            {
                'NIR': scaled.select('SR_B5'),
                'RED': scaled.select('SR_B4'),
                'BLUE': scaled.select('SR_B2')
            }
        ).rename('EVI')

        return image.addBands([ndvi, ndwi, evi])

    # Calculate indices for all images
    with_indices = image_collection.map(calculate_indices)

    # Extract values
    def extract_values(image):
        # Calculate mean values over geometry
        values = image.select(['NDVI', 'NDWI', 'EVI']).reduceRegion(
            reducer=ee.Reducer.mean(),
            geometry=geometry,
            scale=scale,
            maxPixels=1e9
        )

        # Return feature with date and values
        return ee.Feature(None, {
            'date': image.get('system:time_start'),
            'year': image.get('year'),
            'month': image.get('month'),
            'day_of_year': image.get('day_of_year'),
            'decimal_year': image.get('decimal_year'),
            'ndvi': values.get('NDVI'),
            'ndwi': values.get('NDWI'),
            'evi': values.get('EVI'),
            'cloud_cover': image.get('CLOUD_COVER')
        })

    return with_indices.map(extract_values)

Temporal Compositing

Monthly composites:

def create_monthly_composites(collection, start_year, end_year):
    """Create monthly median composites."""

    months = ee.List.sequence(1, 12)
    years = ee.List.sequence(start_year, end_year)

    def create_monthly_composite(year):
        def create_month_composite(month):
            # Filter to specific year and month
            monthly = collection.filter(
                ee.Filter.calendarRange(year, year, 'year')
            ).filter(
                ee.Filter.calendarRange(month, month, 'month')
            )

            # Create composite
            composite = monthly.median()

            # Add date information
            date = ee.Date.fromYMD(year, month, 1)
            return composite.set({
                'year': year,
                'month': month,
                'system:time_start': date.millis(),
                'image_count': monthly.size()
            })

        return months.map(create_month_composite)

    # Create all composites
    composites = years.map(create_monthly_composite).flatten()
    return ee.ImageCollection.fromImages(composites)

Seasonal composites:

def create_seasonal_composites(collection):
    """Create seasonal composites."""

    # Define seasons by day of year
    seasons = {
        'spring': [80, 171],   # March 21 - June 20
        'summer': [172, 264],  # June 21 - September 21
        'fall': [265, 354],    # September 22 - December 20
        'winter': [355, 79]    # December 21 - March 20 (spans years)
    }

    seasonal_composites = []

    for season_name, day_range in seasons.items():
        if season_name == 'winter':
            # Handle winter spanning years
            winter_images = collection.filter(
                ee.Filter.Or(
                    ee.Filter.gte('day_of_year', day_range[0]),
                    ee.Filter.lte('day_of_year', day_range[1])
                )
            )
        else:
            winter_images = collection.filter(
                ee.Filter.And(
                    ee.Filter.gte('day_of_year', day_range[0]),
                    ee.Filter.lte('day_of_year', day_range[1])
                )
            )

        composite = winter_images.median().set('season', season_name)
        seasonal_composites.append(composite)

    return seasonal_composites

Trend Analysis

Linear trend calculation:

def calculate_temporal_trend(collection, band_name):
    """Calculate linear trend for a specific band."""

    # Add time variable (years since start)
    def add_time_band(image):
        time_start = ee.Date(image.get('system:time_start'))
        years_since_start = time_start.difference(ee.Date('2020-01-01'), 'year')
        return image.addBands(ee.Image(years_since_start).rename('time'))

    # Add time band to all images
    collection_with_time = collection.map(add_time_band)

    # Calculate linear regression
    linear_fit = collection_with_time.select(['time', band_name]).reduce(
        ee.Reducer.linearFit()
    )

    return linear_fit

Harmonic analysis for seasonality:

def harmonic_analysis(collection, band_name):
    """Perform harmonic analysis to detect seasonal patterns."""

    def add_harmonic_terms(image):
        # Get time in fractional years
        time_start = ee.Date(image.get('system:time_start'))
        t = time_start.difference(ee.Date('2020-01-01'), 'year')

        # Add harmonic terms
        omega = 2 * 3.14159
        cos_1 = t.multiply(omega).cos().rename('cos_1')
        sin_1 = t.multiply(omega).sin().rename('sin_1')
        cos_2 = t.multiply(omega * 2).cos().rename('cos_2')
        sin_2 = t.multiply(omega * 2).sin().rename('sin_2')

        return image.addBands([
            ee.Image(t).rename('time'),
            cos_1, sin_1, cos_2, sin_2
        ])

    # Add harmonic terms
    harmonic_collection = collection.map(add_harmonic_terms)

    # Perform harmonic regression
    harmonic_fit = harmonic_collection.select([
        'time', 'cos_1', 'sin_1', 'cos_2', 'sin_2', band_name
    ]).reduce(ee.Reducer.linearRegression({
        'numX': 5,
        'numY': 1
    }))

    return harmonic_fit

Change Detection

Simple change detection:

def detect_change(before_collection, after_collection, threshold=0.1):
    """Simple change detection between two time periods."""

    # Create composites
    before_composite = before_collection.median()
    after_composite = after_collection.median()

    # Calculate NDVI for both periods
    before_ndvi = before_composite.normalizedDifference(['SR_B5', 'SR_B4'])
    after_ndvi = after_composite.normalizedDifference(['SR_B5', 'SR_B4'])

    # Calculate change
    change = after_ndvi.subtract(before_ndvi)

    # Classify change
    change_mask = change.abs().gt(threshold)
    increase = change.gt(threshold)
    decrease = change.lt(-threshold)

    return {
        'change': change,
        'change_mask': change_mask,
        'increase': increase,
        'decrease': decrease
    }

LandTrendr for advanced change detection:

def run_landtrendr(collection):
    """Run LandTrendr algorithm for change detection."""

    # Prepare collection for LandTrendr
    def prepare_image(image):
        scaled = image.select('SR_B.').multiply(0.0000275).add(-0.2)
        ndvi = scaled.normalizedDifference(['SR_B5', 'SR_B4']).multiply(1000)
        return image.addBands(ndvi.rename('NDVI')).select('NDVI')

    lt_collection = collection.map(prepare_image)

    # LandTrendr parameters
    lt_params = {
        'maxSegments': 6,
        'spikeThreshold': 0.9,
        'vertexCountOvershoot': 3,
        'preventOneYearRecovery': True,
        'recoveryThreshold': 0.25,
        'pvalThreshold': 0.05,
        'bestModelProportion': 0.75,
        'minObservationsNeeded': 6
    }

    # Run LandTrendr
    lt_result = ee.Algorithms.TemporalSegmentation.LandTrendr(
        timeSeries=lt_collection,
        **lt_params
    )

    return lt_result

Complete Example

Here’s the complete time series analysis example:

  1"""
  2Intermediate Example 1: Time Series Analysis
  3============================================
  4
  5This example demonstrates:
  6- Working with image collections
  7- Filtering by date and location
  8- Time series chart creation
  9- Statistical analysis over time
 10
 11Use case: Analyzing NDVI trends over agricultural areas
 12"""
 13
 14import ee
 15import matplotlib.pyplot as plt
 16import pandas as pd
 17from datetime import datetime, timedelta
 18
 19def analyze_ndvi_time_series(geometry, start_date, end_date):
 20    """
 21    Analyze NDVI time series for a given geometry and date range.
 22    
 23    Args:
 24        geometry: ee.Geometry object defining the area of interest
 25        start_date: Start date as string ('YYYY-MM-DD')
 26        end_date: End date as string ('YYYY-MM-DD')
 27    
 28    Returns:
 29        Dictionary containing time series data and statistics
 30    """
 31    
 32    # Load Landsat 8 Surface Reflectance collection
 33    collection = (ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
 34                 .filterDate(start_date, end_date)
 35                 .filterBounds(geometry)
 36                 .filter(ee.Filter.lt('CLOUD_COVER', 20)))
 37    
 38    print(f"Found {collection.size().getInfo()} images in the collection")
 39    
 40    # Function to calculate NDVI and add time properties
 41    def add_ndvi_and_time(image):
 42        # Calculate NDVI
 43        ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI')
 44        
 45        # Add time properties
 46        date = ee.Date(image.get('system:time_start'))
 47        years = date.difference(ee.Date('1970-01-01'), 'year')
 48        
 49        return (image.addBands(ndvi)
 50                    .set('date', date)
 51                    .set('year', date.get('year'))
 52                    .set('month', date.get('month'))
 53                    .set('day_of_year', date.getRelative('day', 'year'))
 54                    .set('decimal_year', years))
 55    
 56    # Apply NDVI calculation to collection
 57    ndvi_collection = collection.map(add_ndvi_and_time)
 58    
 59    # Create time series by reducing each image to mean NDVI
 60    def extract_ndvi_value(image):
 61        # Calculate mean NDVI over the geometry
 62        ndvi_mean = image.select('NDVI').reduceRegion(
 63            reducer=ee.Reducer.mean(),
 64            geometry=geometry,
 65            scale=30,
 66            maxPixels=1e9
 67        )
 68        
 69        # Return feature with date and NDVI value
 70        return ee.Feature(None, {
 71            'date': image.get('date'),
 72            'decimal_year': image.get('decimal_year'),
 73            'ndvi': ndvi_mean.get('NDVI'),
 74            'year': image.get('year'),
 75            'month': image.get('month'),
 76            'day_of_year': image.get('day_of_year')
 77        })
 78    
 79    # Extract time series data
 80    time_series = ndvi_collection.map(extract_ndvi_value)
 81    
 82    # Convert to pandas DataFrame for analysis
 83    time_series_list = time_series.getInfo()['features']
 84    df_data = []
 85    
 86    for feature in time_series_list:
 87        props = feature['properties']
 88        if props['ndvi'] is not None:  # Filter out null values
 89            df_data.append({
 90                'date': datetime.fromtimestamp(props['date']['value'] / 1000),
 91                'decimal_year': props['decimal_year'],
 92                'ndvi': props['ndvi'],
 93                'year': props['year'],
 94                'month': props['month'],
 95                'day_of_year': props['day_of_year']
 96            })
 97    
 98    df = pd.DataFrame(df_data)
 99    df = df.sort_values('date')
100    
101    # Calculate statistics
102    stats = {
103        'mean_ndvi': df['ndvi'].mean(),
104        'std_ndvi': df['ndvi'].std(),
105        'min_ndvi': df['ndvi'].min(),
106        'max_ndvi': df['ndvi'].max(),
107        'data_points': len(df)
108    }
109    
110    return df, stats
111
112def plot_time_series(df, stats, title="NDVI Time Series"):
113    """
114    Create time series plot with trend analysis.
115    """
116    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
117    
118    # Main time series plot
119    ax1.plot(df['date'], df['ndvi'], 'o-', linewidth=1, markersize=3, alpha=0.7)
120    ax1.set_title(f"{title}\nMean: {stats['mean_ndvi']:.3f} ± {stats['std_ndvi']:.3f}")
121    ax1.set_xlabel('Date')
122    ax1.set_ylabel('NDVI')
123    ax1.grid(True, alpha=0.3)
124    ax1.axhline(y=stats['mean_ndvi'], color='r', linestyle='--', alpha=0.5, label='Mean')
125    ax1.legend()
126    
127    # Monthly aggregation
128    monthly_mean = df.groupby(df['date'].dt.month)['ndvi'].mean()
129    ax2.bar(monthly_mean.index, monthly_mean.values, alpha=0.7, color='green')
130    ax2.set_title('Monthly Average NDVI')
131    ax2.set_xlabel('Month')
132    ax2.set_ylabel('Mean NDVI')
133    ax2.set_xticks(range(1, 13))
134    ax2.grid(True, alpha=0.3)
135    
136    plt.tight_layout()
137    plt.show()
138    
139    return fig
140
141def main():
142    """
143    Main function demonstrating time series analysis.
144    """
145    # Initialize Earth Engine
146    try:
147        ee.Initialize(project='your-project-id')
148        print("✓ Earth Engine initialized successfully!")
149    except Exception as e:
150        print(f"✗ Error initializing Earth Engine: {e}")
151        return
152    
153    # Define area of interest (example: agricultural area in Iowa)
154    geometry = ee.Geometry.Rectangle([-94.6, 41.9, -94.4, 42.1])
155    
156    # Define date range
157    start_date = '2020-01-01'
158    end_date = '2023-12-31'
159    
160    print(f"Analyzing NDVI time series from {start_date} to {end_date}")
161    print("This may take a few minutes...")
162    
163    # Perform analysis
164    df, stats = analyze_ndvi_time_series(geometry, start_date, end_date)
165    
166    # Print results
167    print("\n📊 Time Series Analysis Results:")
168    print(f"Data points: {stats['data_points']}")
169    print(f"Mean NDVI: {stats['mean_ndvi']:.3f}")
170    print(f"Standard deviation: {stats['std_ndvi']:.3f}")
171    print(f"Range: {stats['min_ndvi']:.3f} to {stats['max_ndvi']:.3f}")
172    
173    # Create visualization
174    plot_time_series(df, stats, "Agricultural Area NDVI Time Series")
175    
176    # Seasonal analysis
177    print("\n🌱 Seasonal Analysis:")
178    seasonal_stats = df.groupby(df['date'].dt.month)['ndvi'].agg(['mean', 'std', 'count'])
179    print(seasonal_stats)
180    
181    print("\n✅ Time series analysis completed successfully!")
182
183if __name__ == "__main__":
184    main()

Visualization Techniques

Time series plots:

# Extract data for plotting
time_series_fc = extract_time_series(collection, geometry)
time_series_list = time_series_fc.getInfo()['features']

# Convert to pandas DataFrame
import pandas as pd
from datetime import datetime

data = []
for feature in time_series_list:
    props = feature['properties']
    if props.get('ndvi') is not None:
        data.append({
            'date': datetime.fromtimestamp(props['date']['value'] / 1000),
            'ndvi': props['ndvi'],
            'evi': props['evi'],
            'cloud_cover': props['cloud_cover']
        })

df = pd.DataFrame(data)

# Create plots
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# NDVI time series
ax1.plot(df['date'], df['ndvi'], 'g-', linewidth=1, alpha=0.7)
ax1.set_title('NDVI Time Series')
ax1.set_ylabel('NDVI')
ax1.grid(True, alpha=0.3)

# Monthly averages
monthly_avg = df.groupby(df['date'].dt.month)['ndvi'].mean()
ax2.bar(monthly_avg.index, monthly_avg.values, alpha=0.7)
ax2.set_title('Monthly Average NDVI')
ax2.set_xlabel('Month')
ax2.set_ylabel('NDVI')

plt.tight_layout()
plt.show()

Best Practices

Data Quality: * Apply cloud masking consistently * Filter by metadata quality indicators * Handle missing data appropriately

Temporal Sampling: * Consider sensor revisit frequency * Account for seasonal variations * Use appropriate temporal windows

Scale Considerations: * Match analysis scale to phenomena * Consider mixed pixel effects * Validate with ground truth data

Next: Image Collection Filtering