Demonstrating Land Cover Land Use (LCLU) Classification Methods Using Google Earth Engine and Random Forest and SVM

---

Background

---

The following project will walk through some basic Land Cover Land Use (LCLU) classification methods using Google Earth engine, using Bhopal, India as an example. The city features prominent lakes and has some neighboring croplands.

I will use two supervised classification methods, including Random Forest and Support Vector Machines (SVMs). For the Random Forest classification, I will use LCLU data available through ESRI to train our model. Our SVM model will use some crude demarcations of buildings, water, crops and shrubs as inputs to train the model.

Supervised Classification Using Random Forest

The random forest classifier for this project will be trained on LCLU data available through ESRI. The ESRI land cover data is itself the output of a predictive model, with years of global land cover maps "generated with Impact Observatory’s deep learning AI land classification model, trained on billions of human-labeled image pixels from the National Geographic Society". The ESRI data has an assessed average accuracy of > 75%. For more information on the data set, see here and here.

# Import earth engine libraries, authenticate session and initialize:
import ee

ee.Authenticate() 
ee.Initialize()

To authorize access needed by Earth Engine, open the following URL in a web browser and follow the instructions:

https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=GNXhLPy7GcOemkZTC1bWD3krfuZfDdG-oKR2eX8Zvo4&tc=bcjcp7hV0g_T74cBEwBhH7YirO9suNykJGgo0a37la0&cc=N_x7UwJ2y002OqzGih3kPtofNjDG7M1JHAqYGHrCzEg

The authorization workflow will generate a code, which you should paste in the box below.

Enter verification code: 4/1AbUR2VMHirQWuQKHt745swRkA93pCoPvGvGrFnUKu4qIZdfbBiOT1vhXICY

Successfully saved authorization token.

# Let's define a point of interest southeast of the city center of Bhopal, where the town borders some farmland
# We will also define a region of interest (ROI) as being a buffer of 10km around this point

poi = ee.Geometry.Point([77.511808, 23.170332])
roi = poi.buffer(20000)

# Let's load an image of Bhopal by filtering through the Landsat 8 satellite collection

# We apply a date filter of 1 year to grab all images available in that period
# We further apply a location filter by providing our ROI coordinates for Bhopal to find the images that include it
# Then we sort the images based on the CLOUD_COVER tag, selecting the first image, i.e. the one with the least clouds
# We clip the final image to be the region we defined above

image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
.filterDate('2021-01-01', '2021-12-31') \
.filterBounds(roi) \
.sort('CLOUD_COVER') \
.first() \
.clip(roi)

# Rather than manually creating training data, we'll be using ESRI's 2020 land cover data, also available in EE
# Note that the name of the image collection notes the 10m resolution of the dataset
# We clip the ESRI data to match our image's geometry using .clip() and .geometry()
# .mosaic() composites all of the images in the collection
lc = ee.ImageCollection("projects/sat-io/open-datasets/landcover/ESRI_Global-LULC_10m").mosaic().clip(image.geometry())

# We see that the ESRI LULC data has a single 'band' of data, containing the class labels

# This band is called 'b1' which we will reference later

lc.getInfo()

{'type': 'Image',
 'bands': [{'id': 'b1',
   'data_type': {'type': 'PixelType',
    'precision': 'int',
    'min': 0,
    'max': 255},
   'dimensions': [3, 2],
   'origin': [76, 22],
   'crs': 'EPSG:4326',
   'crs_transform': [1, 0, 0, 0, 1, 0]}],
 'properties': {'system:footprint': {'type': 'Polygon',
   'coordinates': [[[77.511808, 23.350309319029638],
     [77.45705650110803, 23.343136748918305],
     [77.4066797562988, 23.32219219265991],
     [77.3646990829293, 23.289148788785194],
     [77.33445871483706, 23.24664460109604],
     [77.31835882553887, 23.198070174203366],
     [77.31766669640452, 23.14729639240468],
     [77.3324203332756, 23.098365009910154],
     [77.36143077530768, 23.055166897666048],
     [77.40238124673525, 23.021133613452623],
     [77.45201387653591, 22.998966436164196],
     [77.50638843789102, 22.990423752225055],
     [77.56119271486486, 22.996182977979235],
     [77.6120808058517, 23.015787432868606],
     [77.65501394963238, 23.04768215079655],
     [77.68657830633921, 23.089335943084212],
     [77.7042555412508, 23.137440518237405],
     [77.7066250493335, 23.188171533442336],
     [77.69348120547018, 23.237491508625666],
     [77.66585503442744, 23.281470968115077],
     [77.62593692575653, 23.316602336427817],
     [77.57690503001459, 23.340081266530493],
     [77.52267210286468, 23.350032343726497],
     [77.511808, 23.350309319029638]]]}}}

import numpy as np

# Get the 'b1' band as an Earth Engine Image
b1_image = lc.select('b1')

# Convert the 'b1' band to a NumPy array
scale = 30  # Adjust the scale value as needed
b1_array = np.array(b1_image.reduceRegion(
    reducer=ee.Reducer.toList(),
    geometry=lc.geometry(),
    scale=scale,
    maxPixels=1e9
).get('b1').getInfo())

# Print the NumPy array
print(b1_array)

[5 5 5 ... 5 5 5]

# Get the unique values and their counts
unique_values, value_counts = np.unique(b1_array, return_counts=True)

# Print the unique values and their counts
print("Total pixels:",value_counts.sum())
for value, count in zip(unique_values, value_counts):
    print("Value", value, "Count:", count)
    
import matplotlib.pyplot as plt

# Create a bar chart
plt.bar(unique_values, value_counts)

# Add labels and title
plt.xlabel('Unique Values')
plt.ylabel('Count')
plt.title('Value Counts')

# Display the chart
plt.show()

Total pixels: 1503989
Value 1 Count: 45157
Value 2 Count: 62323
Value 3 Count: 2493
Value 4 Count: 3153
Value 5 Count: 760678
Value 6 Count: 353729
Value 7 Count: 276404
Value 8 Count: 52

import numpy as np
import matplotlib.pyplot as plt

# Get the unique values and their counts
unique_values, value_counts = np.unique(b1_array, return_counts=True)

# Print the unique values and their counts
print("Total pixels:", value_counts.sum())
for value, count in zip(unique_values, value_counts):
    print("Value", value, "Count:", count)

Total pixels: 1503989
Value 1 Count: 45157
Value 2 Count: 62323
Value 3 Count: 2493
Value 4 Count: 3153
Value 5 Count: 760678
Value 6 Count: 353729
Value 7 Count: 276404
Value 8 Count: 52

# Define the labels for the values, note that class 9 and 10 ('snow/ice' and 'clouds') are not present above:
lc_labels = ['Water', 'Trees', 'Grass', 'Flooded Vegetation', 'Crops', 'Shrub', 'Built Area', 'Bare Ground']

# Create a bar chart with labeled x-axis
plt.bar(lc_labels, value_counts)

# Add labels and title
plt.xlabel('Land Cover')
plt.ylabel('Count')
plt.title('Value Counts')

# Rotate the x-axis labels for better readability
plt.xticks(rotation=45)

# Display the chart
plt.show()

# The classes are clearly imbalanced, and we should assume the 'bare ground' will not be picked up by our classifier

#######  SETTING UP OUR GROUND TRUTH  #######

# Our ESRI data will serve as our 'ground truth'

# We take our original image, add the 'ground truth' data by adding the ESRI data as a band ('addBands')
# Then we sample it with .sample() to select some number of randomly selected pixels
# For initial explanatory purposes only let's select 10 random pixels

sample = image.addBands(lc).sample(**{
  'region': image.geometry(),
  'numPixels': 10,
  'seed': 0
})

# We now assign an additional datafield to the sample image containing a randomized value using an EE function
sample = sample.randomColumn()

# Based on this randomized value, we will define 80% of our sample image as the training set, and 20% for testing
# We do this by choosing values that were below 0.80, and values above 0.80 in the randomized values column
sample_train = sample.filter('random <= 0.8')
sample_test = sample.filter('random > 0.8')

# We will eventually train our random forest model on the sample_train (i.e. 80% of the sample pixels we picked)
# We can then see how well the model performed by seeing it's performance on the sample_train

