<jupyter_start><jupyter_text>Regression Script LOO CV Partial Least Squares Regression <NAME> 2020-08-05 V0.2Quellen: - https://nirpyresearch.com/partial-least-squares-regression-python/variable importance: https://github.com/scikit-learn/scikit-learn/issues/7050<jupyter_code># Requried Packages
from osgeo import gdal, ogr, gdal_array # I/O image data
import numpy as np # math and array handling
import matplotlib.pyplot as plt # plot figures
import pandas as pd # handling large data as table sheets
from joblib import dump, load
from operator import itemgetter

from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import cross_val_predict

import datetime

import seaborn as sns

# Tell GDAL to throw Python exceptions, and register all drivers
gdal.UseExceptions()
gdal.AllRegister()

settings_sns = {'axes.facecolor': 'white',
                'axes.edgecolor': '0',
                'axes.grid': True,
                'axes.axisbelow': True,
                'axes.labelcolor': '.15',
                'figure.facecolor': 'white',
                'grid.color': '.8',
                'grid.linestyle': '--',
                'text.color': '0',
                'xtick.color': '0',
                'ytick.color': '0',
                'xtick.direction': 'in',
                'ytick.direction': 'in',
                'lines.solid_capstyle': 'round',
                'patch.edgecolor': 'w', 
                'patch.force_edgecolor': True,
                'image.cmap': 'Greys',
                'font.family': ['serif'],
                'font.sans-serif': ['Arial', 'Liberation Sans', 'DejaVu Sans', 'Bitstream Vera Sans', 'sans-serif'],
                'xtick.bottom': True,
                'xtick.top': True,
                'ytick.left': True,
                'ytick.right': True,
                'axes.spines.left': True,
                'axes.spines.bottom': True,
                'axes.spines.right': True,
                'axes.spines.top': True}
# data input

# the remote sensing image you want to classify
img_RS = 'D:/OwnCloud/DissBeyer/tibass/bal/images/007_stack/2017_09_08_Bal_stack_25_bands_40cm.tif'

# training and validation
#    as image in the same extand as your remote sensing image
#    no data pixels = 0 or negative
#    class pixels > 0 and as integer
field = 'D:/OwnCloud/DissBeyer/tibass/bal/vector/2020_09_08_Balje_Messwerte.shp'


# what is the attributes name of your classes in the shape file (field name of the classes)?
attribute = 'PAR'


# save path, predicted image
prediction_map = 'D:/OwnCloud/DissBeyer/tibass/bal/results/2020_08_24_BAL_PLSR_'+attribute+'.tif'

comp_plot = 'D:/OwnCloud/DissBeyer/tibass/bal/results/2020_08_24_BAL_PLSR_'+attribute+'_components.pdf'

validation_plot = 'D:/OwnCloud/DissBeyer/tibass/bal/results/2020_08_24_BAL_PLSR_'+attribute+'_regressionplot.pdf'

results_txt = 'D:/OwnCloud/DissBeyer/tibass/bal/results/2020_08_24_BAL_PLSR_'+attribute+'.txt'

# save path, trained model
#save_model = 'D:\\OwnCloud\\WetScapes\\2020_02_21_Regression_Biomass\\Results\\PLSR_model.joblib'

# laod training data and show all shape attributes

#model_dataset = gdal.Open(model_raster_fname)
shape_dataset = ogr.Open(field)
shape_layer = shape_dataset.GetLayer()

# extract the names of all attributes (fieldnames) in the shape file
attributes = []
ldefn = shape_layer.GetLayerDefn()
for n in range(ldefn.GetFieldCount()):
    fdefn = ldefn.GetFieldDefn(n)
    attributes.append(fdefn.name)
    
# print the attributes
print('Available attributes in the shape file are: {}'.format(attributes))
# prepare results text file:

print('PLSR - LOOCV Regression', file=open(results_txt, "a"))
print('Processing: {}'.format(datetime.datetime.now()), file=open(results_txt, "a"))
print('-------------------------------------------------', file=open(results_txt, "a"))
print('PATHS:', file=open(results_txt, "a"))
print('Image: {}'.format(img_RS), file=open(results_txt, "a"))
print('samples shape: {}'.format(field) , file=open(results_txt, "a"))
print('      choosen attribute: {}'.format(attribute) , file=open(results_txt, "a"))
print('predicted regression image: {}'.format(prediction_map) , file=open(results_txt, "a"))
print('Report text file: {}'.format(results_txt) , file=open(results_txt, "a"))
print('-------------------------------------------------', file=open(results_txt, "a"))

