<jupyter_start><jupyter_text>Training of a sequential neural network for soil moisture prediction by <NAME>, University of Tor Vergata, Rome, Italy   Watch the video tutorial<jupyter_code>from IPython.display import HTML
HTML('<div align="center"><iframe src="https://player.vimeo.com/video/631902326?h=e76118f979" width="640" height="360" frameborder="0" allow="autoplay; fullscreen; picture-in-picture" allowfullscreen align="middle"></iframe></div>')<jupyter_output><empty_output><jupyter_text>Introduction This workflow will guide you through the different steps to train a `sequential neural network model` in order to predict soil moisture based on SAR backscatter and incidence angle information. The model is trained based on VH and VV backscatter information retrieved from an implementation of the [Integral Equation Model modified by Baghdadi (IEM-B)](https://www.mdpi.com/2073-4441/9/1/38).The IEM-B is a well-known inversion model which can simulate the backscattering and polarisation of SAR data based on ground-truth soil moisture values. The input data to the model were soil moisture values from the [International Soil Moisture Network (ISMN)](https://hess.copernicus.org/preprints/hess-2021-2/) and SAR incidence angles from 30 to 32.The `sequential neural network` learns how to invert from SAR backscatter and incidence angle information to predict soil moisture content. Backscatter and incidence angle information from a Sentinel-1 image and the associated soil moisture value from an ISMN station near Bordeaux, France (44.789900 N, -0.576500 W) serve for model validation.  Machine-Learning Algorithm The workflow trains a [sequential neural network with Keras](https://keras.io/api/models/sequential/), which is a linear stack of layers. There are two ways to build Keras models:* in a `sequential` way and* a `functional` way.The `sequential` API allows you to create models in a sequential way, which means you add one layer after the other. A sequential model is applicable for most problems. Limitations of `sequential` models are that they do not allow you to create models that share layers or that have multiple inputs or outputs.  Data This workflow makes use of the .txt file [MOOC_IEM_Data](./MOOC_IEM_Data.txt). The file contains of a total of 268 rows and four columns. The first 264 rows contain the IEM-B model outputs (`VV- and VH-backscatter` values) of a sequence of `soil moisture values` ranging from 0.01 to 0.81 with a step of 0.025 and associated `incidence angles` ([30,32]). The first 264 rows are used to train the sequential model. The last four rows contain `VV- and VH-backscatter values, incidence angles` from Sentinel-1 images and associated measured `soil moisture values` from a location near Bordeaux, France (44.789900 N, -0.576500 W) taken from the International Soil Moisture Network (ISMN).The columns are as follows: * Column 1 - `SAR VH backscattering values in dB`* Column 2 - `SAR VV backscattering values in dB`* Column 3 - `incidence angles in degrees`* Column 4 - `volumetric soil moisture content`  Further resources * [Evaluation of the Oh, Dubois and IEM Backscatter Models Using a Large Dataset of SAR Data and Experimental Soil Measurements](https://www.mdpi.com/2073-4441/9/1/38)* [International Soil Moisture Network (ISMN)](https://hess.copernicus.org/preprints/hess-2021-2/)   Notebook outline* [1 - Load and prepare the input and output data for model training](load_prepare)* [2 - Define and compile a sequential neural network model with Keras](define_seq_model)* [3 - Training (fitting) of the sequential model](sequential_fitting)* [4 - Evaluate model performance](model_evaluation)* [5 - Predicting soil moisture with the trained sequential neural network](soil_moisture_prediction)  Import libraries<jupyter_code>## BEGIN S3FS IMPORT SNIPPET ##
import os, sys
s3_home =  os.getcwd()
try: sys.path.remove(s3_home) # REMOVE THE S3 ROOT FROM THE $PATH
except Exception: pass
current_dir = os.getcwd()
os.chdir('/home/jovyan') # TEMPORARILY MOVE TO ANOTHER DIRECTORY

# BEGIN IMPORTS #
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import LeakyReLU
from tensorflow.keras.optimizers import SGD

from sklearn import preprocessing

from matplotlib import pyplot
import datetime

import numpy as np
import scipy
# END IMPORTS #

os.chdir(current_dir) # GO BACK TO YOUR PREVIOUS DIRECTORY
sys.path.append(s3_home) # RESTORE THE S3 ROOT IN THE $PATH
## END S3FS IMPORT SNIPPET ##<jupyter_output><empty_output><jupyter_text>1. Load and prepare the input and output data for model training The first step is to load the .txt file  [MOOC_IEM_Data.txt](./MOOC_IEM_Data.txt). The file contains 268 rows and four columns. Each row contains a different sample. The dataset contains `SAR VH backscatter in dB`, `SAR VV backscatter in dB`, `incidence angle in degrees` and `volumetric soil moisture content`. The first 264 rows indicate the IEM-B model outputs, while the last four rows hold information from Sentinel-1 images and soil moisture values from a station of the ISMN.With the numpy function `loadtxt`, you can load a text file, by indicating the delimiter used in the file.<jupyter_code>dataset = np.loadtxt('./MOOC_IEM_Data.txt', delimiter=',')
dataset<jupyter_output><empty_output><jupyter_text>We split the loaded dataset in `input (X)` and `output (y)` variables. From the dataset file, the first three columns (`SAR VH backscattering values in dB`, `SAR VV backscattering values in dB` and `incidence angle in degrees`) serve as input (X) values.The fourth column of the dataset file contains values of `volumetric soil moisture content`, which serve as output (y) variables.<jupyter_code>X = dataset[:, 0:3]
y = dataset[:, 3]

X, y<jupyter_output><empty_output><jupyter_text>Normalisation of input variables `(Feature) Normalisation` is a technique often applied in Machine Learning. The goal of normalisation is to change the values of the numeric columns to a common scale, without distorting differences in the ranges of values. For Machine Learning, not every dataset requires normalisation. It is only required when features have different ranges or units. This is the case with the present input data, which are `VV backscatters (in dB)`, `VH backscatters (in dB)` and `incidence angles in degrees`.Scikitlearn's preprocessing library offers the function `minmax_scale()`, which normalizes the input values to a range from 0 to 1. This normalisation process is known as `unity-based normalisation` and is a very common normalisation process in Machine Learning.The normalised values range then from 0 to 1.<jupyter_code>preprocessed_X = preprocessing.minmax_scale(X, feature_range=(0, 1))
preprocessed_X<jupyter_output><empty_output><jupyter_text>Split input data into training and test data In the next step, you want to split the input data into training and test data. As a first step, let us define a variable called `iem_samples` which holds the number of `integral equation model (IEM)` samples (264). This variable helps you to split the input (X) data into training and test data.<jupyter_code>iem_samples = dataset.shape[0]-4
iem_samples<jupyter_output><empty_output><jupyter_text>With the variable defined in the previous step, you can split the normalised input (X) and output (y) data into training and test samples. You want to use the 264 IEM_B samples for training and the four samples based on Sentinel-1 images as test data.<jupyter_code>test_X_S1 = preprocessed_X[iem_samples:]
preprocessed_X_iem = preprocessed_X[:iem_samples]

test_y_S1 = y[iem_samples:]
y_iem = y[:iem_samples]<jupyter_output><empty_output><jupyter_text>In a next step, you want to `random shuffle` the data. By `random shuffling` the data, you de-correlate potential correlation or a human-induced ordering of the input (X) and output (y) data.In order to random shuffle the data, you first create a list of integers with 264 entries with `np.arange()`, indicating the index values of a Python list. Then, you randomly shuffle the list of index values with `np.random.shuffle()`.The randomly shuffled list of index values can then be applied to the training input (`preprocessed_X_iem`) and training output (`y_iem`) data.<jupyter_code>s = np.arange(iem_samples)
np.random.shuffle(s)
preprocessed_X_shuffled = preprocessed_X_iem[s]
y_shuffled = y_iem[s]<jupyter_output><empty_output><jupyter_text>Split training data into training and test data  The next step splits the training data into two subsets: * `train_X` / `train_y` --> 70% of the training data* `test_X` / `test_y` --> 30% of the training data. The first 185 entries (70% of 264) are used for training the model and the remaining data are used for testing. You generate training and testing subsets originating from `preprocessed_X_shuffled` (input) and `y_shuffled` (output) data.<jupyter_code>n_train_perc = 0.7
n_train = round(preprocessed_X_shuffled.shape[0] * n_train_perc)
n_train
train_X  = preprocessed_X_shuffled[:n_train, :]
train_y  = y_shuffled[:n_train] 

test_X = preprocessed_X_shuffled[n_train:, :]
test_y = y_shuffled[n_train:]

n_train, train_X, test_X<jupyter_output><empty_output><jupyter_text>2. Define and compile a sequential neural network model with Keras Now you can start defining the sequential neural network model with `keras.Sequential()`.<jupyter_code>model = keras.Sequential()<jupyter_output><empty_output><jupyter_text>The next step is to add a defined number of layers to the sequential model. The procedure is to define three `dense layers`:* Layer 1: `dense layer` with 8 neurons and 3 input dimensions, because the input data has three variables* Layer 2: `dense layer` with 3 neurons* Layer 3: `dense layer` with 1 neuron, because we want to have one output, soil moisture Between the layers, we provide the network with `non-linear properties` by applying an activation function, e.g. `LeakyReLU` or `sigmod` to each layer.<jupyter_code># Layer 1
model.add(Dense(8, input_dim=3))
model.add(LeakyReLU(alpha=0.1))

# Layer 2
model.add(Dense(3))
model.add(LeakyReLU(alpha=0.1))

# Layer 3
model.add(Dense(1, activation='sigmoid'))<jupyter_output><empty_output><jupyter_text>Once the model is defined, you can compile (configure) the model with `model.compile()` and you can define the following hyperparameters:* `loss='mean_absolute_error'` - MAE is one of many loss options and calculates the distance to the output values which is used to minimize the coefficients* `optimizer=SGD(lr=0.01, momentum=0.9)` - Optimizers are algorithms the network learns from and SGD stands for gradient descent optimizer* `metrics=['mse']` - is used to evaluate how the model is performing**Note related to optimizers:***During the learning process the network is performing a forward pass which provides the network with an amount to which the coefficients shall be corrected - lr=learning rate. Learning rate and momentum are hyperparameters that can be changed*<jupyter_code>model.compile(loss='mean_absolute_error', 
              optimizer=SGD(lr=0.01, momentum=0.9), 
              metrics=['mse'])<jupyter_output><empty_output><jupyter_text>3. Training (Fitting) of the sequential model Now, the data is prepared and the model has been initiated and compiled. The next step is now to train (fit) the model with the defined training data. Before you can start with the model training, it is helpful to define some callbacks that are useful to track model performance and outcomes.We define three callbacks: (i) a `tensorboard`, (ii) an `EarlyStopping` and (iii) a `ModelCheckpoint`. Defining a `tensorboard` is optional, but it is a Tensorflow framework that enables you to easily visualize specific model metrics. You can create a tensorboard with `keras.callbacks.TensorBoard()`.<jupyter_code>log_dir = "../../../../logs/fit/" + datetime.datetime.now().strftime("%Y%m%d-%H%M%S")
tensorboard_callback = keras.callbacks.TensorBoard(log_dir=log_dir, histogram_freq=1)<jupyter_output><empty_output><jupyter_text>Another useful `callback` metric may be `EarlyStopping`, which stops the training if the validation loss (`val_loss`) reaches a minimum (`min`), but waits (`patience`) for 50 epochs.You can define an `EarlyStopping` callback with `keras.callbacks.EarlyStopping()`.<jupyter_code>es = keras.callbacks.EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=50)<jupyter_output><empty_output><jupyter_text>Another common callback used is `ModelCheckpoint`, which safes the keras model. During training the `Validation Mean-Squared-Error (MSE)` is monitored and if it reaches a minimum, the model is safed and only the best model fit is kept.You can define a `ModelCheckpoint` callback with `keras.callbacks.ModelCheckpoint()`.<jupyter_code>mc = keras.callbacks.ModelCheckpoint('../../../../best_model.h5', monitor='val_mse', mode='min', verbose=1, save_best_only=True)<jupyter_output><empty_output><jupyter_text>The three callbacks can now be used during the training (fitting) process of the model.You can train (fit) a model in Keras with `model.fit()` and you need to specify the following parameters:* `input (X)` and `output (y)` data: here specify the input and output data of your model* `validation_data`: here we enter the test data subsets test_X and test_y and our model outputs are validated against these validation data after each epoch (training cycle)* `epochs`: number of training cycles* `batch_size`: defines the size of a training data subset (e.g. 10 samples) after which the weights of the network are updated* `callbacks`: define here the callbacks you would like to make use of during the training process **Note:***The process where the weights are updated based on the batch_size is also called 'gradient update'.* The output of the training process is a `history` object and for this reason, the output object has the name `history`.<jupyter_code>history = model.fit(train_X, train_y, 
                    validation_data=(test_X, test_y), 
                    epochs=100, 
                    batch_size=10,
                    callbacks=[tensorboard_callback, es, mc])