# We currently only selected 10 pixels for demonstration purposes. Here is what the sample_train data will look like:

# It includes 9 of the 10 pixels that had a 'random' value less than 0.8, namely pixels: 0, 1, 2, 3, 5, 6, 7, 8, 9
# We also see the band values

sample_train.getInfo()

{'type': 'FeatureCollection',
 'columns': {'random': '<any>'},
 'features': [{'type': 'Feature',
   'geometry': None,
   'id': '0',
   'properties': {'B1': 150,
    'B10': 3005,
    'B11': 2990,
    'B2': 285,
    'B3': 655,
    'B4': 601,
    'B5': 2929,
    'B6': 1468,
    'B7': 827,
    'b1': 5,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.31952938406342224,
    'sr_aerosol': 224}},
  {'type': 'Feature',
   'geometry': None,
   'id': '1',
   'properties': {'B1': 521,
    'B10': 3124,
    'B11': 3096,
    'B2': 603,
    'B3': 772,
    'B4': 987,
    'B5': 1711,
    'B6': 2474,
    'B7': 2140,
    'b1': 6,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.3620352761013097,
    'sr_aerosol': 96}},
  {'type': 'Feature',
   'geometry': None,
   'id': '2',
   'properties': {'B1': 655,
    'B10': 3030,
    'B11': 3013,
    'B2': 717,
    'B3': 924,
    'B4': 1113,
    'B5': 2999,
    'B6': 3547,
    'B7': 2344,
    'b1': 6,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.03871546282965732,
    'sr_aerosol': 66}},
  {'type': 'Feature',
   'geometry': None,
   'id': '3',
   'properties': {'B1': 631,
    'B10': 3123,
    'B11': 3100,
    'B2': 741,
    'B3': 967,
    'B4': 1260,
    'B5': 2112,
    'B6': 3444,
    'B7': 2719,
    'b1': 6,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.4959765478263076,
    'sr_aerosol': 96}},
  {'type': 'Feature',
   'geometry': None,
   'id': '5',
   'properties': {'B1': 85,
    'B10': 3003,
    'B11': 2987,
    'B2': 206,
    'B3': 529,
    'B4': 421,
    'B5': 3324,
    'B6': 1166,
    'B7': 565,
    'b1': 5,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.6547752781704712,
    'sr_aerosol': 224}},
  {'type': 'Feature',
   'geometry': None,
   'id': '6',
   'properties': {'B1': 466,
    'B10': 3093,
    'B11': 3069,
    'B2': 612,
    'B3': 879,
    'B4': 981,
    'B5': 1652,
    'B6': 1649,
    'B7': 1255,
    'b1': 7,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.7481828715981919,
    'sr_aerosol': 160}},
  {'type': 'Feature',
   'geometry': None,
   'id': '7',
   'properties': {'B1': 284,
    'B10': 3038,
    'B11': 3019,
    'B2': 395,
    'B3': 719,
    'B4': 697,
    'B5': 2764,
    'B6': 1558,
    'B7': 951,
    'b1': 5,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.5949123548364291,
    'sr_aerosol': 160}},
  {'type': 'Feature',
   'geometry': None,
   'id': '8',
   'properties': {'B1': 627,
    'B10': 3132,
    'B11': 3102,
    'B2': 710,
    'B3': 913,
    'B4': 1129,
    'B5': 1906,
    'B6': 2781,
    'B7': 2289,
    'b1': 6,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.09817803297879668,
    'sr_aerosol': 96}},
  {'type': 'Feature',
   'geometry': None,
   'id': '9',
   'properties': {'B1': 331,
    'B10': 3041,
    'B11': 3031,
    'B2': 486,
    'B3': 944,
    'B4': 824,
    'B5': 2657,
    'B6': 1793,
    'B7': 1111,
    'b1': 5,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.5645825697675269,
    'sr_aerosol': 160}}]}

# Our sample_test includes the one pixel that had a random value greater than 0.80 (#4)

sample_test.getInfo()