# Data preparation

# load image data

img_ds = gdal.Open(img_RS, gdal.GA_ReadOnly)

img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount),
               gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType))
for b in range(img.shape[2]):
    img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray()
row = img_ds.RasterYSize
col = img_ds.RasterXSize
band_number = img_ds.RasterCount

print('Image extent: {} x {} (row x col)'.format(row, col))
print('Number of Bands: {}'.format(band_number))


print('Image extent: {} x {} (row x col)'.format(row, col), file=open(results_txt, "a"))
print('Number of Bands: {}'.format(band_number), file=open(results_txt, "a"))
print('---------------------------------------', file=open(results_txt, "a"))
print('TRAINING', file=open(results_txt, "a"))

# Data preparation
# load training data from shape file

#model_dataset = gdal.Open(model_raster_fname)
shape_dataset = ogr.Open(field)
shape_layer = shape_dataset.GetLayer()
mem_drv = gdal.GetDriverByName('MEM')
mem_raster = mem_drv.Create('',img_ds.RasterXSize,img_ds.RasterYSize,1,gdal.GDT_Float32)
mem_raster.SetProjection(img_ds.GetProjection())
mem_raster.SetGeoTransform(img_ds.GetGeoTransform())
mem_band = mem_raster.GetRasterBand(1)
mem_band.Fill(0)
mem_band.SetNoDataValue(0)

att_ = 'ATTRIBUTE='+attribute
# http://gdal.org/gdal__alg_8h.html#adfe5e5d287d6c184aab03acbfa567cb1
# http://gis.stackexchange.com/questions/31568/gdal-rasterizelayer-doesnt-burn-all-polygons-to-raster
err = gdal.RasterizeLayer(mem_raster, [1], shape_layer, None, None, [1],  [att_,"ALL_TOUCHED=TRUE"])
assert err == gdal.CE_None

roi = mem_raster.ReadAsArray()

# Display image and Training data

plt.imshow(img[:, :, 0], cmap=plt.cm.Greys_r)
roi_positions = np.where(roi > 0)
plt.scatter(roi_positions[1], roi_positions[0], marker='x', c='r' )
plt.title('first RS band and sample points')


plt.show()

# Number of training pixels:
n_samples = (roi > 0).sum()
print('We have {n} training samples'.format(n=n_samples))# Subset the image dataset with the training image = X
print('{n} training samples'.format(n=n_samples), file=open(results_txt, "a"))
# Mask the classes on the training dataset = y
# These will have n_samples rows
X = img[roi > 0, :]
y = roi[roi > 0]
features = pd.DataFrame(X)

band_names = []
for i in range(X.shape[1]):
#for i in range(0,2500):
    nband = "Band_"+ str(i+1)
    band_names.append(nband)

features.columns = band_names

print('The shape of our features is:', features.shape)
print('The number of Spectra is:', features.shape[0])
print('The number of bands is:', features.shape[1])

features['value'] = y

features.head()

# Labels are the values we want to predict
labels = np.array(features['value'])

# Remove the labels from the features
# axis 1 refers to the columns
features = features.drop('value', axis = 1)

# Saving feature names for later use
feature_list = list(features.columns)

# Convert to numpy array
features = np.array(features)

print('Training Features Shape: ', features.shape)
print('Training Labels Shape: ', labels.shape)
mse = []
component = np.arange(1, features.shape[1])

for i in component:
    pls = PLSRegression(n_components=i)

    # Cross-validation
    y_cv = cross_val_predict(pls, X, y, cv=10)

    mse.append(mean_squared_error(y, y_cv))

   # comp = 100*(i+1)/40


    # Calculate and print the position of minimum in MSE
msemin = np.argmin(mse)
suggested_comp = msemin+1
print("Suggested number of components: ", suggested_comp)
print("Suggested number of components: {}".format(suggested_comp), file=open(results_txt, "a"))

with plt.style.context(('ggplot')):
    sns.set(rc=settings_sns, font_scale=1.5)
    plt.plot(component, np.array(mse), '-v', color = 'blue', mfc='blue')
    plt.plot(component[msemin], np.array(mse)[msemin], 'P', ms=10, mfc='red')
    plt.xlabel('Anzahl von PLS-Komponenten')
    plt.ylabel('MSE')
    plt.title('PLSR MSE vs. Komponenten')
    plt.xlim(left=-1)
    
    plt.savefig(comp_plot,bbox_inches='tight')

