Timeseries gap filling

Uncomment the following line to install the latest version of geeagri if needed.

In [ ]:

# !pip install -U geeagri

# !pip install -U geeagri

In [ ]:

import ee
import geemap
from geeagri.preprocessing import Sentinel2CloudMask, TemporalInterpolation
from geeagri.extract import extract_timeseries_to_point

import matplotlib.pyplot as plt
import seaborn as sns

from datetime import datetime
from PIL import Image, ImageDraw, ImageFont, ImageSequence

plt.rcParams["font.family"] = "DeJavu Serif"
plt.rcParams["font.serif"] = "Times New Roman"

import ee import geemap from geeagri.preprocessing import Sentinel2CloudMask, TemporalInterpolation from geeagri.extract import extract_timeseries_to_point import matplotlib.pyplot as plt import seaborn as sns from datetime import datetime from PIL import Image, ImageDraw, ImageFont, ImageSequence plt.rcParams["font.family"] = "DeJavu Serif" plt.rcParams["font.serif"] = "Times New Roman"

Initialize a Map object¶

Authenticate and initialize Earth Engine. If it doesn't work, specify a project name

In [ ]:

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

Map = geemap.Map(basemap="SATELLITE")
Map

# ee.Authenticate() # ee.Initialize(project='your-project-id') Map = geemap.Map(basemap="SATELLITE") Map

Import region of interest¶

In [ ]:

bbox = [-98.451233, 38.430732, -98.274765, 38.523996]
region = ee.Geometry.BBox(*bbox)
region_style = {"color": "red", "width": 1}
Map.addLayer(region, region_style, "Region")
Map.centerObject(region, 13)

bbox = [-98.451233, 38.430732, -98.274765, 38.523996] region = ee.Geometry.BBox(*bbox) region_style = {"color": "red", "width": 1} Map.addLayer(region, region_style, "Region") Map.centerObject(region, 13)

Import Cropland Data Layer (CDL)¶

In [ ]:

croplandcover = (
    ee.ImageCollection("USDA/NASS/CDL")
    .filterDate("2020-01-01", "2021-01-01")
    .first()
    .clip(region)
)

cultivated = croplandcover.select("cultivated").eq(2).selfMask()
croplandcover = croplandcover.select("cropland")
Map.addLayer(cultivated, {}, "Cultivated")
Map.addLayer(croplandcover, {}, "Crop Landcover")

croplandcover = ( ee.ImageCollection("USDA/NASS/CDL") .filterDate("2020-01-01", "2021-01-01") .first() .clip(region) ) cultivated = croplandcover.select("cultivated").eq(2).selfMask() croplandcover = croplandcover.select("cropland") Map.addLayer(cultivated, {}, "Cultivated") Map.addLayer(croplandcover, {}, "Crop Landcover")

Load Sentinel-2 image collection and mask clouds and shadows¶

In [ ]:

s2_cloud_masker = Sentinel2CloudMask(
    region=region,
    start_date="2020-01-01",
    end_date="2021-01-01",
    cloud_filter=60,  # maximum scene-level cloudiness allowed (%)
    cloud_prob_threshold=50,  # cloud probability threshold (values above are considered clouds)
    nir_dark_threshold=0.15,  # NIR reflectance threshold (values below considered potential shadows)
    shadow_proj_dist=1,  # maximum distance (km) to search for shadows from clouds.
    buffer=50,  # buffer distance (m) to dilate cloud/shadow masks.
)

s2_cloud_masked = s2_cloud_masker.get_cloudfree_collection()

s2_cloud_masker = Sentinel2CloudMask( region=region, start_date="2020-01-01", end_date="2021-01-01", cloud_filter=60, # maximum scene-level cloudiness allowed (%) cloud_prob_threshold=50, # cloud probability threshold (values above are considered clouds) nir_dark_threshold=0.15, # NIR reflectance threshold (values below considered potential shadows) shadow_proj_dist=1, # maximum distance (km) to search for shadows from clouds. buffer=50, # buffer distance (m) to dilate cloud/shadow masks. ) s2_cloud_masked = s2_cloud_masker.get_cloudfree_collection()

Calculate Normalized Difference Vegetation Index (NDVI)¶

In [ ]:

def calculateNDVI(image):
    ndvi = image.expression(
        "(NIR - Red) / (NIR + Red)",
        {"NIR": image.select("B8"), "Red": image.select("B4")},
    ).copyProperties(image, ["system:time_start"])

    ndvi = ee.Image(ndvi).rename("ndvi").clip(region)

    return ndvi

ndvi_col = s2_cloud_masked.map(calculateNDVI)

def calculateNDVI(image): ndvi = image.expression( "(NIR - Red) / (NIR + Red)", {"NIR": image.select("B8"), "Red": image.select("B4")}, ).copyProperties(image, ["system:time_start"]) ndvi = ee.Image(ndvi).rename("ndvi").clip(region) return ndvi ndvi_col = s2_cloud_masked.map(calculateNDVI)

Timeseries gap filling using Temporal Interpolation¶

In [ ]:

# Initialize a 'TemporalInterpolation' object
temp_interp = TemporalInterpolation(
    image_collection=ndvi_col,  # ee.ImageCollection of NDVI
    window=45,  # Temporal window in days
)

# Get the interpolation collection
ndvi_interpolated = temp_interp.get_interpolated_collection()

# Initialize a 'TemporalInterpolation' object temp_interp = TemporalInterpolation( image_collection=ndvi_col, # ee.ImageCollection of NDVI window=45, # Temporal window in days ) # Get the interpolation collection ndvi_interpolated = temp_interp.get_interpolated_collection()

Extract timelapse for raw and interpolated NDVI image collection¶

In [ ]:

# Function to add timestamp of gif frames
def add_dates_to_gif(
    input_gif: str,
    output_gif: str,
    dates: list[str],
    font_path: str = "DejaVuSans-Bold.ttf",
    font_size: int = 20,
    position: tuple = (20, 20),
    color: str = "red",
):
    """
    Overlay a list of dates onto each frame of an animated GIF.

    Args:
        input_gif (str): Path to the input GIF file.
        output_gif (str): Path where the output GIF will be saved.
        dates (list[str]): List of date strings (one per frame).
        font_path (str, optional): Path to a .ttf font file. Defaults to "DejaVuSans-Bold.ttf".
        font_size (int, optional): Font size. Defaults to 40.
        position (tuple, optional): (x, y) position of the text. Defaults to (20, 20).
        color (str, optional): Text color. Defaults to "black".
    """
    font = ImageFont.truetype(font_path, font_size)

    with Image.open(input_gif) as im:
        frames = []

        for i, frame in enumerate(ImageSequence.Iterator(im)):
            frame = frame.convert("RGB")
            draw = ImageDraw.Draw(frame)

            if i < len(dates):
                draw.text(position, dates[i], fill=color, font=font)

            frames.append(frame)

        # Save the annotated frames as a new GIF
        frames[0].save(
            output_gif,
            save_all=True,
            append_images=frames[1:],
            duration=im.info.get("duration", 200),  # default 200ms if not found
            loop=0,
        )

# Function to add timestamp of gif frames def add_dates_to_gif( input_gif: str, output_gif: str, dates: list[str], font_path: str = "DejaVuSans-Bold.ttf", font_size: int = 20, position: tuple = (20, 20), color: str = "red", ): """ Overlay a list of dates onto each frame of an animated GIF. Args: input_gif (str): Path to the input GIF file. output_gif (str): Path where the output GIF will be saved. dates (list[str]): List of date strings (one per frame). font_path (str, optional): Path to a .ttf font file. Defaults to "DejaVuSans-Bold.ttf". font_size (int, optional): Font size. Defaults to 40. position (tuple, optional): (x, y) position of the text. Defaults to (20, 20). color (str, optional): Text color. Defaults to "black". """ font = ImageFont.truetype(font_path, font_size) with Image.open(input_gif) as im: frames = [] for i, frame in enumerate(ImageSequence.Iterator(im)): frame = frame.convert("RGB") draw = ImageDraw.Draw(frame) if i < len(dates): draw.text(position, dates[i], fill=color, font=font) frames.append(frame) # Save the annotated frames as a new GIF frames[0].save( output_gif, save_all=True, append_images=frames[1:], duration=im.info.get("duration", 200), # default 200ms if not found loop=0, )

In [ ]:

# Get the image dates
dates = ndvi_col.aggregate_array("system:time_start").getInfo()
dates_str = [datetime.fromtimestamp(d / 1000).strftime("%d %b %Y") for d in dates]

video_args = {
    "dimensions": 1000,
    "region": region,
    "framesPerSecond": 2,
    "crs": "EPSG:4326",
    "min": -1,
    "max": 1,
    "palette": [
        "#a50026",
        "#d73027",
        "#f46d43",
        "#fdae61",
        "#fee08b",
        "#d9ef8b",
        "#a6d96a",
        "#66bd63",
        "#1a9850",
        "#006837",
    ],
}

saved_ndvi_gif = "ndvi.gif"
geemap.download_ee_video(ndvi_col, video_args, saved_ndvi_gif)

# add timestamps and plot the timelapse of raw NDVI
add_dates_to_gif(saved_ndvi_gif, saved_ndvi_gif, dates_str)
geemap.show_image(saved_ndvi_gif)

# Get the image dates dates = ndvi_col.aggregate_array("system:time_start").getInfo() dates_str = [datetime.fromtimestamp(d / 1000).strftime("%d %b %Y") for d in dates] video_args = { "dimensions": 1000, "region": region, "framesPerSecond": 2, "crs": "EPSG:4326", "min": -1, "max": 1, "palette": [ "#a50026", "#d73027", "#f46d43", "#fdae61", "#fee08b", "#d9ef8b", "#a6d96a", "#66bd63", "#1a9850", "#006837", ], } saved_ndvi_gif = "ndvi.gif" geemap.download_ee_video(ndvi_col, video_args, saved_ndvi_gif) # add timestamps and plot the timelapse of raw NDVI add_dates_to_gif(saved_ndvi_gif, saved_ndvi_gif, dates_str) geemap.show_image(saved_ndvi_gif)

In [ ]:

saved_ndvi_interpolated_gif = "ndvi_interpolated.gif"
geemap.download_ee_video(ndvi_interpolated, video_args, saved_ndvi_interpolated_gif)

# add timestamps and plot the timelapse of interpolated NDVI
add_dates_to_gif(saved_ndvi_interpolated_gif, saved_ndvi_interpolated_gif, dates_str)
geemap.show_image(saved_ndvi_interpolated_gif)

saved_ndvi_interpolated_gif = "ndvi_interpolated.gif" geemap.download_ee_video(ndvi_interpolated, video_args, saved_ndvi_interpolated_gif) # add timestamps and plot the timelapse of interpolated NDVI add_dates_to_gif(saved_ndvi_interpolated_gif, saved_ndvi_interpolated_gif, dates_str) geemap.show_image(saved_ndvi_interpolated_gif)

Plot the raw and interpolated NDVI for three crop samples¶

In [ ]:

wheat_loc = (-98.39630126953126, 38.46931531751283)
soybean_loc = (-98.34480285644533, 38.50008653288019)
corn_loc = (-98.3169937133789, 38.49645916887533)

wheat_ndvi = extract_timeseries_to_point(
    lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_col, scale=10
)
soybean_ndvi = extract_timeseries_to_point(
    lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_col, scale=10
)
corn_ndvi = extract_timeseries_to_point(
    lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_col, scale=10
)
wheat_ndvi_interpolated = extract_timeseries_to_point(
    lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_interpolated, scale=10
)
soybean_ndvi_interpolated = extract_timeseries_to_point(
    lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_interpolated, scale=10
)
corn_ndvi_interpolated = extract_timeseries_to_point(
    lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_interpolated, scale=10
)

wheat_loc = (-98.39630126953126, 38.46931531751283) soybean_loc = (-98.34480285644533, 38.50008653288019) corn_loc = (-98.3169937133789, 38.49645916887533) wheat_ndvi = extract_timeseries_to_point( lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_col, scale=10 ) soybean_ndvi = extract_timeseries_to_point( lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_col, scale=10 ) corn_ndvi = extract_timeseries_to_point( lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_col, scale=10 ) wheat_ndvi_interpolated = extract_timeseries_to_point( lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_interpolated, scale=10 ) soybean_ndvi_interpolated = extract_timeseries_to_point( lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_interpolated, scale=10 ) corn_ndvi_interpolated = extract_timeseries_to_point( lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_interpolated, scale=10 )

In [ ]:

# Define datasets and colors
datasets = {
    "Wheat": ("#a87000", wheat_ndvi, wheat_ndvi_interpolated),
    "Soybean": ("#267300", soybean_ndvi, soybean_ndvi_interpolated),
    "Corn": ("#ffd400", corn_ndvi, corn_ndvi_interpolated),
}

fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharey=True)

