Geospatial Cloud Computing with the GEE Python API - Part 3

Notebook: https://geemap.org/workshops/Alaska_2024_Part3
Earth Engine: https://earthengine.google.com
Geemap: https://geemap.org

Introduction¶

This notebook contains the materials for the third part of the workshop Geospatial Cloud Computing with the GEE Python API at the University of Alaska Fairbanks.

This workshop provides an introduction to cloud-based geospatial analysis using the Earth Engine Python API. Attendees will learn the basics of Earth Engine data types and how to visualize, analyze, and export Earth Engine data in a Jupyter environment with geemap. In addition, attendees will learn how to develop and deploy interactive Earth Engine web apps with Python. Through practical examples and hands-on exercises, attendees will enhance their learning experience. During each hands-on session, attendees will walk through Jupyter Notebook examples on Google Colab with the instructors. At the end of each session, they will complete a hands-on exercise to apply the knowledge they have learned.

Agenda¶

The workshop is divided into three parts. The third part will cover the following topics:

Image Classification (focused on land cover in Alaska)
Accuracy assessment
Create and export maps
Building interactive web apps

Prerequisites¶

To use geemap and the Earth Engine Python API, you must register for an Earth Engine account and follow the instructions here to create a Cloud Project. Earth Engine is free for noncommercial and research use. To test whether you can use authenticate the Earth Engine Python API, please run this notebook on Google Colab.

Technical requirements¶

Install packages¶

conda create -n gee python=3.11
conda activate gee
conda install -c conda-forge mamba
mamba install -c conda-forge pygis cartopy solara

In [ ]:

# %pip install geemap pygis cartopy solara

Import libraries¶

In [ ]:

import ee
import geemap

In [ ]:

geemap.ee_initialize()

Image Classification¶

North American Land Change Monitoring System (NALCMS)¶

The 2020 North American Land Cover 30-meter dataset was produced as part of the North American Land Change Monitoring System (NALCMS), a trilateral effort between Natural Resources Canada, the United States Geological Survey, and three Mexican organizations.

In [ ]:

m = geemap.Map(center=[40, -100], zoom=4, height=700)
m.add_basemap("Esri.WorldImagery", False)
landcover = ee.Image("USGS/NLCD_RELEASES/2020_REL/NALCMS")

states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.eq("NAME", "Alaska"))

legend_dict = {
    "1 Temperate forest": "033e00",
    "2 Sub-polar forest": "939b71",
    "3 Tropical evergreen forest": "196d12",
    "4 Tropical deciduous forest": "1fab01",
    "5 Temperate deciduous forest": "5b725c",
    "6 Mixed forest": "6b7d2c",
    "7 Tropical shrubland": "b29d29",
    "8 Temperate shrubland": "b48833",
    "9 Tropical grassland": "e9da5d",
    "10 Temperate grassland": "e0cd88",
    "11 Sub-polar shrubland": "a07451",
    "12 Sub-polar grassland": "bad292",
    "13 Sub-polar barren": "3f8970",
    "14 Wetland": "6ca289",
    "15 Cropland": "e6ad6a",
    "16 Barren land": "a9abae",
    "17 Urban and built-up": "db2126",
    "18 Water": "4c73a1",
    "19 Snow and ice ": "fff7fe",
}

palette = list(legend_dict.values())
vis_params = {"palette": palette, "min": 1, "max": 19}

m.add_layer(landcover, vis_params, "NALCMS Land Cover")
m.add_layer(fc, {}, "Alaska", False)
m.center_object(fc, 4)
m.add_legend(title="Land Cover Type", legend_dict=legend_dict)
m

In [ ]:

fig = geemap.image_histogram(
    landcover, region=fc, x_label="Land cover type", y_label="Pixels"
)
fig

In [ ]:

values = list(fig.data[0]["x"])
values

In [ ]:

classes = [list(legend_dict.keys())[int(value) - 1] for value in values]
classes

In [ ]:

fig.update_xaxes(tickvals=values, ticktext=classes)
fig

Unsupervised classification¶

In [ ]:

region = ee.Geometry.BBox(-149.352, 64.5532, -147.0976, 65.1277)
collection = geemap.landsat_timeseries(
    region, start_year=2021, end_year=2021, start_date="06-01", end_date="09-01"
)
image = collection.first()

In [ ]:

m = geemap.Map()
m.add_basemap("Esri.WorldImagery")
vis_params = {"min": 0, "max": 0.3, "bands": ["NIR", "Red", "Green"]}
m.add_layer(image, vis_params, "Image")
m.add_layer(landcover.clip(region), {}, "NALCMS land cover", False)
# m.add_legend(title='Land Cover Type', legend_dict=legend_dict, layer_name='NALCMS land cover')
m.center_object(region, 9)
m

