In [1]:
'''
Name: Josh Virene
Date: 9 March, 2023
Class: NR426: Programming for GIS
Purpose: The purpose of this lab is to accomplish landcover classification for a study area in Big Sky, Montana using
satellite imagery- National Agricultural Imagery Program (NAIP), and the National Landcover Dataset (NLCD)

Program(s): Jupyter Notebook, Google Earth Engine, Python
Data sources: 
1. National Agricultural Imagery Program (NAIP), year: 2019
2. National Landcover Dataset (NLCD), year: 2016
3. Training data- polygons to train the random forest classifier 
(These are generated in the script during the classification stage)

Acknowledgement: The majority of this script is adapted from a script by Ritika et al., and can be found at:
https://github.com/RitikaPrasai/NAIP-aerial-imagery-file/blob/main/Final%20supervised%20classification%20of%20NAIP%20imagery-12%20December.ipynb
'''
Out[1]:
'\nName: Josh Virene\nDate: 9 March, 2023\nClass: NR426: Programming for GIS\nPurpose: The purpose of this lab is to accomplish landcover classification for a study area in Big Sky, Montana using\nsatellite imagery- National Agricultural Imagery Program (NAIP), and the National Landcover Dataset (NLCD)\n\nProgram(s): Jupyter Notebook, Google Earth Engine, Python\nData sources: \n1. National Agricultural Imagery Program (NAIP), year: 2019\n2. National Landcover Dataset (NLCD), year: 2016\n3. Training data- polygons to train the random forest classifier \n(These are generated in the script during the classification stage)\n\nAcknowledgement: The majority of this script is adapted from a script by Ritika et al., and can be found at:\nhttps://github.com/RitikaPrasai/NAIP-aerial-imagery-file/blob/main/Final%20supervised%20classification%20of%20NAIP%20imagery-12%20December.ipynb\n'
In [1]:
# Setup: import modules
try: 
    import ee
    import geemap
    print("Successfully imported the earth engine (ee) and geemap modules")
except: 
    print("Modules not successfully imported")

Successfully imported the earth engine (ee) and geemap modules
In [2]:
# Initialize google earth engine (ee), optional: visualize a blank map 
print("Initializing Earth Engine")
ee.Initialize
print("Successfully initialized")

Map = geemap.Map()

Initializing Earth Engine
Successfully initialized
In [3]:
# Data acquisition: Get NLCD Data
NLCD_Source = 'USGS/NLCD/NLCD2016' # Insert NLCD Imagery pathway
NLCD2016 = ee.Image(NLCD_Source).select('landcover')
try: 
    Map.addLayer(NLCD2016,{},'NLCD_2016',True)
    print("Successfully imported imagery")
except:
    print(f'Failed to import imagery- image source: {NLCD_Source} is not valid')

Successfully imported imagery
In [4]:
Map

Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…
In [5]:
# Select a region of interest over which to classify landcover: 
# Region: Rectangular area south of Big Sky, MT
StudyRegion = ee.Geometry.Rectangle([-111.55543828087393,45.25849168709718,-110.86604619103018, 45.01145511561177])
Map.addLayer(StudyRegion)
print("Created the study region")

Created the study region
In [6]:
# Clip the NLCD Imagery to the study area: ### IMAGE ONE IN LAB WRITE-UP #####
NLCD2016_Clip = NLCD2016.clip(StudyRegion)
NLCD2016_Clip2 = NLCD2016_Clip.geometry()
Map.addLayer(NLCD2016_Clip2,{},'Study_Area') # This is a geometry of the study area- the same as study region
Map.addLayer(NLCD2016_Clip,{},'NLCD_1') # This is NLCD, clipped to the study region
print("Clipped the NLCD imagery to the study region")

Clipped the NLCD imagery to the study region
In [7]:
# Data acquisition: Get the NAIP Imagery- Year: 2019 #### IMAGE TWO IN LAB WRITE-UP #######

NAIP_Source = "USDA/NAIP/DOQQ" # Insert NAIP Imagery pathway
NAIP2019 = ee.ImageCollection(NAIP_Source).filter(ee.Filter.date('2019-01-01', '2019-12-31')).filterBounds(StudyRegion)
NAIP2019_Image = NAIP2019.median()
# .Filterdate specifies the timeframe over which the NAIP imagery was captured
# .filterBounds clips the ImageCollection to the study region defined above
# .median converts the ImageCollection to an image
try: 
    Map.addLayer(NAIP2019_Image,{},'NAIP',True)
    print("Successfully added NAIP Imagery")
except: 
    print(f'Failed to import imagery- image source: {NAIP_Source} is not valid')

Successfully added NAIP Imagery
In [8]:
# Get the class values from the imagery
print("Class and palette information")
print("")
class_values = NLCD2016_Clip.get('landcover_class_values').getInfo()
print(f'Old class values: {class_values}')
print("")
# Create new values 0-19 for the classes: 
number_classes = len(class_values)
new_class_values = list(range(0, number_classes))
print(f'New class values: {new_class_values}')
print("")
# 
palette = NLCD2016_Clip.get('landcover_class_palette').getInfo()
print(f'Palette: {palette}')

Class and palette information