# Raw NDVI
for crop, (color, raw, interp) in datasets.items():
    sns.lineplot(
        data=raw,
        x="time",
        y="ndvi",
        label=f"{crop} NDVI",
        c=color,
        marker="o",
        ax=axes[0],
    )

    # plot NaNs as crosses
    nan_points = raw[raw["ndvi"].isna()]
    axes[0].scatter(
        nan_points["time"],
        [0] * len(nan_points),
        marker="x",
        c=color,
        label=f"{crop} NaN",
    )

axes[0].set_title("Raw NDVI (with gaps)")
axes[0].set_ylabel("NDVI")

# Interpolated NDVI
for crop, (color, raw, interp) in datasets.items():
    sns.lineplot(
        data=interp,
        x="time",
        y="ndvi",
        label=f"{crop} NDVI Interpolated",
        c=color,
        marker="o",
        markersize=5,
        ax=axes[1],
    )

axes[1].set_title("Interpolated NDVI")
axes[1].set_ylabel("NDVI")

plt.tight_layout()
plt.show()

# Define datasets and colors datasets = { "Wheat": ("#a87000", wheat_ndvi, wheat_ndvi_interpolated), "Soybean": ("#267300", soybean_ndvi, soybean_ndvi_interpolated), "Corn": ("#ffd400", corn_ndvi, corn_ndvi_interpolated), } fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharey=True) # Raw NDVI for crop, (color, raw, interp) in datasets.items(): sns.lineplot( data=raw, x="time", y="ndvi", label=f"{crop} NDVI", c=color, marker="o", ax=axes[0], ) # plot NaNs as crosses nan_points = raw[raw["ndvi"].isna()] axes[0].scatter( nan_points["time"], [0] * len(nan_points), marker="x", c=color, label=f"{crop} NaN", ) axes[0].set_title("Raw NDVI (with gaps)") axes[0].set_ylabel("NDVI") # Interpolated NDVI for crop, (color, raw, interp) in datasets.items(): sns.lineplot( data=interp, x="time", y="ndvi", label=f"{crop} NDVI Interpolated", c=color, marker="o", markersize=5, ax=axes[1], ) axes[1].set_title("Interpolated NDVI") axes[1].set_ylabel("NDVI") plt.tight_layout() plt.show()