In [ ]:

training = image.sample(
    **{
        "region": region,
        "scale": 150,
        "numPixels": 5000,
        "seed": 1,
        "geometries": True,
    }
)
m.add_layer(training, {}, "Training samples")

In [ ]:

geemap.ee_to_df(training.limit(5))

In [ ]:

clusterer = ee.Clusterer.wekaXMeans(minClusters=3, maxClusters=6).train(training)

In [ ]:

result = image.cluster(clusterer)
m.add_layer(result.randomVisualizer(), {}, "Clusters")
m.add("layer_manager")
m

In [ ]:

geemap.image_histogram(landcover, region=region, scale=30)

In [ ]:

geemap.image_histogram(result, region=region, scale=300)

Supervised classification¶

In [ ]:

m = geemap.Map()
m.add_basemap("Esri.WorldImagery")
vis_params = {"min": 0, "max": 0.3, "bands": ["NIR", "Red", "Green"]}
m.add_layer(image, vis_params, "Image")
m.add_layer(landcover.clip(region), {}, "NALCMS land cover", False)
# m.add_legend(title='Land Cover Type', legend_dict=legend_dict, layer_name='NALCMS land cover')
m.center_object(region, 9)
m

In [ ]:

points = landcover.sample(
    **{
        "region": region,
        "scale": 150,
        "numPixels": 5000,
        "seed": 1,
        "geometries": True,
    }
)

m.add_layer(points, {}, "Training", False)

In [ ]:

label = "landcover"
features = image.sampleRegions(
    **{"collection": points, "properties": [label], "scale": 150}
)

In [ ]:

geemap.ee_to_df(features.limit(5))

In [ ]:

bands = image.bandNames().getInfo()
params = {
    "features": features,
    "classProperty": label,
    "inputProperties": bands,
}
classifier = ee.Classifier.smileCart(maxNodes=None).train(**params)

In [ ]:

classified = image.classify(classifier).rename("landcover")
m.add_layer(classified.randomVisualizer(), {}, "Classified")
m

In [ ]:

class_values = list(range(1, 20))
class_palette = list(legend_dict.values())
classified = classified.set(
    {"landcover_class_values": class_values, "landcover_class_palette": class_palette}
)

In [ ]:

m.add_layer(classified, {}, "Land cover")
m.add_legend(title="Land cover type", builtin_legend="NLCD")
m

Accuracy assessment¶

In [ ]:

m = geemap.Map()
point = ee.Geometry.Point([-122.4439, 37.7538])

img = (
    ee.ImageCollection("COPERNICUS/S2_SR")
    .filterBounds(point)
    .filterDate("2020-01-01", "2021-01-01")
    .sort("CLOUDY_PIXEL_PERCENTAGE")
    .first()
    .select("B.*")
)

vis_params = {"min": 100, "max": 3500, "bands": ["B11", "B8", "B3"]}

m.centerObject(point, 9)
m.add_layer(img, vis_params, "Sentinel-2")
m

In [ ]:

lc = ee.Image("ESA/WorldCover/v100/2020")
classValues = [10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100]
remapValues = ee.List.sequence(0, 10)
label = "lc"
lc = lc.remap(classValues, remapValues).rename(label).toByte()

In [ ]:

sample = img.addBands(lc).stratifiedSample(
    **{
        "numPoints": 100,
        "classBand": label,
        "region": img.geometry(),
        "scale": 10,
        "geometries": True,
    }
)

In [ ]:

sample = sample.randomColumn()
trainingSample = sample.filter("random <= 0.8")
validationSample = sample.filter("random > 0.8")

In [ ]:

trainedClassifier = ee.Classifier.smileRandomForest(numberOfTrees=10).train(
    **{
        "features": trainingSample,
        "classProperty": label,
        "inputProperties": img.bandNames(),
    }
)

In [ ]:

# print('Results of trained classifier', trainedClassifier.explain().getInfo())

In [ ]:

trainAccuracy = trainedClassifier.confusionMatrix()
trainAccuracy.getInfo()

In [ ]:

trainAccuracy.accuracy().getInfo()

In [ ]:

trainAccuracy.kappa().getInfo()

In [ ]:

validationSample = validationSample.classify(trainedClassifier)
validationAccuracy = validationSample.errorMatrix(label, "classification")
validationAccuracy.getInfo()

In [ ]:

validationAccuracy.accuracy()

In [ ]:

validationAccuracy.producersAccuracy().getInfo()

In [ ]:

validationAccuracy.consumersAccuracy().getInfo()

In [ ]:

import csv

