Image classification with GEE and scikit-learn

So far, we have come a long way. We have learned to work with both, vector and raster data. We have used this knowledge to e.g., extract for specific geo-locations raster-values, to use these data for parameterizing a machine-learning classifier, and to do a subsequent prediction. Now, what this small tutorial does is show you how to do these steps using a combination of functions from the packages scikit-learn and geemap. We will exemplify this on a - admittedly not great - example. Specifically, we will collect some training data from an existing land cover dataset, and predict the same land cover categories on spectra-temporal metrics (STM). Obviously, one would expect the results to be similar, but for the sake of this tutorial, we keep it this way.

We now go step by step over the process with some example code. First, we load some packages and initialize GEE

from osgeo import ogr
from tqdm.notebook import tqdm
import json
import pandas as pd

import ee
import geemap.foliumap as geemap
try:
    ee.Initialize()
except Exception as e:
    ee.Authenticate()
    ee.Initialize()

Load vector data¶

Open them using the `ogr` library¶

First, we will load two vector files, both are located on my personal drive:

One file does only contain one feature, with a squared geometry. This will be our region of interest (ROI), i.e., the boundaries of our classification
A set of points that will serve as our training locations. This is not an uncommon scenario: often you have/receive some shapefiles in the form of points, which you use for your work.

roiSHP = ogr.Open("D:/OneDrive - Conservation Biogeography Lab/_TEACHING/__Classes-Modules_HUB/M8_Geoprocessing-with-python/10_Classification-and-LargeArea/01_ROI_shape.shp")
roiLYR = roiSHP.GetLayer()
roiFEAT = roiLYR[0]
roiGEOM = roiFEAT.geometry()
print(roiGEOM)

POLYGON ((13.2134494446947 52.4151979479578,13.2134494446947 52.4873242002239,13.3018503384978 52.4873242002239,13.3018503384978 52.4151979479578,13.2134494446947 52.4151979479578))

ptsSHP = ogr.Open("D:/OneDrive - Conservation Biogeography Lab/_TEACHING/__Classes-Modules_HUB/M8_Geoprocessing-with-python/10_Classification-and-LargeArea/02_RandomPoints_1000.shp")
ptsLYR = ptsSHP.GetLayer()
ptsLYR.GetFeatureCount()

1000

Convert items into `ee.Geometry()` and `ee.FeatureCollection()`¶