plt.show()
# Define PLS object with optimal number of components
pls_opt = PLSRegression(n_components=msemin+1)

# Fir to the entire dataset
pls_opt.fit(X, y)
y_c = pls_opt.predict(X)

# Cross-validation
y_cv = cross_val_predict(pls_opt, X, y, cv=X.shape[0])

# Calculate scores for calibration and cross-validation
score_c = r2_score(y, y_c)
score_cv = r2_score(y, y_cv)

# Calculate mean squared error for calibration and cross validation
mse_c = mean_squared_error(y, y_c)
mse_cv = mean_squared_error(y, y_cv)


print('R2 calib: %5.3f'  % score_c)
print('R2 LOOCV: %5.3f'  % score_cv)
print('MSE calib: %5.3f' % mse_c)
print('MSE LOOCV: %5.3f' % mse_cv)

print("-----------------------", file=open(results_txt, "a"))
print('R2 calib: %5.3f'  % score_c, file=open(results_txt, "a"))
print('R2 LOOCV: %5.3f'  % score_cv, file=open(results_txt, "a"))
print('MSE calib: %5.3f' % mse_c, file=open(results_txt, "a"))
print('MSE LOOCV: %5.3f' % mse_cv, file=open(results_txt, "a"))

print("-----------------------", file=open(results_txt, "a"))
print("band importance", file=open(results_txt, "a"))

# function for variable importance
def vip(model):
    t = model.x_scores_
    w = model.x_weights_
    q = model.y_loadings_
    p, h = w.shape
    vips = np.zeros((p,))
    s = np.diag(t.T @ t @ q.T @ q).reshape(h, -1)
    total_s = np.sum(s)
    for i in range(p):
        weight = np.array([ (w[i,j] / np.linalg.norm(w[:,j]))**2 for j in range(h) ])
        vips[i] = np.sqrt(p*(s.T @ weight)/total_s)
    return vips

importance = vip(pls_opt)

imp = {}
for i in range(features.shape[1]):
    print('Band {}: {}'.format(i+1, importance[i]))
    print('Band {}: {}'.format(i+1, importance[i]), file=open(results_txt, "a"))
    imp['Band{}'.format(i+1)] = importance[i]



print("-----------------------", file=open(results_txt, "a"))
print("sorted band importance", file=open(results_txt, "a"))
sorted_imp = dict(sorted(imp.items(), key=itemgetter(1), reverse=True))

for i in sorted_imp:
    print('{}: {}'.format(i, sorted_imp[i]))

for i in sorted_imp:
    print('{}: {}'.format(i, sorted_imp[i]), file=open(results_txt, "a"))
print("-----------------------", file=open(results_txt, "a"))
# Plot regression and figures of merit
rangey = max(y) - min(y)
rangex = max(y_c) - min(y_c)

# Fit a line to the CV vs response
z = np.polyfit(y, y_c, 1)
with plt.style.context(('ggplot')):
    sns.set(rc=settings_sns, font_scale=1.5)
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.scatter(y_c, y, c='red', edgecolors='k')
    #Plot the best fit line
    ax.plot(np.polyval(z,y), y, c='blue', linewidth=1)
    #Plot the ideal 1:1 line
    ax.plot(y, y, color='green', linewidth=1)
    plt.title(attribute)
    plt.xlabel('vorhergesagte Werte')
    plt.ylabel('gemessene Werte')
    plt.legend(['Regressionsgerade \nr$^{2}$ (CV): '+str(round(score_cv,2)),'1:1-Linie','Messwerte'])
    plt.savefig(validation_plot, bbox_inches='tight')

    plt.show()

# Calculate the absolute errors
errors = abs(y_cv-y)
# Print out the mean absolute error (mae)

# Print out the mean absolute error (mae)


print('-------------')
print('n of the test data: {}'.format(len(labels)))
print('Mean of the variable: {:.2f}'.format(np.mean(labels)))
print('Standard deviation of the variable: {:.2f}'.format(np.std(labels)))
print('-------------')
print('Mean Absolute Error: {:.2f}'.format(np.mean(errors)))

print('-------------', file=open(results_txt, "a"))
print('n of the test data: {}'.format(len(labels)), file=open(results_txt, "a"))
print('Mean of the variable: {:.2f}'.format(np.mean(labels)), file=open(results_txt, "a"))
print('Standard deviation of the variable: {:.2f}'.format(np.std(labels)), file=open(results_txt, "a"))
print('-------------', file=open(results_txt, "a"))
print('Mean Absolute Error: {:.2f}'.format(np.mean(errors)), file=open(results_txt, "a"))

mse = mean_squared_error(y, y_cv)
print('Mean squared error: {:.2f}'.format(mse))
print('RMSE: {:.2f}'.format(np.sqrt(mse)))
print('RPD: {:.2f} | How often does RMSE of Prediction fit in the Standard Deviation of the samples'.format(np.std(labels)/np.sqrt(mse)))

print('Mean squared error: {:.2f}'.format(mse), file=open(results_txt, "a"))
print('RMSE: {:.2f}'.format(np.sqrt(mse)), file=open(results_txt, "a"))
print('RPD: {:.2f}| How often does RMSE of Prediction fit in the Standard Deviation of the samples'.format(np.std(labels)/np.sqrt(mse)), file=open(results_txt, "a"))


'''
To put our predictions in perspective, we can calculate an accuracy using
the mean average percentage error subtracted from 100 %.
'''
# Calculate mean absolute percentage error (MAPE)
mape = 100 * (errors / labels)
# Calculate and display accuracy
accuracy = 100 - np.mean(mape)
print('mean absolute percentage error (MAPE): {:.2f} %'.format(np.mean(mape)))
print('accuracy (100 % - mape): {:.2f} %'.format(accuracy))


print('-------------')
# The coefficient of determination: 1 is perfect prediction
print('Coefficient of determination r²: {:.2f}'.format(r2_score(y, y_cv)))


print('mean absolute percentage error (MAPE): {:.2f} %'.format(np.mean(mape)), file=open(results_txt, "a"))
print('accuracy (100 % - mape): {:.2f} %'.format(accuracy), file=open(results_txt, "a"))


print('-------------', file=open(results_txt, "a"))
# The coefficient of determination: 1 is perfect prediction
print('Coefficient of determination r²: {:.2f}'.format(r2_score(y, y_cv)), file=open(results_txt, "a"))

# Predicting the rest of the image

# Take our full image and reshape into long 2d array (nrow * ncol, nband) for classification
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :].reshape(new_shape)

print('Reshaped from {o} to {n}'.format(o=img.shape, n=img_as_array.shape))

img_as_array = np.nan_to_num(img_as_array)
prediction_ = pls_opt.predict(img_as_array)
prediction = prediction_.reshape(img[:, :, 0].shape)
print('Reshaped back to {}'.format(prediction.shape))
# generate mask image from red band
mask = np.copy(img[:,:,0])
mask[mask > 0.0] = 1.0 # all actual pixels have a value of 1.0

# plot mask

plt.imshow(mask)
# mask classification an plot

prediction_ = prediction*mask

plt.subplot(121)
plt.imshow(prediction, cmap=plt.cm.Spectral)
plt.title('classification unmasked')

plt.subplot(122)
plt.imshow(prediction_, cmap=plt.cm.Spectral)
plt.title('classification masked')

plt.show()
# Save regression

cols = img.shape[1]
rows = img.shape[0]

driver = gdal.GetDriverByName("gtiff")
outdata = driver.Create(prediction_map, cols, rows, 1, gdal.GDT_Float32)
outdata.SetGeoTransform(img_ds.GetGeoTransform())##sets same geotransform as input
outdata.SetProjection(img_ds.GetProjection())##sets same projection as input
outdata.GetRasterBand(1).WriteArray(prediction_)
outdata.FlushCache() ##saves to disk!!
print('Image saved to: {}'.format(prediction_map))
# Display images


plt.subplot(121)
plt.imshow(img[:, :, 0], cmap=plt.cm.Greys_r)
plt.title('RS image - first band')

plt.subplot(122)
plt.imshow(prediction_, cmap=plt.cm.Spectral, vmin= y.min(), vmax=y.max())
plt.title('Prediction of ' + attribute)
plt.colorbar()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17565663036515056, 'fraction_numerical': 0.016976297245355543, 'mean_word_length': 4.233322158900436, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 12, 'https://': 2, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '180655', 'license': 'MIT', 'n_tokens_mistral': 5838, 'n_tokens_neox': 5529, 'n_words': 1274, 'path': 'PLSR_Regression_LOOCV.ipynb', 'repo_name': 'florianbeyer/PhD_Florian_Beyer', 'size': '15617'}