with open("training.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerows(trainAccuracy.getInfo())

with open("validation.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerows(validationAccuracy.getInfo())

In [ ]:

imgClassified = img.classify(trainedClassifier)

In [ ]:

classVis = {
    "min": 0,
    "max": 10,
    "palette": [
        "006400",
        "ffbb22",
        "ffff4c",
        "f096ff",
        "fa0000",
        "b4b4b4",
        "f0f0f0",
        "0064c8",
        "0096a0",
        "00cf75",
        "fae6a0",
    ],
}
m.add_layer(lc, classVis, "ESA Land Cover", False)
m.add_layer(imgClassified, classVis, "Classified")
m.add_layer(trainingSample, {"color": "black"}, "Training sample")
m.add_layer(validationSample, {"color": "white"}, "Validation sample")
m.add_legend(title="Land Cover Type", builtin_legend="ESA_WorldCover")
m.centerObject(img)
m

Exercise 1 - Unsupervised classification¶

Perform an unsupervised classification of a Sentinel-2 imagery for your preferred area. Relevant Earth Engine assets:

ee.ImageCollection("COPERNICUS/S2_HARMONIZED")

In [ ]:

Create and export maps¶

In [ ]:

m = geemap.Map(center=(41.0462, -109.7424), zoom=6)

dem = ee.Image("USGS/SRTMGL1_003")
landsat7 = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003").select(
    ["B1", "B2", "B3", "B4", "B5", "B7"]
)

vis_params = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}