{'type': 'FeatureCollection',
 'columns': {'random': '<any>'},
 'features': [{'type': 'Feature',
   'geometry': None,
   'id': '4',
   'properties': {'B1': 140,
    'B10': 3011,
    'B11': 2996,
    'B2': 255,
    'B3': 562,
    'B4': 479,
    'B5': 2729,
    'B6': 1220,
    'B7': 687,
    'b1': 5,
    'pixel_qa': 322,
    'radsat_qa': 0,
    'random': 0.8979374428053774,
    'sr_aerosol': 224}}]}

# If we look at the 'b1' values above, storing the ESRI ground truth values, we only have classes 5, 6, and 7 found
# in the randomly selected 10 pixels. Obviously we need to choose a realistic sample to capture available classes

## Let's find out how many pixels are in our image to begin with

## We do this with the following code snippet. reduceRegion() is used to calculate the pixels.
# The parameter "reducer" specifies a way to calculate pixel values. We've defined it as '.count' to count the pixels
# We define the geometry to count the pixels to be our roi (which in this case matches the image since we clipped it)
# The scale defines the pixel size in meters and the 'maxPixels' being set very large at 1e13 is to prevent an error
# if there are more pixels than the maxPixels

pixel_count = image.reduceRegion(
    reducer=ee.Reducer.count(),
    geometry=roi,
    scale=image.select(0).projection().nominalScale(),
    maxPixels=1e13
).getInfo()

# Print the pixel count
print("Pixel Count:", pixel_count)

Pixel Count: {'B1': 1377374, 'B10': 1377374, 'B11': 1377374, 'B2': 1377374, 'B3': 1377374, 'B4': 1377374, 'B5': 1377374, 'B6': 1377374, 'B7': 1377374, 'pixel_qa': 1377374, 'radsat_qa': 1377374, 'sr_aerosol': 1377374}

# There are about 1.4 million pixels, and we don't need to build our model on the entirety of this.

# Let's select 100,000 random pixels which should be large enough to include the various classes present
sample = image.addBands(lc).sample(**{
  'region': image.geometry(),
  'numPixels': 100000,
  'seed': 10
})

# We assign our randomized value column
sample = sample.randomColumn()

# We define 80% of our sample image as the training set, and 20% for testing
sample_train = sample.filter('random <= 0.8')
sample_test = sample.filter('random > 0.8')

sample_train.aggregate_stats(label).getInfo()

{'max': 8,
 'mean': 5.35656836461126,
 'min': 1,
 'sample_sd': 1.2972588830668113,
 'sample_var': 1.6828806096957505,
 'sum': 429570,
 'sum_sq': 2435978,
 'total_count': 80195,
 'total_sd': 1.297250794888451,
 'total_var': 1.6828596248387184,
 'valid_count': 80195,
 'weight_sum': 80195,
 'weighted_sum': 429570}

sample_test.aggregate_stats(label).getInfo()

{'max': 8,
 'mean': 5.350998003992016,
 'min': 1,
 'sample_sd': 1.3122484938388472,
 'sample_var': 1.7219961095823229,
 'sum': 107234,
 'sum_sq': 608316,
 'total_count': 20040,
 'total_sd': 1.3122157526995104,
 'total_var': 1.721910181632743,
 'valid_count': 20040,
 'weight_sum': 20040,
 'weighted_sum': 107234}

##### Let's now train our random forest classifier on the train data #####

# We define the bands in the original image we want to use as our 'input' variables to determine the class:
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11']

# Our label, i.e. the 'class' we are determining is stored as "b1" in the sample_train variable we created
label = 'b1'

# We will use 5 trees initially to see results

## Train the model ##
rf_classifier = ee.Classifier.smileRandomForest(5).train(**{
  'features': sample_train,
  'classProperty': label,
  'inputProperties': bands
})

# Get a confusion matrix and overall accuracy for the classifier based on the train data which it was built off of:
train_accuracy = rf_classifier.confusionMatrix()
print('Training error matrix:')
for x in train_accuracy.getInfo():
    print(x)
    
print('\nTraining overall accuracy', train_accuracy.accuracy().getInfo())