Old class values: [11, 12, 21, 22, 23, 24, 31, 41, 42, 43, 51, 52, 71, 72, 73, 74, 81, 82, 90, 95]

New class values: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]

Palette: ['476ba1', 'd1defa', 'decaca', 'd99482', 'ee0000', 'ab0000', 'b3aea3', '68ab63', '1c6330', 'b5ca8f', 'a68c30', 'ccba7d', 'e3e3c2', 'caca78', '99c247', '78ae94', 'dcd93d', 'ab7028', 'bad9eb', '70a3ba']
In [9]:
# remap the NAIP Imagery, reflecting the new classes
nlcd_2016 = NLCD2016_Clip.remap(class_values, new_class_values).select(['remapped'], ['landcover'])
nlcd_2016 = nlcd_2016.set('landcover_class_values', new_class_values)
nlcd_2016 = nlcd_2016.set('landcover_class_palette', palette)

Map.addLayer(nlcd_2016, {}, 'NLCD_2',True)
print("Successfully added the reclassified NLCD layer")

Successfully added the reclassified NLCD layer
In [10]:
# Create a training dataset by adding random points to the map

num_pixels = 5000 # Changing this number will change the number of points used for the sample data
seed = 0 # The starting point in a random number generating algorithm

# pixel threshold
if num_pixels >20000:
    print(f"WARNING: the script is set to run the model with {num_pixels} pixels, which may cause the map to generate slowly")
else: 
    print(f"The model will run with {num_pixels} pixels")

points = NLCD2016_Clip.sample(
    **{
        'region': NAIP2019.geometry(),
        'scale': 30,
        'numPixels': num_pixels,
        'seed': seed,
        'geometries': True,  # Set this to False to ignore geometries
    }
)

Map.addLayer(points, {}, 'training', False)
print("Successfully created the training dataset")

The model will run with 5000 pixels
Successfully created the training dataset
In [11]:
# Random Forest Classification

# 1) Set imagery bands; these are N - Near Infrared, R- Red, G- Green, B- Blue
classifierbands = classifierbands = ['R', 'G', 'B', 'N']
# 2) Create a training dataset
    # Make a label for the landcover classes
Class = 'landcover'
    # Training
Training = NAIP2019_Image.select(classifierbands).sampleRegions(**{
    'collection' : points,
    'properties' : [Class],
    'scale' : 30
})
print("Successfully Trained the Model")

# 3) # Split the data into a testing dataset and a validation dataset
sample = Training.randomColumn()
split = 0.7 # Split = 0.7 means 70% of points are training data, the remaining 30% are testing data
training_points = sample.filter(ee.Filter.lt('random', split))
validation_points = sample.filter(ee.Filter.gte('random', split))
# This step is necessary for performing accuracy assessment
print("Successfully split the data into training and validation data")
print(f"There are roughly {split*num_pixels:.2f} training points and {(1-split)*num_pixels:.2f} validation points")

Successfully Trained the Model
Successfully split the data into training and validation data
There are roughly 3500.00 training points and 1500.00 validation points
In [12]:
# 4) Run the classifier, create a map from the classified data
trees = 574 # this is the number of trees used in the classifier
RF_classifier = ee.Classifier.smileRandomForest(trees).train(training_points,Class,classifierbands)
Classified_map = NAIP2019_Image.select(classifierbands).classify(RF_classifier)
print("Successfully created a classified map layer")

Successfully created a classified map layer
In [13]:
# Attach the NLCD class values and palette to the map    ### IMAGE THREE IN LAB WRITE-UP ###
landcover_final = Classified_map.set('classification_class_values', class_values)
landcover_final = landcover_final.set('classification_class_palette',palette)
In [14]:
# With the NLCD classes, add the final map to the environment
try: 
    Map.addLayer(landcover_final, {}, 'Land cover')
    print("Successfully created a classified landcover map using the 2019 NAIP Imagery")
except: 
    print("Map failed to generate, check the number of pixels used in the points layer")
    if num_pixels < 20000:
        print("Map failed to generate, check the arguments in Map.addLayer or the random forest classifier (RF_Classifier")
# Attach a legend to the map 
Map.add_legend(builtin_legend='NLCD')

Successfully created a classified landcover map using the 2019 NAIP Imagery
In [15]:
# End of the random forest classification process
In [16]:
# Accuracy Assessment
In [21]:
# Training Data
train_accuracy = RF_classifier.confusionMatrix()
# train_accuracy.accuracy().getInfo()
In [18]:
# Validation Data
validated = validation_points.classify(RF_classifier)
test_accuracy = validated.errorMatrix('landcover', 'classification')
test_accuracy.accuracy().getInfo()
Out[18]:
0.7522522522522522
In [19]:
test_model_accuracy = test_accuracy.accuracy().getInfo()
print(f'Model accuracy value: {test_model_accuracy}')
if test_model_accuracy < 0.7:
    print("Warning, this model has accuracy that is below 70%, check parameters (i.e., number of points, covariates used, and others to increase performance of the model")
else: 
    print("The model accuracy is above 70%, so it is performing reasonably well")

Model accuracy value: 0.7522522522522522
The model accuracy is above 70%, so it is performing reasonably well
In [20]:
Map

Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…
In [ ]: