<jupyter_start><jupyter_text>Phenology Sentinel 2  Import Python Packages<jupyter_code>import datacube
import geopandas as gpd
from odc.ui import with_ui_cbk
import xarray as xr
import seaborn as sns
import hdstats
from deafrica_tools.bandindices import calculate_indices, dualpol_indices
from deafrica_tools.spatial import xr_rasterize
import pandas as pd<jupyter_output><empty_output><jupyter_text>Import the Datacube Configuration<jupyter_code>dc = datacube.Datacube()<jupyter_output><empty_output><jupyter_text>AoI shapefile<jupyter_code>gdf = gpd.read_file('data/sample_fields.shp')
gdf.plot()
gdf.explore()
x = (gdf.total_bounds[0] , gdf.total_bounds[2] )
y  = (gdf.total_bounds[1] , gdf.total_bounds[3] )
ds = load_ard(dc=dc,
            products=['s2_l2a'],
            x= x,
            y= y,
            time = ("2021-01-01", "2021-12-31"), # specifiy time_extent
            output_crs = "EPSG:32632",
            measurements = ["blue","green","red", "nir_1"],
            resolution = (-10,10),
            group_by = "solar_day",
            mask_pixel_quality=True,
            data_coverage = 100,
            min_gooddata=0.90,    
             )<jupyter_output><empty_output><jupyter_text>Calculate NDVI<jupyter_code>ds = calculate_indices(ds, index=['NDVI'], collection='s2')<jupyter_output><empty_output><jupyter_text>Rasterize Fields<jupyter_code>field_raster = xr_rasterize(gdf, ds, attribute_col="field_id")
field_raster.plot()
field_1<jupyter_output><empty_output><jupyter_text>Extract NDVI for one field<jupyter_code>field_1 = ds.NDVI.where(field_raster == 1)
field_1.plot(col="time",col_wrap=5,cmap='RdYlGn')
field_1.mean(['x', 'y']).plot.line('b-^', figsize=(11,4))<jupyter_output><empty_output><jupyter_text>Extract NDVI for all fields<jupyter_code>data_frame = pd.DataFrame({"time": ds.time.values})
data_frame = data_frame.set_index('time')

for i in range(1,8):
    data_frame['field_' + str(i)] = ds.NDVI.where(field_raster == i).mean(['x', 'y']).values

data_frame
sns.set_style('darkgrid')
sns.set_context("poster", font_scale = .7)

ax = data_frame.plot(figsize=[15,7], linewidth=1)<jupyter_output><empty_output><jupyter_text>Calculate phenology statistics using xr_phenology The DE Africa function xr_phenology can calculate a number of land-surface phenology statistics that together describe the characteristics of a plant’s lifecycle. The function can calculate the following statistics on either a zonal timeseries (like the one above), or on a per-pixel basis:SOS = DOY of start of seasonPOS = DOY of peak of seasonEOS = DOY of end of seasonvSOS = Value at start of seasonvPOS = Value at peak of seasonvEOS = Value at end of seasonTrough = Minimum value of seasonLOS = Length of season (DOY)AOS = Amplitude of season (in value units)ROG = Rate of greeningROS = Rate of senescence<jupyter_code>from deafrica_tools.temporal import xr_phenology

basic_pheno_stats = ['SOS','vSOS','POS','vPOS','EOS','vEOS','Trough','LOS','AOS','ROG','ROS']
method_sos = 'first'
method_eos = 'last'

stats=xr_phenology(
    ds.NDVI,
    method_sos=method_sos,
    method_eos=method_eos,
    stats=basic_pheno_stats,
    verbose=False
)
pheno_results = stats
import matplotlib.pyplot as plt
#select the year to plot

cbar_size = 0.5
phen = pheno_results

#mask again with crop-mask
phen = phen.where(field_raster != 0)

# set up figure
fig, ax = plt.subplots(nrows=5,
                       ncols=2,
                       figsize=(18, 25),
                       sharex=True,
                       sharey=True)

# set colorbar size
cbar_size = 0.7

# set aspect ratios
for a in fig.axes:
    a.set_aspect('equal')

# start of season
phen.SOS.plot(ax=ax[0, 0],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 0].set_title('Start of Season (DOY)')

phen.vSOS.plot(ax=ax[0, 1],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[0, 1].set_title('NDVI at SOS')

# peak of season
phen.POS.plot(ax=ax[1, 0],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 0].set_title('Peak of Season (DOY)')

phen.vPOS.plot(ax=ax[1, 1],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[1, 1].set_title('NDVI at POS')

# end of season
phen.EOS.plot(ax=ax[2, 0],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[2, 0].set_title('End of Season (DOY)')

phen.vEOS.plot(ax=ax[2, 1],
               cmap='YlGn',
               vmax=0.8,
               cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[2, 1].set_title('NDVI at EOS')

# Length of Season
phen.LOS.plot(ax=ax[3, 0],
              cmap='magma_r',
              vmax=300,
              vmin=0,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[3, 0].set_title('Length of Season (DOY)')

# Amplitude
phen.AOS.plot(ax=ax[3, 1],
              cmap='YlGn',
              vmax=0.8,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[3, 1].set_title('Amplitude of Season')

# rate of growth
phen.ROG.plot(ax=ax[4, 0],
              cmap='coolwarm_r',
              vmin=-0.02,
              vmax=0.02,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[4, 0].set_title('Rate of Growth')

# rate of Sensescence
phen.ROS.plot(ax=ax[4, 1],
              cmap='coolwarm_r',
              vmin=-0.02,
              vmax=0.02,
              cbar_kwargs=dict(shrink=cbar_size, label=None))
ax[4, 1].set_title('Rate of Senescence')
plt.suptitle('Phenology')
plt.tight_layout();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13155631986242478, 'fraction_numerical': 0.024591573516766982, 'mean_word_length': 3.2483564645726806, 'pattern_counts': {'":': 1, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '304328', 'license': 'MIT', 'n_tokens_mistral': 2280, 'n_tokens_neox': 2150, 'n_words': 412, 'path': '3_Data_Cube_Basics/notebooks/phenology_example/S2_phenology_1.ipynb', 'repo_name': 'eo2cube/datacube_eagles', 'size': '5815'}