Moving window smoothing
Uncomment the following line to install the latest version of geeagri if needed.
In [ ]:
Copied!
# !pip install -U geeagri
# !pip install -U geeagri
In [ ]:
Copied!
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 [ ]:
Copied!
# 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 [ ]:
Copied!
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 [ ]:
Copied!
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 [ ]:
Copied!
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 [ ]:
Copied!
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 [ ]:
Copied!
# 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 [ ]:
Copied!
# 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 [ ]:
Copied!
# 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 [ ]:
Copied!
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 [ ]:
Copied!
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 [ ]:
Copied!
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()