Training error matrix:
[0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 2281, 0, 0, 6, 56, 2, 22, 0]
[0, 3, 2771, 0, 0, 33, 411, 17, 0]
[0, 0, 1, 99, 0, 9, 8, 9, 0]
[0, 22, 2, 0, 123, 16, 1, 5, 0]
[0, 22, 47, 5, 2, 39831, 340, 511, 0]
[0, 8, 240, 1, 2, 412, 18132, 99, 0]
[0, 14, 31, 2, 1, 953, 189, 13453, 0]
[0, 0, 0, 0, 0, 0, 0, 1, 2]

Training overall accuracy 0.9563189725045202

import numpy as np
from sklearn.metrics import confusion_matrix, classification_report

# Or we can use scikit-learn to also pull up a classification report for more comprehensive and class specific metrics

# Classify the training sample using the trained Random Forest classifier
sample_training_classification = sample_train.classify(rf_classifier)

# Extract the true labels and predicted labels from the training sample as arrays
true_labels_train = sample_training_classification.aggregate_array(label).getInfo()
predicted_labels_train = sample_training_classification.aggregate_array('classification').getInfo()

# Compute the confusion matrix
confusion_matrix_train = confusion_matrix(true_labels_train, predicted_labels_train)

# Compute the classification report
classification_report_train = classification_report(true_labels_train, predicted_labels_train)

# Print the confusion matrix and classification report for training data
print('\nConfusion Matrix (Training Data):\n')
print(confusion_matrix_train)
print('\nClassification Report (Training Data):\n')
print(classification_report_train)

Confusion Matrix (Training Data):

[[ 2281     0     0     6    56     2    22     0]
 [    3  2771     0     0    33   411    17     0]
 [    0     1    99     0     9     8     9     0]
 [   22     2     0   123    16     1     5     0]
 [   22    47     5     2 39831   340   511     0]
 [    8   240     1     2   412 18132    99     0]
 [   14    31     2     1   953   189 13453     0]
 [    0     0     0     0     0     0     1     2]]

Classification Report (Training Data):

              precision    recall  f1-score   support

           1       0.97      0.96      0.97      2367
           2       0.90      0.86      0.88      3235
           3       0.93      0.79      0.85       126
           4       0.92      0.73      0.81       169
           5       0.96      0.98      0.97     40758
           6       0.95      0.96      0.95     18894
           7       0.95      0.92      0.94     14643
           8       1.00      0.67      0.80         3

    accuracy                           0.96     80195
   macro avg       0.95      0.86      0.90     80195
weighted avg       0.96      0.96      0.96     80195

# The accuracy on the train is 96%, and we see that it sometimes confuses crops with shrubs and built areas

print(lc_labels)

['Water', 'Trees', 'Grass', 'Flooded Vegetation', 'Crops', 'Shrub', 'Built Area', 'Bare Ground']

## How about our test/validation data results?

# Let's use the random forest classifier we built and apply it on our sample_test data now then see the accuracy:
sample_validation_classification = sample_test.classify(rf_classifier)

# Extract the true labels and predicted labels from the validation sample as arrays
validation_true_labels = sample_validation_classification.aggregate_array(label).getInfo()
validation_predicted_labels = sample_validation_classification.aggregate_array('classification').getInfo()

# Compute the confusion matrix using scikit-learn
validation_confusion_matrix = confusion_matrix(validation_true_labels, validation_predicted_labels)

# Compute the classification report
validation_classification_report = classification_report(validation_true_labels, validation_predicted_labels, zero_division=0)

# Print the confusion matrix and classification report
print('\nConfusion Matrix (Test Data):\n')
print(validation_confusion_matrix)
print('\nClassification Report (Test Data):\n')
print(validation_classification_report)

Confusion Matrix (Test Data):

[[ 535    4    0    5   42    5   21    0]
 [   1  346    0    0   59  416   21    0]
 [   0    0    3    0    4    8    7    0]
 [  24    1    0   12    5    1    1    0]
 [  20   33    3    6 9155  330  555    0]
 [   6  272    4    0  383 3961  110    0]
 [  10   38    5    0  865  145 2617    0]
 [   0    0    0    0    0    0    1    0]]

Classification Report (Test Data):

              precision    recall  f1-score   support

           1       0.90      0.87      0.89       612
           2       0.50      0.41      0.45       843
           3       0.20      0.14      0.16        22
           4       0.52      0.27      0.36        44
           5       0.87      0.91      0.89     10102
           6       0.81      0.84      0.83      4736
           7       0.79      0.71      0.75      3680
           8       0.00      0.00      0.00         1

    accuracy                           0.83     20040
   macro avg       0.57      0.52      0.54     20040
weighted avg       0.83      0.83      0.83     20040

# We have an 83% accuracy for our validation data.

# Let's redo our random forest classifier with 30 trees:

## Train the model ##
rf_classifier = ee.Classifier.smileRandomForest(30).train(**{
  'features': sample_train,
  'classProperty': label,
  'inputProperties': bands
})


# Classify the training sample using the trained Random Forest classifier
sample_training_classification = sample_train.classify(rf_classifier)

# Extract the true labels and predicted labels from the training sample as arrays
true_labels_train = sample_training_classification.aggregate_array(label).getInfo()
predicted_labels_train = sample_training_classification.aggregate_array('classification').getInfo()

# Compute the confusion matrix
confusion_matrix_train = confusion_matrix(true_labels_train, predicted_labels_train)

# Compute the classification report
classification_report_train = classification_report(true_labels_train, predicted_labels_train)

# Print the confusion matrix and classification report for training data
print('\nConfusion Matrix (Training Data):\n')
print(confusion_matrix_train)
print('\nClassification Report (Training Data):\n')
print(classification_report_train)


## Apply the random forest classifier we built on our sample_test data now then see the accuracy:
sample_validation_classification = sample_test.classify(rf_classifier)

# Extract the true labels and predicted labels from the validation sample as arrays
validation_true_labels = sample_validation_classification.aggregate_array(label).getInfo()
validation_predicted_labels = sample_validation_classification.aggregate_array('classification').getInfo()

# Compute the confusion matrix using scikit-learn
validation_confusion_matrix = confusion_matrix(validation_true_labels, validation_predicted_labels)

# Compute the classification report
validation_classification_report = classification_report(validation_true_labels, validation_predicted_labels, zero_division=0)

# Print the confusion matrix and classification report
print('\nConfusion Matrix (Test Data):\n')
print(validation_confusion_matrix)
print('\nClassification Report (Test Data):\n')
print(validation_classification_report)

Confusion Matrix (Training Data):

[[ 2334     0     0     0    23     1     9     0]
 [    1  3105     0     0    12   113     4     0]
 [    0     0   115     0     6     3     2     0]
 [    7     0     0   154     8     0     0     0]
 [    2     1     0     0 40579   113    63     0]
 [    1    23     0     0   116 18719    35     0]
 [    4     1     0     0   332    43 14263     0]
 [    0     0     0     0     0     0     0     3]]

Classification Report (Training Data):

              precision    recall  f1-score   support

           1       0.99      0.99      0.99      2367
           2       0.99      0.96      0.98      3235
           3       1.00      0.91      0.95       126
           4       1.00      0.91      0.95       169
           5       0.99      1.00      0.99     40758
           6       0.99      0.99      0.99     18894
           7       0.99      0.97      0.98     14643
           8       1.00      1.00      1.00         3

    accuracy                           0.99     80195
   macro avg       0.99      0.97      0.98     80195
weighted avg       0.99      0.99      0.99     80195


Confusion Matrix (Test Data):

[[ 526    0    0    2   53    8   23    0]
 [   3  332    0    0   52  427   29    0]
 [   0    0    0    0    7    8    7    0]
 [  25    1    0   10    6    2    0    0]
 [   7   12    1    3 9243  330  506    0]
 [   3  184    0    1  280 4136  132    0]
 [   6   13    2    1  756  117 2785    0]
 [   0    0    0    0    0    0    1    0]]

Classification Report (Test Data):

              precision    recall  f1-score   support

           1       0.92      0.86      0.89       612
           2       0.61      0.39      0.48       843
           3       0.00      0.00      0.00        22
           4       0.59      0.23      0.33        44
           5       0.89      0.91      0.90     10102
           6       0.82      0.87      0.85      4736
           7       0.80      0.76      0.78      3680
           8       0.00      0.00      0.00         1

    accuracy                           0.85     20040
   macro avg       0.58      0.50      0.53     20040