m.add_layer(
    landsat7,
    {"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 2},
    "landsat",
)
m.add_layer(dem, vis_params, "dem", True, 1)
m

In [ ]:

try:
    m.layer_to_image("dem", output="dem.tif", crs="EPSG:3857", region=None, scale=None)
    m.layer_to_image("dem", output="dem.jpg", scale=500)
    geemap.show_image("dem.jpg")
except:
    pass

In [ ]:

try:
    m.layer_to_image("landsat", output="landsat.tif")
    geemap.geotiff_to_image("landsat.tif", output="landsat.jpg")
    geemap.show_image("landsat.jpg")
except:
    pass

In [ ]:

from geemap import cartoee
import matplotlib.pyplot as plt

Plotting single-band images¶

In [ ]:

import matplotlib.pyplot as plt
from geemap import cartoee

In [ ]:

srtm = ee.Image("CGIAR/SRTM90_V4")

# define bounding box [east, south, west, north] to request data
region = [180, -60, -180, 85]
vis = {"min": 0, "max": 3000}

In [ ]:

fig = plt.figure(figsize=(15, 9))

# use cartoee to get a map
ax = cartoee.get_map(srtm, region=region, vis_params=vis)

# add a color bar to the map using the visualization params we passed to the map
cartoee.add_colorbar(
    ax, vis, loc="bottom", label="Elevation (m)", orientation="horizontal"
)

# add grid lines to the map at a specified interval
cartoee.add_gridlines(ax, interval=[60, 30], linestyle=":")

# add coastlines using the cartopy api
ax.coastlines(color="red")

plt.show()

In [ ]:

fig = plt.figure(figsize=(15, 7))

cmap = "terrain"

ax = cartoee.get_map(srtm, region=region, vis_params=vis, cmap=cmap)
cartoee.add_colorbar(
    ax, vis, cmap=cmap, loc="right", label="Elevation (m)", orientation="vertical"
)

cartoee.add_gridlines(ax, interval=[60, 30], linestyle="--")
ax.coastlines(color="red")
ax.set_title(label="Global Elevation Map", fontsize=15)

plt.show()

In [ ]:

cartoee.savefig(fig, fname="srtm.jpg", dpi=300, bbox_inches="tight")

Plotting multi-band images¶

In [ ]:

image = ee.Image("LANDSAT/LC08/C01/T1_SR/LC08_044034_20140318")
vis = {"bands": ["B5", "B4", "B3"], "min": 0, "max": 5000, "gamma": 1.3}

In [ ]:

fig = plt.figure(figsize=(15, 10))

ax = cartoee.get_map(image, vis_params=vis)
cartoee.pad_view(ax)
cartoee.add_gridlines(ax, interval=0.5, xtick_rotation=0, linestyle=":")
ax.coastlines(color="yellow")

plt.show()

In [ ]:

fig = plt.figure(figsize=(15, 10))

region = [-121.8025, 37.3458, -122.6265, 37.9178]
ax = cartoee.get_map(image, vis_params=vis, region=region)
cartoee.add_gridlines(ax, interval=0.15, xtick_rotation=0, linestyle=":")
ax.coastlines(color="yellow")

plt.show()

Using custom projections¶

The PlateCarree projection¶

In [ ]:

ocean = (
    ee.ImageCollection("NASA/OCEANDATA/MODIS-Terra/L3SMI")
    .filter(ee.Filter.date("2018-01-01", "2018-03-01"))
    .median()
    .select(["sst"], ["SST"])
)

In [ ]:

visualization = {"bands": "SST", "min": -2, "max": 30}
bbox = [180, -88, -180, 88]

In [ ]:

fig = plt.figure(figsize=(15, 10))

ax = cartoee.get_map(ocean, cmap="plasma", vis_params=visualization, region=bbox)
cb = cartoee.add_colorbar(ax, vis_params=visualization, loc="right", cmap="plasma")

ax.set_title(label="Sea Surface Temperature", fontsize=15)

ax.coastlines()
plt.show()

In [ ]:

cartoee.savefig(fig, "SST.jpg", dpi=300)

Custom projections¶

In [ ]:

import cartopy.crs as ccrs

In [ ]:

fig = plt.figure(figsize=(15, 10))

projection = ccrs.Mollweide(central_longitude=-180)
ax = cartoee.get_map(
    ocean, vis_params=visualization, region=bbox, cmap="plasma", proj=projection
)
cb = cartoee.add_colorbar(
    ax, vis_params=visualization, loc="bottom", cmap="plasma", orientation="horizontal"
)
ax.set_title("Mollweide projection")
ax.coastlines()

plt.show()

In [ ]:

fig = plt.figure(figsize=(15, 10))

projection = ccrs.Robinson(central_longitude=-180)
ax = cartoee.get_map(
    ocean, vis_params=visualization, region=bbox, cmap="plasma", proj=projection
)
cb = cartoee.add_colorbar(
    ax, vis_params=visualization, loc="bottom", cmap="plasma", orientation="horizontal"
)
ax.set_title("Robinson projection")
ax.coastlines()

plt.show()

In [ ]:

fig = plt.figure(figsize=(15, 10))

projection = ccrs.InterruptedGoodeHomolosine(central_longitude=-180)
ax = cartoee.get_map(
    ocean, vis_params=visualization, region=bbox, cmap="plasma", proj=projection
)
cb = cartoee.add_colorbar(
    ax, vis_params=visualization, loc="bottom", cmap="plasma", orientation="horizontal"
)
ax.set_title("Goode homolosine projection")
ax.coastlines()

plt.show()

In [ ]:

fig = plt.figure(figsize=(15, 10))

projection = ccrs.EqualEarth(central_longitude=-180)
ax = cartoee.get_map(
    ocean, vis_params=visualization, region=bbox, cmap="plasma", proj=projection
)
cb = cartoee.add_colorbar(
    ax, vis_params=visualization, loc="right", cmap="plasma", orientation="vertical"
)
ax.set_title("Equal Earth projection")
ax.coastlines()

plt.show()

In [ ]:

fig = plt.figure(figsize=(11, 10))

projection = ccrs.Orthographic(-130, -10)
ax = cartoee.get_map(
    ocean, vis_params=visualization, region=bbox, cmap="plasma", proj=projection
)
cb = cartoee.add_colorbar(
    ax, vis_params=visualization, loc="right", cmap="plasma", orientation="vertical"
)
ax.set_title("Orographic projection")
ax.coastlines()

plt.show()

Exercise 2 - Creating NDVI maps¶

Create and export a global NDVI map using MODIS data. Relevant Earth Engine assets:

ee.ImageCollection("MODIS/MCD43A4_006_NDVI")

In [ ]:

Building interactive web apps¶

In [ ]:

import ee
import geemap
import solara


class Map(geemap.Map):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.add_ee_data()

    def add_ee_data(self):
        years = ["2001", "2004", "2006", "2008", "2011", "2013", "2016", "2019"]

        def getNLCD(year):
            dataset = ee.ImageCollection("USGS/NLCD_RELEASES/2019_REL/NLCD")
            nlcd = dataset.filter(ee.Filter.eq("system:index", year)).first()
            landcover = nlcd.select("landcover")
            return landcover

        collection = ee.ImageCollection(ee.List(years).map(lambda year: getNLCD(year)))
        labels = [f"NLCD {year}" for year in years]
        self.ts_inspector(
            left_ts=collection,
            right_ts=collection,
            left_names=labels,
            right_names=labels,
        )
        self.add_legend(
            title="NLCD Land Cover Type",
            builtin_legend="NLCD",
            height="460px",
            add_header=False,
        )


@solara.component
def Page():
    with solara.Column(style={"min-width": "500px"}):
        Map.element(
            center=[40, -100],
            zoom=4,
            height="800px",
        )


Page()

solara run ./pages

In [ ]:

# geemap.get_ee_token()

Exercise 3 - Deploying an Earth Engine app on Hugging Face.¶

Follow the instructions here to deploy an Earth Engine web app on Hugging Face.

In [ ]: