Part 14: Spectral and Temporal Analysis (2)
Chapter 2: Time Series Analysis of NDVI for Sample Points
2.1 Input Manually Sampled Points
# Define sample points for phenology analysis
sugarcane = ee.Feature(
geom=ee.Geometry.MultiPoint([[107.566527, 22.333124], [107.552301, 22.337652], [107.589958, 22.332473]]),
opt_properties={'LandCover': 'Sugarcane'}
)
forest = ee.Feature(
ee.Geometry.MultiPoint([[107.550482, 22.332677], [107.58648, 22.32705]]),
opt_properties={'LandCover': 'Eucalyptus'}
)
# Combine samples into a feature collection
samples = ee.FeatureCollection([sugarcane.buffer(30), forest.buffer(30)])2.2 Display Sample Points
# Define display styles for each land cover type
class_styles = ee.Dictionary({
'Eucalyptus': {'color': '#00FF00', 'pointSize': 5, 'pointShape': 'circle'},
'Sugarcane': {'color': '#EEB422', 'pointSize': 6, 'pointShape': 'square'}
})
# Apply styles to the samples
samples = samples.map(lambda f: f.set('style', class_styles.get(f.get('LandCover'))))
# Style the FeatureCollection according to each feature's "style" property
samples_VisCustom = samples.style(
styleProperty='style',
neighborhood=8 # maximum "pointSize" + "width" among features
)
# Display the samples on the map
Map = geemap.Map(basemap='HYBRID')
Map.addLayer(samples_VisCustom, None, 'samples_all_VisCustom')
Map.centerObject(samples, 14)
Map2.3 Get Landsat Images for the Time Range
imgs = pyeepkg.Landsat().get_readyDat_Landsat5789_C2L2SR(
bound=samples.geometry().buffer(5000).bounds(),
yearS=2018,
yearE=2021,
clip_with_bound=True
)2.4 Moving Window Smoothing Function
# Define the smoothing function
def avg_smoothing_imgs(imgs, days_window_radius, bands):
""" Smoothing the imageCollection using mean value
Args:
imgs (ee.ImageCollection): ImageCollection to be processed
days_window_radius (int): Days window to smooth the time series
bands (list): Bands to be smoothed, e.g. ['LSWI']
Returns:
ee.ImageCollection: Smoothed time series
"""
join = ee.Join.saveAll(matchesKey='images')
days_millis = days_window_radius * 24 * 60 * 60 * 1000
filter_days_differ = ee.Filter.maxDifference(
difference=days_millis,
leftField='system:time_start',
rightField='system:time_start'
)
smoothingJoin = join.apply(
primary=imgs,
secondary=imgs,
condition=filter_days_differ
)
band_types = imgs.first().select(bands).bandTypes()
new_bnames = [band + 'm' for band in bands]
def avg_imgs(image):
col = ee.ImageCollection.fromImages(image.get('images'))
return image.addBands(col.select(bands).mean().cast(band_types).rename(new_bnames))
return ee.ImageCollection(smoothingJoin).map(avg_imgs)2.5 Calculate NDVI and Apply Smoothing
# Function to add NDVI band to images
def addNDVI(image):
return image.addBands(image.normalizedDifference(['nir', 'red']).rename('NDVI'))
# Apply NDVI calculation and smoothing
with_ndvi = imgs.map(addNDVI)
with_ndvim = avg_smoothing_imgs(with_ndvi, 32, ['NDVI'])2.6 Extract Sample Pixel Values from All Images
def get_time_series(img):
img = ee.Image(img)
dat_regions = img.reduceRegions(
collection=samples,
reducer=ee.Reducer.mean(),
scale=30
)
return dat_regions.map(lambda f: f.set({'t': img.get('system:time_start')}))2.7 Iteratively Extract Sample Pixel Values
fc_dat = with_ndvim.sort('CLOUD_COVER_LAND').select(['NDVI', 'NDVIm']).map(get_time_series)
lst_fc_dat = ee.FeatureCollection(fc_dat).toList(500).map(lambda l: ee.FeatureCollection(l).toList(10))
fc_dat_ = ee.FeatureCollection(lst_fc_dat.flatten())
df_samples = geemap.ee_to_pandas(fc_dat_)
df_samples.drop('style', axis=1, inplace=True)
df_samples
df_samples.to_excel('../Data/ndvim.xlsx')2.8 Add Image Timestamp Field
# Convert milliseconds to datetime
df_samples['dt'] = pd.to_datetime(df_samples['t'], unit='ms')
df_samples.sort_values('t')2.9 Plot NDVI Time Series for Two Tree Species
import seaborn as sns
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(16, 4))
sns.set_theme(style='ticks', font_scale=1.4)
# Plot original NDVI data
g = sns.lineplot(x='dt', y='NDVI', hue='LandCover', data=df_samples, style='LandCover', markers=["o", "d"], palette=['#54b345', '#BB9727'], sizes=20)
# Modify plot parameters
g.set(ylabel='NDVI', ylim=(0, 1), xlabel='')
g.spines['right'].set_visible(False)
g.spines['top'].set_visible(False)
ax.fill_between(['2020-01-01', '2020-05-01'], 0, 1, facecolor='#ff2500', alpha=0.1)
g.get_legend().set_title('')
g.grid(ls='-.')
fig.savefig('../Figures/Pheno.tif', dpi=400, bbox_inches='tight')2.10 Plot Original and Smoothed NDVI Time Series for Two Tree Species
from matplotlib.lines import Line2D
fig, ax = plt.subplots(figsize=(16, 4))
sns.set_theme(style='ticks', font_scale=1.5)
# Plot smoothed NDVI time series
g = sns.lineplot(x='dt', y='NDVIm', hue='LandCover', data=df_samples, style='LandCover', markers=["o", "d"], palette=['#F08080', '#A2CD5A'], sizes=10)
# Plot original NDVI time series
g = sns.scatterplot(x='dt', y='NDVI', hue='LandCover', data=df_samples, style='LandCover', markers=["o", "d"], palette=['#FF0000', '#008B00'], sizes=20)
# Modify plot parameters
g.set(ylabel='NDVI', ylim=(0, 1), xlabel='')
g.spines['right'].set_visible(False)
g.spines['top'].set_visible(False)
ax.fill_between(['2018-01-01', '2018-06-01'], 0, 1, facecolor='#87CEFA', alpha=0.2)
ax.fill_between(['2019-01-01', '2019-06-15'], 0, 1, facecolor='#87CEFA', alpha=0.2)
ax.fill_between(['2020-01-01', '2020-05-01'], 0, 1, facecolor='#87CEFA', alpha=0.2)
g.grid(ls='-.')
# Custom legend
legend_elements = [
Line2D([], [], marker='o', color='#FF0000', label='$Sugarcane_{original-value}$', ls=''),
Line2D([], [], marker='d', color='#008B00', label='$Eucalyptus_{original-value}$', ls=''),
Line2D([], [], marker='o', color='#F08080', label='$Sugarcane_{moving-avg}$', ls='-'),
Line2D([], [], marker='d', color='#A2CD5A', label='$Eucalyptus_{moving-avg}$', ls='-'),
]
g.legend(ncol=4, handles=legend_elements, prop={'size': 14})
fig.savefig('../Figures/Pheno_with_Fit.tif', dpi=400, bbox_inches='tight')Last updated