weighted avg       0.84      0.85      0.85     20040

# We see improvement in the test data but it's marginal.
# Let's go ahead and classify the entire image with our new model

classified_image = image.classify(rf_classifier)

# Import the Folium library.
import folium
from IPython.display import display

# Define a method for displaying Earth Engine image tiles onto the folium map.
def addLayer (self, ee_image_object, vis_params, name):
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
    tiles = map_id_dict['tile_fetcher'].url_format,
    attr = 'Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
    name = name,
    overlay = True,
    control = True
  ).add_to(self)
    
# Apply the defined EE addLayer method to folium
folium.Map.addLayer = addLayer

# Instantiate the map object with Bhopal in the center
Map = folium.Map(location=[23.170332,77.511808], zoom_start = 10)

# Create a google map basemap
basemaps = {'Google Satellite Hybrid': folium.TileLayer(
tiles = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
attr = 'Google',
name = 'Google Satellite Hybrid',
overlay = True,
control = True
)}

# Apply a google satellite image as a basemap to the Map
basemaps['Google Satellite Hybrid'].add_to(Map)

# Display
display(Map)

Make this Notebook Trusted to load map: File -> Trust Notebook

##### Add the original landsat image to our map:

# We need to define the bands for visualizing the satellite image colors: RGB, corresponding with B4, B3, and B2)
visParamsTrue = {'bands': ['B4',  'B3',  'B2'], 'min': 0, 'max': 3000, 'gamma': 1.4}
Map.addLayer(image, visParamsTrue, "Landsat 2021")


##### Add the ESRI LULC image to our map:

# Define a dictionary with the names of the ESRI LULC classes and colors for visualization

labels_and_colors = {
  "names": [
    "Water",
    "Trees",
    "Grass",
    "Flooded Vegetation",
    "Crops",
    "Scrub/Shrub",
    "Built Area",
    "Bare Ground",
    "Snow/Ice",
    "Clouds"
  ],
  
  "colors": [
    "#1A5BAB", # blue
    "#358221", # green
    "#31D431", # bright green
    "#2B9974", # teal
    "#FFDB5C", # yellow
    "#EECCA8", # tan
    "#ED022A", # red
    "#EDE9E4", # light gray/white
    "#FFFFFF", # white
    "#C8C8C8" # dark gray
  ]}

Map.addLayer(lc, {'min':1, 'max':10, 'palette':labels_and_colors['colors']}, 'ESRI LULC')

##### Add our classified image to our map above:
# We use the same color schema just defined above for the ESRI data layer since we want it to match
Map.addLayer(classified_image, {'min':1, 'max':10, 'palette':labels_and_colors['colors']},'Landsat 2021 - Classified')

# Add a layer control panel to the map.
Map.add_child(folium.LayerControl(position='topleft'))

# Display the map:
display(Map)

Make this Notebook Trusted to load map: File -> Trust Notebook

Our RF model, trained on the ESRI data, is able to produce a fairly robust classification that is very close to the ESRI classifications for Bhopal. The LULC classes fairly accurately classify the satellite image, correctly identifying the location of the lakes, river, large croplands, built out city, and forested areas.

Supervised Classification Using SVM

For our SVM model, we will use some custom-demarcated polygons of some classes and see how the model performs after being trained on them. These demarcated regions will be somewhat crude in design, but will serve to illustrate the process of building simple input data for training the model. We will also run our RVM classifier on a larger region than the random forest model.

landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1")

# Pull composite landsat image for 2020 without cloud coverage
# clip to a bhopal region we define

bhopal_region = ee.Geometry.Point([77.41, 23.26]).buffer(50000)

selected_image = ee.Algorithms.Landsat.simpleComposite(**{
     'collection': landsat8.filterDate('2020-01-01', '2020-12-31'),
      'asFloat': True
}).clip(bhopal_region)

# Here are the simple polygons that I've selected to represent examples of the classes:

