In this chapter, we will employ machine learning techniques to classify a scene using satellite imagery. Specifically, we will utilize scikit-learn to implement two distinct classifiers and subsequently compare their results. To begin, we need to import the following modules.
Code
from datetime import datetime, timedeltaimport xarray as xrimport pystac_clientimport stackstacimport odc.stacimport rioxarrayimport geopandas as gpdfrom odc.geo.geobox import GeoBoxfrom dask.diagnostics import ProgressBarfrom rasterio.enums import Resamplingfrom shapely.geometry import Polygon, mappingimport cmcrameri as cmcimport numpy as npimport matplotlib.pyplot as pltimport seaborn as snsimport pandas as pdfrom pathlib import Path# Scikit Learnfrom sklearn.naive_bayes import GaussianNBfrom sklearn.metrics import confusion_matrixfrom sklearn.metrics import classification_reportfrom sklearn.ensemble import RandomForestClassifierfrom sklearn.model_selection import train_test_splitimport matplotlib.colors as colors
Before we start, we need to load the data. We will use odc-stac to obtain data from Earth Search by Element 84. Here we define the area of interest and the time frame, aswell as the EPSG code and the resolution.
1.1.1 Searching in the Catalog
The module odc-stac provides access to free, open source satelite data. To retrieve the data, we must define several parameters that specify the location and time period for the satellite data. Additionally, we must specify the data collection we wish to access, as multiple collections are available. In this example, we will use multispectral imagery from the Sentinel-2 satellite.
We will now focus on the area south-east of Vienna, where the Nationalpark Donauauen is situated. The time frame we are interested in is the beginning of May 2024. After passing these parameters to the stac-catalog we have found 10 scenes that we can use for our analysis.
1.1.2 Loading the Data
Now we will load the data directly into an xarray dataset, which we can use to perform computations on the data. xarray is a powerful library for working with multi-dimensional arrays, making it well-suited for handling satellite data.
Here’s how we can load the data using odc-stac and xarray:
With the image data now in our possession, we can proceed with computations and visualizations.
First, we define a mask to exclude cloud cover and areas with missing data. Subsequently, we create a composite median image, where each pixel value represents the median value across all the scenes we have identified. This approach helps to eliminate clouds and outliers present in some of the images, thereby providing a clearer and more representative visualization of the scene.
Code
# define a mask for valid pixels (non-cloud)def is_valid_pixel(data):# include only vegetated, not_vegitated, water, and snowreturn ((data >3) & (data <7)) | (data==11)ds_odc['valid'] = is_valid_pixel(ds_odc.scl)#ds_odc.valid.sum("time").plot()def avg(ds):return (ds / ds.max() *2)# compute the masked medianrgb_median = ( ds_odc[['red', 'green', 'blue']] .where(ds_odc.valid) .to_dataarray(dim="band") .transpose(..., "band") .median(dim="time") .astype(int))rgb_comp = avg(rgb_median)plot = rgb_comp.plot.imshow(rgb="band", figsize=(8, 8))plot.axes.set_title(f"RGB - Median Composite\n{start_date.strftime('%d.%m.%Y')} - {end_date.strftime('%d.%m.%Y')}")plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-2593019971002186.0..2.0].
1.2.2 False Color Image
In addition to the regular RGB Image, we can swap any of the bands from the visible spectrum with any other bands. In this specific case the red band has been changed to the near infrared band. This allows us to see vegetated areas more clearly, since they now appear in a bright red color. This is due to the fact that plants absorb regular red light while reflecting near infrared light (NASA 2020).
Code
# compute the false color imagefc_median = ( ds_odc[['nir', 'green', 'blue']] .where(ds_odc.valid) .to_dataarray(dim="band") .transpose(..., "band") .median(dim="time") .astype(int))fc_comp = avg(fc_median)plot = fc_comp.plot.imshow(rgb="band", figsize=(8, 8))plot.axes.set_title(f"False Color - Median Composite\n{start_date.strftime('%d.%m.%Y')} - {end_date.strftime('%d.%m.%Y')}")plt.show()
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-2368915381239187.5..2.0].
1.2.3 NDVI Image
To get an first impression of the data, we can calculate the NDVI (Normalized Difference Vegetation Index) and plot it. The NDVI is calculated by useing the following formula. (Rouse et al. 1974)
\[
NDVI = \frac{NIR - Red}{NIR + Red}
\]
This gives us a good overview of the vegetation in the area. The values can range from -1 to 1 where the following meanings are associated with these values:
-1 to 0 indicate dead plants or inanimate objects
0 to 0.33 are unhealthy plants
0.33 to 0.66 are moderatly healthy plants
0.66 to 1 are very healthy plants
Code
# Normalized Difference Vegetation Index (NDVI)def normalized_difference(a, b):return (a - b*1.) / (a + b)ndvi = normalized_difference(ds_odc.nir, ds_odc.red)ndvi.median(dim="time").plot.imshow(cmap='cmc.cork', vmin=-1, vmax=1).axes.set_title('NDVI')plt.show()
1.3 Classification
In this chapter, we will classify the satellite data to identify forested areas within the scene. By using supervised machine learning techniques, we can train classifiers to distinguish between forested and non-forested regions based on the training data we provide. We will explore two different classifiers and compare their performance in accurately identifying forest areas.
1.3.1 Regions of Interest
Since this is a supervised classification, we need to have some training data. Therefore we need to define areas or regions, which we are certain represent the feature which we are classifiying. In this case we are interested in forested areas and regions that are definitly not forested. These regions will be used to train our classifiers.
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Got range [-2593019971002186.0..2.0].
1.3.2 Data Preparation
In addition to the Regions of Interest we will extract the specific bands from the loaded dataset that we intend to use for the classification, which are the red, green, blue and near-infrared bands, although other bands can also be utilized. Using these bands, we will create both a training and a testing dataset. The training dataset will be used to train the classifier, while the testing dataset will be employed to evaluate its performance.
Code
# Classifiying dataset (only necessary bands)bands = ['red', 'green', 'blue', 'nir']ds_class = ( ds_odc[bands] .where(ds_odc.valid) .median(dim="time"))ds_class = avg(ds_class)ds_class = ds_class.fillna(0)def clip_array(ds:xr.Dataset, polygons): clipped = ds.rio.clip(polygons, invert=False, all_touched=False, drop=True) clipped_nan = clipped.where(clipped == ds)return clipped_nan# Dictionaries with Dataarrays, each clipped by a Polygondata_dict_feat = {idx: clip_array(ds_class, polygon) for idx, polygon in forest_areas.items()}data_dict_nonfeat = {idx: clip_array(ds_class, polygon) for idx, polygon in nonforest_areas.items()}# Reshape the polygon dataarrays to get a tuple (one value per band) of pixel valuesfeat_data = [xarray.to_array().values.reshape(len(bands),-1).T for xarray in data_dict_feat.values()] # replaced median_data_dict_feat with data_dict_featnonfeat_data = [xarray.to_array().values.reshape(len(bands),-1).T for xarray in data_dict_nonfeat.values()] # replaced median_data_dict_feat with data_dict_feat# The rows of the different polygons are concatenated to a single array for further processingfeat_values = np.concatenate(feat_data)nonfeat_values = np.concatenate(nonfeat_data)# Drop Nan ValuesX_feat_data = feat_values[~np.isnan(feat_values).any(axis=1)]X_nonfeat_data = nonfeat_values[~np.isnan(nonfeat_values).any(axis=1)]# Creating Output Vector (1 for pixel is features; 0 for pixel is not feature)y_feat_data = np.ones(X_feat_data.shape[0])y_nonfeat_data = np.zeros(X_nonfeat_data.shape[0])# Concatnate all Classes for training X = np.concatenate([X_feat_data, X_nonfeat_data])y = np.concatenate([y_feat_data, y_nonfeat_data])# Split into Training and Testing Data.X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=42)
Now that we have prepared the training and testing data, we will create an image array of the actual scene that we intend to classify. This array will serve as the input for our classification algorithms, allowing us to apply the trained classifiers to the entire scene and identify the forested and non-forested areas accurately.
Now that we have prepared all the needed data, we can begin the actual classification process.
We will start with a Naive Bayes classifier. First, we will train the classifier using our training dataset. Once trained, we will apply the classifier to the actual image to identify the forested and non-forested areas.
Code
# Naive Bayes initialization and trainingnb = GaussianNB()nb_test = nb.fit(X_train, y_train)nb_predict = nb.predict(X_test)# Prediction on imagenb_predict_img = nb.predict(X_image_data)nb_predict_img = nb_predict_img.reshape(ds_class.sizes['latitude'], ds_class.sizes['longitude'])# Adding the Naive Bayes Prediction to the datasetds_class['NB-forest'] = xr.DataArray(nb_predict_img, dims=['latitude', 'longitude'], coords={'longitude': ds_class['longitude'], 'latitude': ds_class['latitude']})
To evaluate the effectiveness of the classification, we will plot the image predicted by the classifier. Additionally, we will examine the Classification Report and the Confusion Matrix to gain further insights into the classifier’s performance.
To ensure our results are robust, we will explore an additional classifier. In this section, we will use the Random Forest classifier. The procedure for using this classifier is the same as before: we will train the classifier using our training dataset and then apply it to the actual image to classify the scene.
Code
# Random Forest initialization and trainingrf = RandomForestClassifier(n_estimators=100)rf_test = rf.fit(X_train, y_train)rf_predict = rf.predict(X_test)# Prediction on imagerf_predict_img = rf.predict(X_image_data)rf_predict_img = rf_predict_img.reshape(ds_class.sizes['latitude'], ds_class.sizes['longitude'])# Adding the Random Forest Prediction to the datasetds_class['RF-forest'] = xr.DataArray(rf_predict_img, dims=['latitude', 'longitude'], coords={'longitude': ds_class['longitude'], 'latitude': ds_class['latitude']})plot = ds_class['RF-forest'].plot.imshow(cmap=cmap_green, cbar_kwargs={'ticks': [0.25,0.75]})cbar = plot.colorbarcbar.set_ticklabels(['non-forest', 'forest'])plot.axes.set_title('Random Forest Classification')plt.show()# Print the Classification reportprint("RANDOM FOREST: \n "+ classification_report(y_test, rf_predict))# Print the confusion matrixcon_mat_rf = pd.DataFrame(confusion_matrix(y_test, rf_predict), index=['Actual Negative', 'Actual Positive'], columns=['Predicted Negative', 'Predicted Positive'])display(con_mat_rf)
We can already see from the classification reports and the confusion matrices that the Random Forest classifier has outperformed the Naive Bayes classifier. This is particularly evident from the lower values in the secondary diagonal, indicating minimal False Positives and False Negatives. It appears that the Naive Bayes classifier is more sensitive to False Positives, resulting in a higher rate of incorrect classifications.
1.3.5 Comparison of the Classificators
To gain a more in-depth understanding of the classifiers’ performance, we will compare their results. Specifically, we will identify the areas where both classifiers agree and the areas where they disagree. This comparison will provide valuable insights into the strengths and weaknesses of each classifier, allowing us to better assess their effectiveness in identifying forested and non-forested regions.
The areas where both classifiers agree include the larger forested regions, such as the Nationalpark Donau-Auen and the Leithagebirge. Additionally, both classifiers accurately identified the urban areas of Vienna and correctly excluded them from being classified as forested.
Code
# Plot only one class, either None (0), Naive Bayes (1), Random Forest (2), or Both (3)fig, axs = plt.subplots(2,2, figsize=(8,8))ax = axs.ravel()for i inrange(4): ax[i].imshow(double_clf==i, cmap='cmc.oleron_r', interpolation='none') category = ['by None', 'only by Naive Bayes', 'only by Random Forest', 'by Both'][i] title ='Areas classified '+ category ax[i].set_title(title) ax[i].set_axis_off()plt.tight_layout()
When plotting the classified areas individually, we observe that the Random Forest classifier mistakenly identified the Danube River as a forested area. Conversely, the Naive Bayes classifier erroneously classified a significant amount of cropland as forest.
Finally, by analyzing the proportion of forested areas within the scene, we find that approximately 18% of the area is classified as forest, while around 66% is classified as non-forest. The remaining areas, which include water bodies and cropland, fall into less clearly defined categories.
The accompanying bar chart illustrates the distribution of these classifications, highlighting the percentage of forested areas, non-forested areas, and regions classified by only one of the two classifiers. This visual representation helps to quantify the areas of agreement and disagreement between the classifiers, providing a clearer picture of their performance.
Code
counts = {}for num inrange(0,4): num_2_class = {0: 'None', 1: 'Naive Bayes', 2: 'Random Forest', 3: 'Both'} counts[num_2_class[num]] =int((double_clf == num).sum().values)class_counts_df = pd.DataFrame(list(counts.items()), columns=['Class', 'Count'])class_counts_df['Percentage'] = (class_counts_df['Count'] / class_counts_df['Count'].sum())*100ax = class_counts_df.plot.bar(x='Class', y='Percentage', rot=0, color='darkgreen', ylim=(0,100), title='Classified Areas per Classificator (%)')# Annotate the bars with the percentage valuesfor p in ax.patches: ax.annotate(f'{p.get_height():.1f}%', (p.get_x() + p.get_width() /2., p.get_height()), ha='center', va='center', xytext=(0, 9), textcoords='offset points')
1.4 Conclusion
In this chapter, we utilized machine learning to classify satellite imagery into forested and non-forested areas, comparing Naive Bayes and Random Forest classifiers. The Random Forest classifier generally outperformed Naive Bayes, with fewer errors in classification, although it misclassified the Danube River as forested, while Naive Bayes incorrectly identified cropland as forest. The analysis, supported by the bar chart, revealed that about 18% of the scene was classified as forest, 66% as non-forest, and the remainder included ambiguous categories. This comparison highlights the strengths and limitations of each classifier, underscoring the need for careful selection and evaluation of classification methods.
Rouse, John Wilson, Rüdiger H Haas, John A Schell, Donald W Deering, et al. 1974. “Monitoring Vegetation Systems in the Great Plains with ERTS.”NASA Spec. Publ 351 (1): 309.
