Alaska_2024_Part3.ipynb

Kernel: Python 3

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

# %pip install geemap pygis cartopy solara

Import libraries

import ee
import geemap

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.

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

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

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

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

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

Unsupervised classification

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()

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

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

geemap.ee_to_df(training.limit(5))

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

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

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

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

Supervised classification

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

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

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

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

geemap.ee_to_df(features.limit(5))

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

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

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}
)

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

Accuracy assessment

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

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()

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

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

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

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

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

trainAccuracy.accuracy().getInfo()

trainAccuracy.kappa().getInfo()

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

validationAccuracy.accuracy()

validationAccuracy.producersAccuracy().getInfo()

validationAccuracy.consumersAccuracy().getInfo()

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())

imgClassified = img.classify(trainedClassifier)

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")

Create and export maps

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

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

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

from geemap import cartoee
import matplotlib.pyplot as plt

Plotting single-band images

import matplotlib.pyplot as plt
from geemap import cartoee

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}

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()

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()

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

Plotting multi-band images

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

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()

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

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

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

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()

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

Custom projections

import cartopy.crs as ccrs

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()

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()

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()

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()

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")

Building interactive web apps

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

# 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.