Uncomment the following line to install geemap if needed.
# !pip install geemap
Machine Learning with Earth Engine - Accuracy Assessment¶
Supervised classification algorithms available in Earth Engine¶
Source: https://developers.google.com/earth-engine/classification
The Classifier package handles supervised classification by traditional ML algorithms running in Earth Engine. These classifiers include CART, RandomForest, NaiveBayes and SVM. The general workflow for classification is:
- Collect training data. Assemble features which have a property that stores the known class label and properties storing numeric values for the predictors.
- Instantiate a classifier. Set its parameters if necessary.
- Train the classifier using the training data.
- Classify an image or feature collection.
- Estimate classification error with independent validation data.
The training data is a FeatureCollection with a property storing the class label and properties storing predictor variables. Class labels should be consecutive, integers starting from 0. If necessary, use remap() to convert class values to consecutive integers. The predictors should be numeric.
To assess the accuracy of a classifier, use a ConfusionMatrix. The sample() method generates two random samples from the input data: one for training and one for validation. The training sample is used to train the classifier. You can get resubstitution accuracy on the training data from classifier.confusionMatrix(). To get validation accuracy, classify the validation data. This adds a classification property to the validation FeatureCollection. Call errorMatrix() on the classified FeatureCollection to get a confusion matrix representing validation (expected) accuracy.
import ee
import geemap
Create an interactive map¶
Map = geemap.Map()
Map
Add data to the map¶
Let's add the USGS National Land Cover Database, which can be used to create training data with class labels.
NLCD2016 = ee.Image('USGS/NLCD/NLCD2016').select('landcover')
Map.addLayer(NLCD2016, {}, 'NLCD 2016')
Load the NLCD metadata to find out the Landsat image IDs used to generate the land cover data.
NLCD_metadata = ee.FeatureCollection("users/giswqs/landcover/NLCD2016_metadata")
Map.addLayer(NLCD_metadata, {}, 'NLCD Metadata')
# point = ee.Geometry.Point([-122.4439, 37.7538]) # Sanfrancisco, CA
# point = ee.Geometry.Point([-83.9293, 36.0526]) # Knoxville, TN
point = ee.Geometry.Point([-88.3070, 41.7471]) # Chicago, IL
metadata = NLCD_metadata.filterBounds(point).first()
region = metadata.geometry()
metadata.get('2016on_bas').getInfo()
doy = metadata.get('2016on_bas').getInfo().replace('LC08_', '')
doy
ee.Date.parse('YYYYDDD', doy).format('YYYY-MM-dd').getInfo()
start_date = ee.Date.parse('YYYYDDD', doy)
end_date = start_date.advance(1, 'day')
image = (
ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterBounds(point)
.filterDate(start_date, end_date)
.first()
.select('B[1-7]')
.clip(region)
)
vis_params = {'min': 0, 'max': 3000, 'bands': ['B5', 'B4', 'B3']}
Map.centerObject(point, 8)
Map.addLayer(image, vis_params, "Landsat-8")
Map
nlcd_raw = NLCD2016.clip(region)
Map.addLayer(nlcd_raw, {}, 'NLCD')
Prepare for consecutive class labels¶
In this example, we are going to use the USGS National Land Cover Database (NLCD) to create label dataset for training.
First, we need to use the remap() function to turn class labels into consecutive integers.
raw_class_values = nlcd_raw.get('landcover_class_values').getInfo()
print(raw_class_values)
n_classes = len(raw_class_values)
new_class_values = list(range(0, n_classes))
new_class_values
class_palette = nlcd_raw.get('landcover_class_palette').getInfo()
print(class_palette)
nlcd = nlcd_raw.remap(raw_class_values, new_class_values).select(
['remapped'], ['landcover']
)
nlcd = nlcd.set('landcover_class_values', new_class_values)
nlcd = nlcd.set('landcover_class_palette', class_palette)
Map.addLayer(nlcd, {}, 'NLCD')
Map
Make training data¶
# Make the training dataset.
points = nlcd.sample(
**{
'region': region,
'scale': 30,
'numPixels': 5000,
'seed': 0,
'geometries': True, # Set this to False to ignore geometries
}
)
Map.addLayer(points, {}, 'training', False)
print(points.size().getInfo())
print(points.first().getInfo())
Split training and testing¶
# Use these bands for prediction.
bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7']
# This property of the table stores the land cover labels.
label = 'landcover'
# Overlay the points on the imagery to get training.
sample = image.select(bands).sampleRegions(
**{'collection': points, 'properties': [label], 'scale': 30}
)
# Adds a column of deterministic pseudorandom numbers.
sample = sample.randomColumn()
split = 0.7
training = sample.filter(ee.Filter.lt('random', split))
validation = sample.filter(ee.Filter.gte('random', split))
training.first().getInfo()
validation.first().getInfo()
Train the classifier¶
In this examples, we will use random forest classification.
classifier = ee.Classifier.smileRandomForest(10).train(training, label, bands)
Classify the image¶
# Classify the image with the same bands used for training.
result = image.select(bands).classify(classifier)
# # Display the clusters with random colors.
Map.addLayer(result.randomVisualizer(), {}, 'classfied')
Map
Render categorical map¶
To render a categorical map, we can set two image properties: classification_class_values and classification_class_palette. We can use the same style as the NLCD so that it is easy to compare the two maps.
class_values = nlcd.get('landcover_class_values').getInfo()
print(class_values)
class_palette = nlcd.get('landcover_class_palette').getInfo()
print(class_palette)
landcover = result.set('classification_class_values', class_values)
landcover = landcover.set('classification_class_palette', class_palette)
Map.addLayer(landcover, {}, 'Land cover')
Map
Visualize the result¶
print('Change layer opacity:')
cluster_layer = Map.layers[-1]
cluster_layer.interact(opacity=(0, 1, 0.1))
Add a legend to the map¶
Map.add_legend(builtin_legend='NLCD')
Map
Accuracy assessment¶
Training dataset¶
confusionMatrix() computes a 2D confusion matrix for a classifier based on its training data (ie: resubstitution error). Axis 0 of the matrix correspond to the input classes (i.e., reference data), and axis 1 to the output classes (i.e., classification data). The rows and columns start at class 0 and increase sequentially up to the maximum class value
train_accuracy = classifier.confusionMatrix()
train_accuracy.getInfo()
Overall Accuracy essentially tells us out of all of the reference sites what proportion were mapped correctly. The overall accuracy is usually expressed as a percent, with 100% accuracy being a perfect classification where all reference site were classified correctly. Overall accuracy is the easiest to calculate and understand but ultimately only provides the map user and producer with basic accuracy information.
train_accuracy.accuracy().getInfo()
The Kappa Coefficient is generated from a statistical test to evaluate the accuracy of a classification. Kappa essentially evaluates how well the classification performed as compared to just randomly assigning values, i.e. did the classification do better than random. The Kappa Coefficient can range from -1 t0 1. A value of 0 indicated that the classification is no better than a random classification. A negative number indicates the classification is significantly worse than random. A value close to 1 indicates that the classification is significantly better than random.
train_accuracy.kappa().getInfo()
Producer's Accuracy is the map accuracy from the point of view of the map maker (the producer). This is how often are real features on the ground correctly shown on the classified map or the probability that a certain land cover of an area on the ground is classified as such. The Producer's Accuracy is complement of the Omission Error, Producer's Accuracy = 100%-Omission Error. It is also the number of reference sites classified accurately divided by the total number of reference sites for that class.
train_accuracy.producersAccuracy().getInfo()
The Consumer's Accuracy is the accuracy from the point of view of a map user, not the map maker. the User's accuracy essentially tells use how often the class on the map will actually be present on the ground. This is referred to as reliability. The User's Accuracy is complement of the Commission Error, User's Accuracy = 100%-Commission Error. The User's Accuracy is calculating by taking the total number of correct classifications for a particular class and dividing it by the row total.
train_accuracy.consumersAccuracy().getInfo()
Validation dataset¶
validated = validation.classify(classifier)
validated.first().getInfo()
errorMatrix computes a 2D error matrix for a collection by comparing two columns of a collection: one containing the actual values, and one containing predicted values.
test_accuracy = validated.errorMatrix('landcover', 'classification')
test_accuracy.getInfo()
test_accuracy.accuracy().getInfo()
test_accuracy.kappa().getInfo()
test_accuracy.producersAccuracy().getInfo()
test_accuracy.consumersAccuracy().getInfo()
Download confusion matrix¶
import csv
import os
out_dir = os.path.join(os.path.expanduser('~'), 'Downloads')
training_csv = os.path.join(out_dir, 'train_accuracy.csv')
testing_csv = os.path.join(out_dir, 'test_accuracy.csv')
with open(training_csv, "w", newline="") as f:
writer = csv.writer(f)
writer.writerows(train_accuracy.getInfo())
with open(testing_csv, "w", newline="") as f:
writer = csv.writer(f)
writer.writerows(test_accuracy.getInfo())
Reclassify land cover map¶
landcover = landcover.remap(new_class_values, raw_class_values).select(
['remapped'], ['classification']
)
landcover = landcover.set('classification_class_values', raw_class_values)
landcover = landcover.set('classification_class_palette', class_palette)
Map.addLayer(landcover, {}, 'Final land cover')
Map
Export the result¶
Export the result directly to your computer:
import os
out_dir = os.path.join(os.path.expanduser('~'), 'Downloads')
out_file = os.path.join(out_dir, 'landcover.tif')
geemap.ee_export_image(landcover, filename=out_file, scale=900)
Export the result to Google Drive:
geemap.ee_export_image_to_drive(
landcover, description='landcover', folder='export', scale=900
)
Created: 2020-08-12