# A polygon of some buildings in the city center:
bhopal_buildings = ee.Geometry.Polygon([[77.4128293, 23.2626651],[77.4234723, 23.2654249],[77.4355744, 23.2582493],[77.4269055, 23.2516252],[77.4129151, 23.2624286]])

# Using a sample of a lake in Bhopal:
bhopal_water = ee.Geometry.Polygon([[77.3382499, 23.2541704], [77.3491504, 23.2539338], [77.3486355, 23.2478616], [77.3379066, 23.2478616], [77.3381641, 23.2541704]])

# Using a sample of crops by Bhopal:
bhopal_crops = ee.Geometry.Polygon([[77.5181188, 23.1943395], [77.5202109, 23.1952862], [77.5212301, 23.1940634], [77.5208868, 23.1914500], [77.5173249, 23.1916768], [77.5181402, 23.1943099]])

# Using a sample of forest shrubs in Bhopal:
bhopal_shrubs = ee.Geometry.Polygon([[77.6637996, 22.9319460], [77.6778329, 22.9314322], [77.6716531, 22.9236063], [77.6628983, 22.9236853], [77.6638425, 22.9317484]])

# We make a feature collection from the custom geometries we made for buildings, water, crops and shrubs,
# along with identifying the class number for each:
training_polygons = ee.FeatureCollection([
  ee.Feature(bhopal_buildings, {'class': 0}),
  ee.Feature(bhopal_water, {'class': 1}),
  ee.Feature(bhopal_crops, {'class': 2}),
  ee.Feature(bhopal_shrubs, {'class': 3})
])

# Pull the pixels in the sample polygons we made for training the model, and define 'properties' as 'class',
# since that's where we encoded our classes.
# Note that the 'scale' is 30 because that is the resolution for landsat images
training_pixels = selected_image.sampleRegions(**{
    'collection': training_polygons,
    'properties': ['class'],
    'scale': 30
})

# Create an SVM classifier with basic parameters
svm = ee.Classifier.libsvm(**{
    'kernelType': 'RBF',
    'gamma': 0.5,
    'cost': 10
})

# Bands we'll use for our prediction
bands = ['B1', 'B2', 'B3', 'B4', 'B5','B6','B7', 'B10', 'B11']

# Train model
trained_svm = svm.train(training_pixels, 'class', bands)

svm_classified_image = selected_image.classify(trained_svm)

# Apply the pre-defined EE addLayer method to folium
folium.Map.addLayer = addLayer

# Instantiate the folium map object with Bhopal in the center
Map = folium.Map(location=[23.170332,77.511808], zoom_start = 9)

# Create a google map basemap
basemaps = {'Google Satellite Hybrid': folium.TileLayer(
tiles = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
attr = 'Google',
name = 'Google Satellite Hybrid',
overlay = True,
control = True
)}

# Apply a google satellite image as a basemap to the Map
basemaps['Google Satellite Hybrid'].add_to(Map)

# Add the selected satellite image layer to the map object.
Map.addLayer(selected_image, {'bands': ['B4', 'B3', 'B2'], 'max': 0.5, 'gamma': 2}, "Bhopal Region")


# Create colors for classes:
labels_and_colors = {
  "names": [
    "Built Area (0)",
    "Water (1)",
    "Crops (2)",
    "Shrub (3)",
  ],
  
  "colors": [
    "#ED022A",
    "#1A5BAB",
    "#FFDB5C",
    "#44911A",
  ]}

# Set min to 0 and max to 3, to show all 4 classes (0,1,2,3)
Map.addLayer(svm_classified_image, {'min': 0, 'max': 3, 'palette': labels_and_colors['colors']}, 'Bhopal Classified')
    
# Add a layer control panel to the map.
Map.add_child(folium.LayerControl(position = 'topleft'))

# Display
display(Map)

Make this Notebook Trusted to load map: File -> Trust Notebook

Even though our input data was incredibly crude, our SVM classifier was still able to identify bodies of water and some of the building concentrations in the city. Using better input data will obviously result in a much better classification, and image segmentation (e.g., using segment-geospatial: https://samgeo.gishub.org/) can be used to produce such input data.