Harmonic regression

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

In [ ]:

# !pip install -U geeagri

# !pip install -U geeagri

Import libraries¶

In [ ]:

import ee
import geemap
from geeagri.analysis import HarmonicRegression
from geeagri.preprocessing import Sentinel2CloudMask
from geeagri.extract import extract_timeseries_to_point

import matplotlib.pyplot as plt
import seaborn as sns

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

import ee import geemap from geeagri.analysis import HarmonicRegression from geeagri.preprocessing import Sentinel2CloudMask from geeagri.extract import extract_timeseries_to_point import matplotlib.pyplot as plt import seaborn as sns 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 Harmonics¶

In [ ]:

# Fit a 3rd order harmonic
harmonics = HarmonicRegression(
    image_collection=ndvi_col, ref_date="2020-01-01", band_name="ndvi", order=3, omega=1
)

# Get the harmonic coefficients
harmonic_coeffs = harmonics.get_harmonic_coeffs()

# Get the phase amplitude in an HSV image
phase_amplitude_hsv = harmonics.get_phase_amplitude(
    harmonic_coeffs,
    cos_band="ndvi_cos_1",
    sin_band="ndvi_sin_1",
    hsv=True,
    stretch_factor=5,
)
Map.addLayer(phase_amplitude_hsv, {}, "Phase (hue), Amplitude (sat), NDVI (val)")

# Get the raw phase and amplitude image
phase_amplitude = harmonics.get_phase_amplitude(
    harmonic_coeffs,
    cos_band="ndvi_cos_1",
    sin_band="ndvi_sin_1",
    hsv=False,
)
phase_vis = {
    "min": -3.10,
    "max": 3.10,
    "palette": [
        "#d73027",
        "#f46d43",
        "#fdae61",
        "#fee08b",
        "#d9ef8b",
        "#a6d96a",
        "#66bd63",
        "#1a9850",
    ],
}
amplitude_vis = {
    "min": 0,
    "max": 0.4,
    "palette": [
        "#ffffe5",
        "#f7fcb9",
        "#d9f0a3",
        "#addd8e",
        "#78c679",
        "#41ab5d",
        "#238443",
        "#005a32",
    ],
}

Map.addLayer(phase_amplitude.select("phase"), phase_vis, "Phase (-π to π)")
Map.addLayer(phase_amplitude.select("amplitude"), amplitude_vis, "Amplitude")

# Get fitted harmonics
fitted_harmonics = harmonics.get_fitted_harmonics(harmonic_coeffs)

# Fit a 3rd order harmonic harmonics = HarmonicRegression( image_collection=ndvi_col, ref_date="2020-01-01", band_name="ndvi", order=3, omega=1 ) # Get the harmonic coefficients harmonic_coeffs = harmonics.get_harmonic_coeffs() # Get the phase amplitude in an HSV image phase_amplitude_hsv = harmonics.get_phase_amplitude( harmonic_coeffs, cos_band="ndvi_cos_1", sin_band="ndvi_sin_1", hsv=True, stretch_factor=5, ) Map.addLayer(phase_amplitude_hsv, {}, "Phase (hue), Amplitude (sat), NDVI (val)") # Get the raw phase and amplitude image phase_amplitude = harmonics.get_phase_amplitude( harmonic_coeffs, cos_band="ndvi_cos_1", sin_band="ndvi_sin_1", hsv=False, ) phase_vis = { "min": -3.10, "max": 3.10, "palette": [ "#d73027", "#f46d43", "#fdae61", "#fee08b", "#d9ef8b", "#a6d96a", "#66bd63", "#1a9850", ], } amplitude_vis = { "min": 0, "max": 0.4, "palette": [ "#ffffe5", "#f7fcb9", "#d9f0a3", "#addd8e", "#78c679", "#41ab5d", "#238443", "#005a32", ], } Map.addLayer(phase_amplitude.select("phase"), phase_vis, "Phase (-π to π)") Map.addLayer(phase_amplitude.select("amplitude"), amplitude_vis, "Amplitude") # Get fitted harmonics fitted_harmonics = harmonics.get_fitted_harmonics(harmonic_coeffs)

Plot the fitted harmonics 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_fitted = extract_timeseries_to_point(
    lat=wheat_loc[1], lon=wheat_loc[0], image_collection=fitted_harmonics, scale=10
)
soybean_ndvi_fitted = extract_timeseries_to_point(
    lat=soybean_loc[1], lon=soybean_loc[0], image_collection=fitted_harmonics, scale=10
)
corn_ndvi_fitted = extract_timeseries_to_point(
    lat=corn_loc[1], lon=corn_loc[0], image_collection=fitted_harmonics, 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_fitted = extract_timeseries_to_point( lat=wheat_loc[1], lon=wheat_loc[0], image_collection=fitted_harmonics, scale=10 ) soybean_ndvi_fitted = extract_timeseries_to_point( lat=soybean_loc[1], lon=soybean_loc[0], image_collection=fitted_harmonics, scale=10 ) corn_ndvi_fitted = extract_timeseries_to_point( lat=corn_loc[1], lon=corn_loc[0], image_collection=fitted_harmonics, 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_fitted,
    x="time",
    y="fitted",
    label="Wheat NDVI Fitted",
    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_fitted,
    x="time",
    y="fitted",
    label="Soybean NDVI Fitted",
    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_fitted,
    x="time",
    y="fitted",
    label="Corn NDVI Fitted",
    c="#ffd400",
    marker="o",
    markersize=5,
)

plt.title("3rd-Order Harmonic Fit 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_fitted, x="time", y="fitted", label="Wheat NDVI Fitted", 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_fitted, x="time", y="fitted", label="Soybean NDVI Fitted", 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_fitted, x="time", y="fitted", label="Corn NDVI Fitted", c="#ffd400", marker="o", markersize=5, ) plt.title("3rd-Order Harmonic Fit of NDVI Time Series") plt.ylabel("NDVI") plt.show()