<jupyter_start><jupyter_text>Introduction - Dense Timeseries Change DetectionChange detection is a *hot* topic within the Earth Observation (EO) community, particular since the release of the landsat archive and more recently the availablity of the  ESA Sentinel-1 and 2 data.There are several approaches to change detection and these can be categorised as:- Map-to-Map- Image-to-Image- Map-to-Image- Dense TimeseriesThese approaches have different advances and disadvantages and this tutorial will demonstrate the implementation of a method for each type. What is Dense Timeseries Analysis?Dense timeseries methods are a **form of image-to-image methods** in that image pixel values are being compared. However, the assumptions and properties of these methods can be quite different. For example, rather than seasonability being noise causing false changes to be found it is a single which is used to identify change. These methods are typically more **computationally complex**, and therefore **time consuming** when compared to other methods, while also **requiring a singnificant amount of data** to work. These methods will usually fit some form of model to the data and predict image pixel values.  What are the steps:This kind of analysis is complex and can be very computationally intensive - hense the very small subset of data we are using. The first step is to fit a model to each pixel in the timeseries. Then we can   1. Define Imports<jupyter_code>import os


from rsgislib.timeseries import modelfitting
import rsgislib.imagecalc
import rsgislib.imagecalc.calcindices
import rsgislib.imageutils
import rsgislib.rastergis

# Imports for data visualisation
import rsgislib.tools.plotting
# Import the matplotlib plotting module and set it to work within
# a notebook environment.
%matplotlib inline
import matplotlib
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>2. Define Input Files and Directories<jupyter_code># A JSON file which specifies the order and acquasition dates
# of the input images
sen2_imgs_lut = "05_timeseries_imgs.json"

timeseries_chg_dir = "timeseries_chg"
if not os.path.exists(timeseries_chg_dir):
    os.mkdir(timeseries_chg_dir)<jupyter_output><empty_output><jupyter_text>3. Fit the timeseries model<jupyter_code>coeffs_img = os.path.join(timeseries_chg_dir, "timeseries_mdl_coeffs.kea")

# All the input images need to have the same number of image bands
# we will however only use a subset of the bands to fit the model
# in this case bands 3 (Red), 7 (NIR) and 9 (SWIR1).
modelfitting.get_ST_model_coeffs(
    sen2_imgs_lut, coeffs_img, bands=[3, 7, 9], num_processes=1
)<jupyter_output><empty_output><jupyter_text>4. Predict an Output imageUse the model we have fitted we can now use that model to predict an output image.<jupyter_code>pred_refl_test_img = os.path.join(timeseries_chg_dir, "sen2_pred_img_test.kea")

modelfitting.predict_for_date("2017-07-15", coeffs_img, pred_refl_test_img)<jupyter_output><empty_output><jupyter_text>5. Visualise Predicted Image<jupyter_code>ref_refl_img = "../data/chg_data/timeseries_sen2/SEN2_20170715_lats77lone395_T37MEM_ORB049_utm37s_vmsk_mclds_topshad_rad_srefdem_stdsref.kea"

# Get image data to visualise
img_ref_data, img_ref_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    ref_refl_img, bands=[7, 9, 3]
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
img_ref_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(img_ref_data, no_data_val=0.0)

# Get image data to visualise
img_pred_data, img_pred_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    pred_refl_test_img, bands=[2, 3, 1]
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
img_pred_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(img_pred_data, no_data_val=0.0)


# Create the matplotlib figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 10))

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax1.imshow(img_ref_data_strch, extent=img_ref_coords)

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax2.imshow(img_pred_data_strch, extent=img_pred_coords)<jupyter_output>Image Data Size: 169 x 140
Image Data Size: 169 x 140<jupyter_text>6. Predict an Image into the Future.In this case we are predicting an image into the future, as the model is fitted to data from 2016 to 2018 but the image we are predicting will be for 2019-06-10.<jupyter_code>pred_refl_img = os.path.join(timeseries_chg_dir, "sen2_pred_img_20190610.kea")

modelfitting.predict_for_date("2019-06-10", coeffs_img, pred_refl_img)
# Get image data to visualise
img_pred_data, img_pred_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    pred_refl_img, bands=[2, 3, 1]
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
img_pred_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(img_pred_data, no_data_val=0.0)


# Create the matplotlib figure
fig, ax = plt.subplots(figsize=(10, 10))

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax.imshow(img_pred_data_strch, extent=img_ref_coords)<jupyter_output>Image Data Size: 169 x 140<jupyter_text>7. Calculate Change  7.1 Define Input Image<jupyter_code>sen2_img_20190610 = "../data/chg_data/timeseries_sen2/SEN2_20190610_lats77lone395_T37MEM_ORB049_utm37s_vmsk_mclds_topshad_rad_srefdem_stdsref.kea"<jupyter_output><empty_output><jupyter_text>7.2 Calculate NDVI<jupyter_code>pred_ndvi_img = os.path.join(timeseries_chg_dir, "sen2_pred_img_20190610_ndvi.kea")
rsgislib.imagecalc.calcindices.calc_ndvi(pred_refl_img, 1, 2, pred_ndvi_img)

sen2_ndvi_img = os.path.join(timeseries_chg_dir, "sen2_img_20190610_ndvi.kea")
rsgislib.imagecalc.calcindices.calc_ndvi(sen2_img_20190610, 3, 7, sen2_ndvi_img)<jupyter_output>Image: timeseries_chg/sen2_pred_img_20190610.kea
	 Variable 'red' is band 1
	 Variable 'nir' is band 2
New image width = 169 height = 140 bands = 1



Calculating Image Pyramids.