history<jupyter_output><empty_output><jupyter_text>4. Evaluate model performance After the training process, you want to evaluate the performance of the model. Since we defined the `ModelCheckpoint` callback, you can load the output of the `ModelCheckPoint` with `keras.models.load_model()`. Note that the output is different to the `history` object during model training.<jupyter_code>model = keras.models.load_model('./best_model.h5')
model<jupyter_output><empty_output><jupyter_text>With `model.evaluate()` you can evaluate the model predictions based on the test data subsets for X (`test_X`) and y (`test_y`) we have created before. The output returns the mean-squared-error (MSE) between predictions and test data. The aim is to optimize (minimize) the `MSE`<jupyter_code>_, test_mse = model.evaluate(test_X, test_y)

print('Mse: %.4f' % test_mse)<jupyter_output><empty_output><jupyter_text>Another useful evaluation is to plot the `loss` and `validation_loss` metrics during training:* `loss`: mean absolute error during training (blue line)* `validation_loss`: mean absolute error while validating the training outputs with test data (orange line)Both lines (blue and orange) show a decreasing trend during the training process, which indicates a successful training process. The decreasing trend of the orange line shows also that the model is capable of decreasing the mean absolute error when validated with test data.<jupyter_code>pyplot.plot(history.history['loss'], label='MAE - training')
pyplot.plot(history.history['val_loss'], label='MAE - validation with test data')
pyplot.legend()<jupyter_output><empty_output><jupyter_text>5. Predicting soil moisture with the trained sequential neural network As a final step, you can now use the trained sequential neural network and predict soil moisture, based on the test input `test_X_S1`, which has `SAR VH and VV backscattering values in db` and `incidence angle in degrees` of four Sentinel-1 samples.The code below goes through the four samples of the test input variable and based on the trained model, predicts the soil moisture based on the Sentinel-1 samples. The prediction is done with the function `model.predict()`. The predicted soil moisture values are stored in a list (`predicted`).At the end, the input and output soil moisture values are printed together with the predicted soil moisture.<jupyter_code>predicted = []
for index, value in enumerate(test_X_S1):
    pred = np.expand_dims(value, axis=0)
    predicted.append(model.predict(pred))
    
    print("Input values: " + str(value) + ",  Output soil moisture: " + str(test_y_S1[index]))
    print("Predicted soil moisture: " + str(model.predict(pred)) + "\n ")<jupyter_output><empty_output><jupyter_text>Evaluate model predictions A common tool to evaluate the accuracy of predicted values versus validation values is the visualisation of both as a scatter plot. Such a scatter plot shows the `original vs the predicted` soil moisture values as scatters. It a great graphical representation showing whether the trained sequential neural network model has a tendency to over- or underpredict `soil moisture information`. The result below shows you that the trained model has a tendency to over-predict soil moisture values.**Note**: for more accurate results, you would need more than only four point locations to evaluate the predicted values against. Further, you can modify the hyperparameters in order to improve the model's accuracy.<jupyter_code>fig, ax = pyplot.subplots()
