Moving window smoothing

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, MovingWindowSmoothing
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, MovingWindowSmoothing 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)

Apply moving window smoothing¶

In [ ]:

# Initialize a 'MovingWindowSmoothing' object
smoother = MovingWindowSmoothing(
    image_collection=ndvi_col,  # ee.ImageCollection of NDVI
    window=15,  # Temporal window in days
    reducer="MEDIAN",  # Reducer for smoothing ("MEAN" or "MEDIAN")
)

# Get the smoothed collectiom
ndvi_smoothed = smoother.get_smoothed_collection()

# Initialize a 'MovingWindowSmoothing' object smoother = MovingWindowSmoothing( image_collection=ndvi_col, # ee.ImageCollection of NDVI window=15, # Temporal window in days reducer="MEDIAN", # Reducer for smoothing ("MEAN" or "MEDIAN") ) # Get the smoothed collectiom ndvi_smoothed = smoother.get_smoothed_collection()

Extract timelapse for raw and smoothed 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_smoothed_gif = "ndvi_smoothed.gif"
geemap.download_ee_video(ndvi_smoothed, video_args, saved_ndvi_smoothed_gif)

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

saved_ndvi_smoothed_gif = "ndvi_smoothed.gif" geemap.download_ee_video(ndvi_smoothed, video_args, saved_ndvi_smoothed_gif) # add timestamps and plot the timelapse of smoothed NDVI add_dates_to_gif(saved_ndvi_smoothed_gif, saved_ndvi_smoothed_gif, dates_str) geemap.show_image(saved_ndvi_smoothed_gif)

Plot the raw and smoothed NDVI for three crop samples¶

In [ ]:

wheat_loc = (-98.39630126953126, 38.46931531751283)
soybean_loc = (-98.34480285644533, 38.50022087618732)
corn_loc = (-98.31716537475587, 38.496190467979176)

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_smoothed = extract_timeseries_to_point(
    lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_smoothed, scale=10
)
soybean_ndvi_smoothed = extract_timeseries_to_point(
    lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_smoothed, scale=10
)
corn_ndvi_smoothed = extract_timeseries_to_point(
    lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_smoothed, scale=10
)

wheat_loc = (-98.39630126953126, 38.46931531751283) soybean_loc = (-98.34480285644533, 38.50022087618732) corn_loc = (-98.31716537475587, 38.496190467979176) 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_smoothed = extract_timeseries_to_point( lat=wheat_loc[1], lon=wheat_loc[0], image_collection=ndvi_smoothed, scale=10 ) soybean_ndvi_smoothed = extract_timeseries_to_point( lat=soybean_loc[1], lon=soybean_loc[0], image_collection=ndvi_smoothed, scale=10 ) corn_ndvi_smoothed = extract_timeseries_to_point( lat=corn_loc[1], lon=corn_loc[0], image_collection=ndvi_smoothed, scale=10 )

In [ ]:

plt.figure(figsize=(12, 4))

sns.lineplot(
    data=wheat_ndvi,
    x="time",
    y="ndvi",
    label="Wheat NDVI",
    alpha=0.5,
    c="#a87000",
)
sns.lineplot(
    data=wheat_ndvi_smoothed,
    x="time",
    y="ndvi_median",
    label="Wheat NDVI Smoothed",
    c="#a87000",
    marker="o",
    markersize=5,
)

sns.lineplot(
    data=soybean_ndvi,
    x="time",
    y="ndvi",
    label="Soybean NDVI",
    alpha=0.5,
    c="#267300",
)
sns.lineplot(
    data=soybean_ndvi_smoothed,
    x="time",
    y="ndvi_median",
    label="Soybean NDVI Smoothed",
    c="#267300",
    marker="o",
    markersize=5,
)

sns.lineplot(
    data=corn_ndvi, x="time", y="ndvi", label="Corn NDVI", alpha=0.5, c="#ffd400"
)
sns.lineplot(
    data=corn_ndvi_smoothed,
    x="time",
    y="ndvi_median",
    label="Corn NDVI Smoothed",
    c="#ffd400",
    marker="o",
    markersize=5,
)

plt.title("Moving Window Smoothing of NDVI Time Series")
plt.ylabel("NDVI")
plt.show()

plt.figure(figsize=(12, 4)) sns.lineplot( data=wheat_ndvi, x="time", y="ndvi", label="Wheat NDVI", alpha=0.5, c="#a87000", ) sns.lineplot( data=wheat_ndvi_smoothed, x="time", y="ndvi_median", label="Wheat NDVI Smoothed", c="#a87000", marker="o", markersize=5, ) sns.lineplot( data=soybean_ndvi, x="time", y="ndvi", label="Soybean NDVI", alpha=0.5, c="#267300", ) sns.lineplot( data=soybean_ndvi_smoothed, x="time", y="ndvi_median", label="Soybean NDVI Smoothed", c="#267300", marker="o", markersize=5, ) sns.lineplot( data=corn_ndvi, x="time", y="ndvi", label="Corn NDVI", alpha=0.5, c="#ffd400" ) sns.lineplot( data=corn_ndvi_smoothed, x="time", y="ndvi_median", label="Corn NDVI Smoothed", c="#ffd400", marker="o", markersize=5, ) plt.title("Moving Window Smoothing of NDVI Time Series") plt.ylabel("NDVI") plt.show()