Next, we will convert these items into GEE-Objects. The ROI we convert it into an ee.Geometry(), the points into an ee.FeatureCollection(). Both items we have already learned how to do (check the intro chapter on GEE of this course

geomJSON = json.loads(roiGEOM.ExportToJson())
geomCOORD = geomJSON['coordinates']
geomEE = ee.Geometry.Polygon(coords=geomCOORD)
geomEE

# Create an empty list that we populate
featList = []
# Iterate over all features in `ptsLYR`
for feat in tqdm(ptsLYR):
    pid = feat.GetField('Id') # Get the ID value
    # Get the geometry
    g = feat.geometry()
    gjson = json.loads(g.ExportToJson())
    gcoord = gjson['coordinates']
    gEE = ee.Geometry.Point(coords=gcoord)
    # Create ee.Feature(), append to list
    eeFeat = ee.Feature(gEE, {"UniqueID": pid})
    featList.append(eeFeat)
# Convert the list into ee.FeatureCollection()
fc = ee.FeatureCollection(ee.List(featList))
fc

Prepare STM¶

From our chapters on image classifications and image workflows in GEE we learned what STMs are and how we can use them. Here, we repeat many of these steps, while at the same time providing some smaller functions that you can incorporate in your library.

Specifically, we take all Landsat 8 and Landsat 9 imagery of the year 2023 with less than 70% cloud cover. We merge the images of these two satellites into one collection, mask and scale the data, and calculate basic STMs.

def ScaleMask_l89(img):
    ## Scale the bands
    refl   = img.select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']).multiply(0.0000275).add(-0.2)
    img_SR = ee.Image(refl).addBands(img.select(['QA_PIXEL']))
    ## mask cloud
    # Get the QA band and apply the bits
    qa = ee.Image(img_SR).select(['QA_PIXEL'])
    dilatedCloud = qa.bitwiseAnd(1 << 1)
    cirrusMask   = qa.bitwiseAnd(1 << 2)
    cloud        = qa.bitwiseAnd(1 << 3)
    cloudShadow  = qa.bitwiseAnd(1 << 4)
    snow         = qa.bitwiseAnd(1 << 5)
    water        = qa.bitwiseAnd(1 << 7)
    # Apply the bits
    clear = dilatedCloud.Or(cloud).Or(cloudShadow).eq(0)
    cfmask = (clear.eq(0).where(water.gt(0), 1)
        .where(snow.gt(0), 3)
        .where(dilatedCloud.Or(cloud).gt(0), 4)
        .where(cloudShadow.gt(0), 2)
        .rename('cfmask'))
    img_SR_cl = ee.Image(img_SR).select(['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']).mask(cfmask.lte(2))
    return img_SR_cl

# Get the collections, merge them
l8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")\
    .filterDate('2023-03-01', '2023-11-30').filter(ee.Filter.lt('CLOUD_COVER', 70)).filterBounds(geomEE)
l9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")\
    .filterDate('2022-03-01', '2023-11-30').filter(ee.Filter.lt('CLOUD_COVER', 70)).filterBounds(geomEE)
lAll = l8.merge(l9)
# Apply masking and filtering
lALL_masked = lAll.map(ScaleMask_l89)
# Define reducers
mean = ee.Reducer.mean().unweighted()
sd = ee.Reducer.stdDev().unweighted()
percentiles = ee.Reducer.percentile([10, 25, 50, 75, 90]).unweighted()
allMetrics = mean.combine(sd, sharedInputs=True).combine(percentiles, sharedInputs=True)
# Apply reducers
stm = lALL_masked.reduce(allMetrics)

We can check the outcome by having a quick look at the band names. I do that here by getting the band names for all metrics calculated for the blue band:

bn = stm.bandNames().getInfo()
bn[0:7]

['SR_B2_mean',
 'SR_B2_stdDev',
 'SR_B2_p10',
 'SR_B2_p25',
 'SR_B2_p50',
 'SR_B2_p75',
 'SR_B2_p90']

Prepare Land Cover variables¶

That seemed to have worked well. In essence, we have now our independent variables prepared. As dependent variable - again: thematically not the most suitable case - we extract the label of the ESA World Cover, and add it as an additional band to the SMT image:

# Is an image collection from which we take the first image
esaWC21 = ee.ImageCollection("ESA/WorldCover/v200").select(['Map']).toBands()
# Add to the median composite
allRas = esaWC21.addBands(stm)

Extract values at the point locations¶

Now that we have everything ready, we can extract the data for our point locations and store them on our computers. We use the function.sampleRegions() to do that. This function works well over smaller areas; for larger areas we need some workaround (see next chapter).

vals = allRas.sampleRegions(collection=fc, properties=['UniqueID'], scale=30, tileScale=4, geometries=False).getInfo()
#vals

the code box below then allows you to convert the json dictionary into a pandas data frame. It is probably not the smartest option/solution, but it does the job. Feel free to write us with a better option:

flag = 0
featureValues = vals['features']
for f in featureValues:         
 # Chekc here always the right order of the indices
    prop = f['properties']
    if flag == 0:
        out_pd = pd.DataFrame(prop, index=[0])
        flag = 1
    else:
        out_pd = pd.concat([out_pd, pd.DataFrame(prop, index=[0])])
out_pd.head(10)

Parameterization in `scikit-learn` and conversion into GEE-classifier¶

The data are ready. We can now move on with the parameterization and the prediction. We will use some pretty cool tools from the geemap, including the upload of trained trees beack into GEE. We use a random forest classifier to parameterize a model, similar to what we did in the machine learning session session.

from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestClassifier

# Create the two arrays
y = out_pd['2021_Map']
X = out_pd.drop(['2021_Map', 'UniqueID'], axis=1)
# Create a classifier instance
RF = RandomForestClassifier(random_state=1)
param_grid = {'n_estimators': [10, 25, 50, 75, 100],# 'numberOfTrees
              'max_features': [5, 10, 'sqrt'],# 'variablesPerSplit'
              'min_samples_split': [5, 10, 15], # minLeafPopulation
              'max_depth': [5, 10, 20]} # maxNodes
RF_cv = GridSearchCV(RF, param_grid=param_grid, cv=3, refit=True, n_jobs=4).fit(X, y)
best_mod = RF_cv.best_estimator_

D:\miniforge3\envs\py310\lib\site-packages\sklearn\model_selection\_split.py:725: UserWarning: The least populated class in y has only 2 members, which is less than n_splits=3.
  warnings.warn(

Now, we apply a pretty cool part, implemented in the geemap package. What the two functions below are doing is that they (1) convert the random forest model into strings (i.e., the trees), (2) convert it into a feature collection, before (3) converting it into a ee.Classifer. For very large models with very long variable names the maximum size of 10MB per model might not be sufficient, but for the problem here it is fine.

from geemap import ml
feature_names = X.columns.tolist() # get the names of the variables
best_mod_str = ml.rf_to_strings(best_mod, feature_names)
ee_classifier = ml.strings_to_classifier(best_mod_str)

ee_classifier

Prediction in GEE and export¶

Now we have (1) the STM sitting in form of a container in GEE and (2) a model which we have trained. All what is needed is to bring these two together through a prediction. This is very easily done through one line of code. Remember, that until you hit getInfo() or export the dataset to your drive or into an asset, GEE does not a tiny thing of work. Some examples for exporting data were already provided in the GEE raster chapter, so we only repeat the code here without muhc explanation.

classification = stm.classify(ee_classifier)

task = ee.batch.Export.image.toDrive(
    image=classification,
    description='Berlin_Class_Export',
    folder='your_folder',
    fileNamePrefix= 'Berlin_Class_Export',
    region=geomEE,
    scale=30)
task.start()

palette = {'min': 10,'max': 100,
           'palette': ['006400','ffbb22','ffff4c','f096ff','fa0000','b4b4b4','f0f0f0','0064c8','0096a0','00cf75','fae6a0']}
ctr = geomEE.centroid().getInfo()['coordinates']
Map = geemap.Map(center=(ctr[1], ctr[0]), zoom=11)

Map.addLayer(stm, {"bands": ['SR_B7_mean', 'SR_B5_mean', 'SR_B3_mean'], "min": 0.05, "max": 0.55, "gamma": 1.5}, 'image')
Map.addLayer(classification, palette, 'classification')

Map

And that is it. We essentially have used GEE as a pure processing and data engine. The advantages of this approach are manyfold:

We can have now basically a classification as a pre-processing step before doing some more advanced techniques offline which we cannot in GEE, such as image segmentation. Think about, for example, a probability image [of a class] which you want to segment using a random walker segmentation. In GEE alone you cannot do this, but you also don’t want to (a) download all the images, (b) run the pre-processing, and (c) calculate the probabilities before running the segmentation.
You can also think about everything else until the model prediction - so, harvesting GEE for its data to run all sorts of tests (e.g., different STMs, different images included in a growing season, etc.
Sometimes you want to keep you workflow, and maybe even dynamically adjust your training sample, e.g., by using different buffer sizes. You can do all the vector-based manipulation offline and only use GEE for the bits of data extraction and processing.

Image classification with GEE and scikit-learn

Load vector data¶

Open them using the ogr library¶

Convert items into ee.Geometry() and ee.FeatureCollection()¶

Prepare STM¶

Prepare Land Cover variables¶

Extract values at the point locations¶

Parameterization in scikit-learn and conversion into GEE-classifier¶

Prediction in GEE and export¶

Further reading¶

Open them using the `ogr` library¶

Convert items into `ee.Geometry()` and `ee.FeatureCollection()`¶

Parameterization in `scikit-learn` and conversion into GEE-classifier¶