Tuesday 14 March 2023

Simple classification using Google Earth Engine Python API!

Hello,

I like documenting some of code in my blog. Today I would like to share a simple classification using the python API of Google Earth Egnine. The results are not great but the approach is working.

 

I divided the code into two parts. In the first part, I define three polygons for deriving testing data for three classes: Urban, Green and Water and then classify a Landsat image. For the 2nd part, the visualisation, I modified the folium example provided here: https://colab.research.google.com/github/giswqs/qgis-earthengine-examples/blob/master/Folium/ee-api-folium-setup.ipynb#scrollTo=VP33EU7hJBq8 (Accessed Feb 2023). 


 First part:



import sys
import csv
import ee
ee.Authenticate()
ee.Initialize()

# Get the polygon defining Washington 
countyData = ee.FeatureCollection('TIGER/2018/Counties')
countyConnect = countyData.filter(ee.Filter.eq('STATEEF','11'))
countyCpmmectDiss = countyConnect.union(100)


# Load US county dataset.
countyData = ee.FeatureCollection('TIGER/2018/Counties')
# Filter the counties that are in Connecticut (more on filters later).
washington = countyData.filter(ee.Filter.eq('STATEFP', '11'))

Urban = ee.FeatureCollection(
        [ee.Feature(
            ee.Geometry.Polygon(
                [[[-77.0535926206179, 38.89903486157564],
                  [-77.0535926206179, 38.895561277858654],
                  [-77.04672616553978, 38.895427675091845],
                  [-77.04655450416283, 38.899569243990044]]]),
            {
              "class": 0,
              "system:index": "0"
            })])

Water = ee.FeatureCollection(
        [ee.Feature(
            ee.Geometry.Polygon(
                [[[-77.05302814295763, 38.886078313198695],
                  [-77.05611804774279, 38.88340585335565],
                  [-77.04890826991075, 38.880198768850136]]]),
            {
              "class": 1,
              "system:index": "0"
            })])


Green = ee.FeatureCollection(
        [ee.Feature(
            #pGreenPoint,
            ee.Geometry.Polygon(
                [[[-77.03282, 38.87528],
                  [-77.02759, 38.87545],
                  [-77.0296 , 38.87254]]]),
            {
              "class": 2,
              "system:index": "0"
            })])

alltdata = Green.merge(Water)
alltdata = alltdata.merge(Urban)

bands = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7', 'ST_B10']


# Get a Landsat image and select Washington to avoid processing too many images
image = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') \
  .filterBounds(washington) \
  .filterDate('2021-03-01', '2021-07-01')\
  .median()


label = "class"
training = image.select(bands).sampleRegions(
  collection = alltdata,
  properties = [label],
  scale      = 10
)


# Train classifier with 20 trees
classifier = ee.Classifier.smileCart(20).train(training, label,bands)
#print(classifier.confusionMatrix().accuracy())
# Classify the landsat images 
classified = image.select(bands).classify(classifier)

print("TestFeatureVector: Processing Complete")
 

Second part:

#visualisation
import folium
import pandas as pd
import numpy as np
from folium import plugins

# Add custom basemaps to folium
basemaps = {
    'Google Maps': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Maps',
        overlay = True,
        control = True
    ),
    'Google Satellite': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Satellite',
        overlay = True,
        control = True
    ),
    'Google Terrain': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=p&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Terrain',
        overlay = True,
        control = True
    ),
    'Google Satellite Hybrid': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Satellite',
        overlay = True,
        control = True
    ),
    'Esri Satellite': folium.TileLayer(
        tiles = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
        attr = 'Esri',
        name = 'Esri Satellite',
        overlay = True,
        control = True
    )
}

# Define a method for displaying Earth Engine image tiles on a folium map.
def add_ee_layer(self, ee_object, vis_params, name):
    
    try:    
        # display ee.Image()
        if isinstance(ee_object, ee.image.Image):    
            map_id_dict = ee.Image(ee_object).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
            ).add_to(self)
        # display ee.ImageCollection()
        elif isinstance(ee_object, ee.imagecollection.ImageCollection):    
            ee_object_new = ee_object.mosaic()
            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
            ).add_to(self)
        # display ee.Geometry()
        elif isinstance(ee_object, ee.geometry.Geometry):    
            folium.GeoJson(
            data = ee_object.getInfo(),
            name = name,
            overlay = True,
            control = True
        ).add_to(self)
        # display ee.FeatureCollection()
        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):  
            ee_object_new = ee.Image().paint(ee_object, 0, 2)
            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
        ).add_to(self)

    except:
        print("Could not display {}".format(name))
    
# Add EE drawing method to folium.
folium.Map.add_ee_layer = add_ee_layer


my_map = folium.Map(location=[-8571600, 4579425], zoom_start=5, height=400)
# Add custom basemaps
basemaps['Google Maps'].add_to(my_map)
basemaps['Google Satellite Hybrid'].add_to(my_map)
# Add Land Mask to the map

# Set visualization parameters.
classVis = {
  'min': 0,
  'max': 10,
  'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']}



#my_map.add_ee_layer(s2median , {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 2500, 'gamma': 1.1}, 'median'      )
#my_map.add_ee_layer(ppBuffer, {}, 'ppBuffer'      )
my_map.add_ee_layer(Urban      , {}, 'Urban'        )
my_map.add_ee_layer(Water      , {}, 'Water'        )
my_map.add_ee_layer(Green      , {}, 'Green'        )

my_map.add_ee_layer(classified,classVis,'classified')



# Add a layer control panel to the map.
my_map.add_child(folium.LayerControl())
plugins.Fullscreen().add_to(my_map)

# Add a layer control panel to the map.
my_map.add_child(folium.LayerControl())

# Add fullscreen button
plugins.Fullscreen().add_to(my_map)

# Display the map.
display(my_map)
 


Here is an example of the output: 

For the first part, I used this javascript tutorial to start with (https://www.youtube.com/watch?v=Md5_qkT-vC8&ab_channel=NoelGorelick), which was very useful but the final code is actually very differnt. I predominanly kept the study area, which is Washington.

 I hope you find this tutorial useful to help you start working with your own data. :)  


Acknowledgements: This work is done as part of my postdoctdoctoral research at University of Cambridge, funded by UKRI Future Leader's fellowship of Dr Emily Lines (FLF LCAG/439).