ax.scatter(test_y_S1, predicted)
ax.plot([test_y.min(), test_y.max()], 'k--', lw=4)
ax.set_xlabel('Original')
ax.set_ylabel('Predicted')<jupyter_output><empty_output><jupyter_text>As a final step, you can calculate the `Pearson correlation coefficient` between original (to-predict) and predicted soil moisture values. The coefficient value can range between `[-1,1]` and is a measure about the strength of the correlation between original and predicted soil moisture.<jupyter_code>test_y_nparray = np.array(test_y_S1)
test_y_flattened = test_y_nparray.flatten()
predicted_nparray = np.array(predicted)
predicted_flattened = predicted_nparray.flatten()

corr, _ = scipy.stats.pearsonr(test_y_flattened, predicted_flattened)
print('Pearson correlation: %.3f' % corr)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08098305635251449, 'fraction_numerical': 0.012991934174200184, 'mean_word_length': 5.52560932532674, 'pattern_counts': {'":': 0, '<': 94, '<?xml version=': 0, '>': 96, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, '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': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '356863', 'license': 'MIT', 'n_tokens_mistral': 5449, 'n_tokens_neox': 5108, 'n_words': 2239, 'path': '3_land/3D1_soil_moisture_estimation/3D_soil_moisture_estimation.ipynb', 'repo_name': 'AIAstronomy/ai4EM_MOOC', 'size': '18489'}