Image: ../data/chg_data/timeseries_sen2/SEN2_20190610_lats77lone395_T37MEM_ORB049_utm37s_vmsk_mclds_topshad_rad_srefdem_stdsref.kea
	 Variable 'red' is band 3
	 Variable 'nir' is band 7
New image width = 169 height = 140 bands = 1



Calculating Image Pyramids.<jupyter_text>7.3 Visualise NDVI Images<jupyter_code># Get image data to visualise
img_sen2_data, img_sen2_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    sen2_ndvi_img,
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
img_sen2_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(img_sen2_data, no_data_val=0.0)

# Get image data to visualise
img_pred_data, img_pred_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    pred_ndvi_img,
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
img_pred_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(img_pred_data, no_data_val=0.0)


# Create the matplotlib figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 10))

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax1.imshow(img_sen2_data_strch, extent=img_ref_coords)
ax1.set_title("Sentinel-2 NDVI")

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax2.imshow(img_pred_data_strch, extent=img_pred_coords)
ax2.set_title("Predicted NDVI")<jupyter_output>Image Data Size: 169 x 140
Image Data Size: 169 x 140<jupyter_text>7.4 Calculate NDVI Difference<jupyter_code>ndvi_diff_img = os.path.join(timeseries_chg_dir, "ndvi_diff_act_pred.kea")

band_defns = []
band_defns.append(rsgislib.imagecalc.BandDefn("ndvi_pred", pred_ndvi_img, 1))
band_defns.append(rsgislib.imagecalc.BandDefn("ndvi_act", sen2_ndvi_img, 1))

exp = "ndvi_pred-ndvi_act"

rsgislib.imagecalc.band_math(
    ndvi_diff_img, exp, "KEA", rsgislib.TYPE_32FLOAT, band_defns
)

rsgislib.imageutils.pop_img_stats(
    ndvi_diff_img, use_no_data=False, no_data_val=0, calc_pyramids=True
)<jupyter_output>Image: timeseries_chg/sen2_img_20190610_ndvi.kea
	 Variable 'ndvi_act' is band 1
Image: timeseries_chg/sen2_pred_img_20190610_ndvi.kea
	 Variable 'ndvi_pred' is band 1
New image width = 169 height = 140 bands = 1



Calculating Image Pyramids.<jupyter_text>7.5 Visualise NDVI Difference<jupyter_code># Get image data to visualise
ndvi_diff_data, ndvi_diff_coords = rsgislib.tools.plotting.get_gdal_raster_mpl_imshow(
    ndvi_diff_img,
)

# Stretch the image data so is visualises appropriately and usefully
# In this case the resulting values will be between 0-1.
ndvi_diff_data_strch = rsgislib.tools.plotting.cumulative_stretch_np_arr(ndvi_diff_data, no_data_val=0.0)


# Create the matplotlib figure
fig, ax = plt.subplots(figsize=(10, 10))

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax.imshow(ndvi_diff_data_strch, extent=ndvi_diff_coords)<jupyter_output>Image Data Size: 169 x 140<jupyter_text>7.6 Calculate Change Features<jupyter_code>change_feats_img = os.path.join(timeseries_chg_dir, "change_feats.kea")

band_defns = []
band_defns.append(rsgislib.imagecalc.BandDefn("ndvi_diff", ndvi_diff_img, 1))

exp = "ndvi_diff>0.3?1:ndvi_diff<-0.3?2:0"

rsgislib.imagecalc.band_math(
    change_feats_img, exp, "KEA", rsgislib.TYPE_8UINT, band_defns
)

rsgislib.rastergis.pop_rat_img_stats(
    change_feats_img, add_clr_tab=True, calc_pyramids=True, ignore_zero=True
)<jupyter_output>Image: timeseries_chg/ndvi_diff_act_pred.kea
	 Variable 'ndvi_diff' is band 1
New image width = 169 height = 140 bands = 1

Get Image Min and Max.

Get Image Histogram.

Adding Histogram and Colour Table to image file
Calculating Image Pyramids.<jupyter_text>7.7 Apply Colours and Names<jupyter_code>class_info_dict = dict()
class_info_dict[1] = {
    "classname": "Pos NDVI Chng",
    "red": 255,
    "green": 0,
    "blue": 0,
}
class_info_dict[2] = {
    "classname": "Neg NDVI Chng",
    "red": 0,
    "green": 0,
    "blue": 255,
}

rsgislib.rastergis.set_class_names_colours(change_feats_img, "class_names", class_info_dict)<jupyter_output><empty_output><jupyter_text>7.8 Visualise Change Features<jupyter_code># Get the image data using the get_gdal_thematic_raster_mpl_imshow function.
(
    img_msk_data,
    img_msk_coords,
    lgd_msk_patches,
) = rsgislib.tools.plotting.get_gdal_thematic_raster_mpl_imshow(
    change_feats_img,
)

# Create the matplotlib figure
fig, ax = plt.subplots(figsize=(10, 10))

# Use the imshow function to display the image data within the plot
# the extent option defines the x and y axis values.
ax.imshow(img_msk_data, extent=img_msk_coords)<jupyter_output>Image Data Size: 169 x 140
{'fraction_non_alphanumeric': 0.07510233137568963, 'fraction_numerical': 0.03488165153941983, 'mean_word_length': 6.002492211838006, 'pattern_counts': {'":': 8, '<': 50, '<?xml version=': 0, '>': 50, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '458054', 'license': 'MIT', 'n_tokens_mistral': 4338, 'n_tokens_neox': 4028, 'n_words': 1092, 'path': '06_change_detection/05_dense_timeseries_analysis.ipynb', 'repo_name': 'remotesensinginfo/rsgislib-tutorials', 'size': '11238'}