Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
giswqs
GitHub Repository: giswqs/geemap
 Path: blob/master/docs/workshops/SciPy_2023.ipynb
2313 views
Kernel: Python 3

image

An Introduction to Cloud-Based Geospatial Analysis with Earth Engine and Geemap

Introduction (10 mins)

This notebook is for the tutorial presented at the SciPy 2023 Conference by Qiusheng Wu and Steve Greenberg. Check out this link for more information about the tutorial.

Abstract

This tutorial is an introduction to cloud-based geospatial analysis with Earth Engine and the geemap Python package. We will cover the basics of Earth Engine data types and how to visualize, analyze, and export Earth Engine data in a Jupyter environment using geemap. We will also demonstrate how to develop and deploy interactive Earth Engine web apps. Throughout the session, practical examples and hands-on exercises will be provided to enhance learning. The attendees should have a basic understanding of Python and Jupyter Notebooks. Familiarity with Earth science and geospatial datasets is not required, but will be useful.

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.

Prior Python Programming Level of Knowledge Expected

The attendees are expected to have a basic understanding of Python and Jupyter Notebook. Familiarity with Earth science and geospatial datasets is not necessary, but it will be helpful.

G4G Summit Tickets

We offer two Geo for Good (G4G) Summit tickets to attendees of the EE SciPy workshop. The G4G Summit will be held on October 11-12, 2023, at the Google campus in Mountain View, CA. The link to express your interest in the tickets will be provided during the workshop.

Introduction to Earth Engine and geemap (15 mins)

Earth Engine is free for noncommercial and research use. For more than a decade, Earth Engine has enabled planetary-scale Earth data science and analysis by nonprofit organizations, research scientists, and other impact users.

With the launch of Earth Engine for commercial use, commercial customers will be charged for Earth Engine services. However, Earth Engine will remain free of charge for noncommercial use and research projects. Nonprofit organizations, academic institutions, educators, news media, Indigenous governments, and government researchers are eligible to use Earth Engine free of charge, just as they have done for over a decade.

The geemap Python package is built upon the Earth Engine Python API and open-source mapping libraries. It allows Earth Engine users to interactively manipulate, analyze, and visualize geospatial big data in a Jupyter environment. Since its creation in April 2020, geemap has received over 2,700 GitHub stars and is being used by over 900 projects on GitHub.

Google Colab and Earth Engine Python API authentication (5 mins)

image

Install geemap

Uncomment the following line to install geemap if you are running this notebook in Google Colab.

# %pip install geemap[workshop]

Import libraries

import ee
import geemap

Authenticate and initialize Earth Engine

You will need to create a Google Cloud Project and enable the Earth Engine API for the project. You can find detailed instructions here.

geemap.ee_initialize()

Creating interactive maps

Let's create an interactive map using the ipyleaflet plotting backend. The geemap.Map class inherits the ipyleaflet.Map class. Therefore, you can use the same syntax to create an interactive map as you would with ipyleaflet.Map.

Map = geemap.Map()

To display it in a Jupyter notebook, simply ask for the object representation:

Map

To customize the map, you can specify various keyword arguments, such as center ([lat, lon]), zoom, width, and height. The default width is 100%, which takes up the entire cell width of the Jupyter notebook. The height argument accepts a number or a string. If a number is provided, it represents the height of the map in pixels. If a string is provided, the string must be in the format of a number followed by px, e.g., 600px.

Map = geemap.Map(center=[40, -100], zoom=4, height=600)
Map

To hide a control, set control_name to False, e.g., draw_ctrl=False.

Map = geemap.Map(data_ctrl=False, toolbar_ctrl=False, draw_ctrl=False)
Map

Adding basemaps

There are several ways to add basemaps to a map. You can specify the basemap to use in the basemap keyword argument when creating the map. Alternatively, you can add basemap layers to the map using the add_basemap method. Geemap has hundreds of built-in basemaps available that can be easily added to the map with only one line of code.

Create a map by specifying the basemap to use as follows. For example, the Esri.WorldImagery basemap represents the Esri world imagery basemap.

Map = geemap.Map(basemap="Esri.WorldImagery")
Map

You can add as many basemaps as you like to the map. For example, the following code adds the OpenTopoMap basemap to the map above:

Map.add_basemap("OpenTopoMap")

Print out the first 10 basemaps:

basemaps = list(geemap.basemaps.keys())
len(geemap.basemaps)

basemaps[:10]

Using Earth Engine data (30 mins)

Earth Engine data types (Image, ImageCollection, Geometry, Feature, FeatureCollection)

Earth Engine objects are server-side objects rather than client-side objects, which means that they are not stored locally on your computer. Similar to video streaming services (e.g., YouTube, Netflix, and Hulu), which store videos/movies on their servers, Earth Engine data are stored on the Earth Engine servers. We can stream geospatial data from Earth Engine on-the-fly without having to download the data just like we can watch videos from streaming services using a web browser without having to download the entire video to your computer.

  • Image: the fundamental raster data type in Earth Engine.

  • ImageCollection: a stack or time-series of images.

  • Geometry: the fundamental vector data type in Earth Engine.

  • Feature: a Geometry with attributes.

  • FeatureCollection: a set of features.

Image

Raster data in Earth Engine are represented as Image objects. Images are composed of one or more bands and each band has its own name, data type, scale, mask and projection. Each image has metadata stored as a set of properties.

Loading Earth Engine images

image = ee.Image("USGS/SRTMGL1_003")
image

Visualizing Earth Engine images

Map = geemap.Map(center=[21.79, 70.87], zoom=3)
image = ee.Image("USGS/SRTMGL1_003")
vis_params = {
    "min": 0,
    "max": 6000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
Map.addLayer(image, vis_params, "SRTM")
Map

ImageCollection

An ImageCollection is a stack or sequence of images. An ImageCollection can be loaded by passing an Earth Engine asset ID into the ImageCollection constructor. You can find ImageCollection IDs in the Earth Engine Data Catalog.

Loading image collections

For example, to load the image collection of the Sentinel-2 surface reflectance:

collection = ee.ImageCollection("COPERNICUS/S2_SR")

Visualizing image collections

To visualize an Earth Engine ImageCollection, we need to convert an ImageCollection to an Image by compositing all the images in the collection to a single image representing, for example, the min, max, median, mean or standard deviation of the images. For example, to create a median value image from a collection, use the collection.median() method. Let's create a median image from the Sentinel-2 surface reflectance collection:

Map = geemap.Map()
collection = ee.ImageCollection("COPERNICUS/S2_SR")
image = collection.median()

vis = {
    "min": 0.0,
    "max": 3000,
    "bands": ["B4", "B3", "B2"],
}

Map.setCenter(83.277, 17.7009, 12)
Map.addLayer(image, vis, "Sentinel-2")
Map

Filtering image collections

Map = geemap.Map()
collection = (
    ee.ImageCollection("COPERNICUS/S2_SR")
    .filterDate("2021-01-01", "2022-01-01")
    .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", 5))
)
image = collection.median()

vis = {
    "min": 0.0,
    "max": 3000,
    "bands": ["B4", "B3", "B2"],
}

Map.setCenter(83.277, 17.7009, 12)
Map.addLayer(image, vis, "Sentinel-2")
Map

FeatureCollection

A FeatureCollection is a collection of Features. A FeatureCollection is analogous to a GeoJSON FeatureCollection object, i.e., a collection of features with associated properties/attributes. Data contained in a shapefile can be represented as a FeatureCollection.

Loading feature collections

The Earth Engine Data Catalog hosts a variety of vector datasets (e.g,, US Census data, country boundaries, and more) as feature collections. You can find feature collection IDs by searching the data catalog. For example, to load the TIGER roads data by the U.S. Census Bureau:

Map = geemap.Map()
fc = ee.FeatureCollection("TIGER/2016/Roads")
Map.setCenter(-73.9596, 40.7688, 12)
Map.addLayer(fc, {}, "Census roads")
Map

Filtering feature collections

Map = geemap.Map()
states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.eq("NAME", "Louisiana"))
Map.addLayer(fc, {}, "Louisiana")
Map.centerObject(fc, 7)
Map

feat = fc.first()
feat.toDictionary()

Map = geemap.Map()
states = ee.FeatureCollection("TIGER/2018/States")
fc = states.filter(ee.Filter.inList("NAME", ["California", "Oregon", "Washington"]))
Map.addLayer(fc, {}, "West Coast")
Map.centerObject(fc, 5)
Map

region = Map.user_roi
if region is None:
    region = ee.Geometry.BBox(-88.40, 29.88, -77.90, 35.39)

fc = ee.FeatureCollection("TIGER/2018/States").filterBounds(region)
Map.addLayer(fc, {}, "Southeastern U.S.")
Map.centerObject(fc, 6)

Visualizing feature collections

Map = geemap.Map(center=[40, -100], zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
Map.addLayer(states, {}, "US States")
Map

Map = geemap.Map(center=[40, -100], zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
style = {"color": "0000ffff", "width": 2, "lineType": "solid", "fillColor": "FF000080"}
Map.addLayer(states.style(**style), {}, "US States")
Map

Map = geemap.Map(center=[40, -100], zoom=4)
states = ee.FeatureCollection("TIGER/2018/States")
vis_params = {
    "color": "000000",
    "colorOpacity": 1,
    "pointSize": 3,
    "pointShape": "circle",
    "width": 2,
    "lineType": "solid",
    "fillColorOpacity": 0.66,
}
palette = ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"]
Map.add_styled_vector(
    states, column="NAME", palette=palette, layer_name="Styled vector", **vis_params
)
Map

Earth Engine Data Catalog

The Earth Engine Data Catalog hosts a variety of geospatial datasets. As of March 2023, the catalog contains over 1,000 datasets with a total size of over 80 petabytes. Some notable datasets include: Landsat, Sentinel, MODIS, NAIP, etc. For a complete list of datasets in CSV or JSON formats, see the Earth Engine Datasets List.

Searching for datasets

The Earth Engine Data Catalog is searchable. You can search datasets by name, keyword, or tag. For example, enter "elevation" in the search box will filter the catalog to show only datasets containing "elevation" in their name, description, or tags. 52 datasets are returned for this search query. Scroll down the list to find the NASA SRTM Digital Elevation 30m dataset. On each dataset page, you can find the following information, including Dataset Availability, Dataset Provider, Earth Engine Snippet, Tags, Description, Code Example, and more (see {numref}ch03_gee_srtm). One important piece of information is the Image/ImageCollection/FeatureCollection ID of each dataset, which is essential for accessing the dataset through the Earth Engine JavaScript or Python APIs.

Map = geemap.Map()
Map

dataset_xyz = ee.Image("USGS/SRTMGL1_003")
Map.addLayer(dataset_xyz, {}, "USGS/SRTMGL1_003")

Map = geemap.Map()
dem = ee.Image("USGS/SRTMGL1_003")
vis_params = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
Map.addLayer(dem, vis_params, "SRTM DEM")
Map

Using the datasets module

from geemap.datasets import DATA

Map = geemap.Map(center=[40, -100], zoom=4)
dataset = ee.Image(DATA.USGS_GAP_CONUS_2011)
Map.addLayer(dataset, {}, "GAP CONUS")
Map

from geemap.datasets import get_metadata

get_metadata(DATA.USGS_GAP_CONUS_2011)

Converting Earth Engine JavaScripts to Python

Find some Earth Engine JavaScript code that you want to convert to Python. For example, you can grab some sample code from the Earth Engine Documentation.

Map = geemap.Map()
Map

# Load an image.
image = ee.Image("LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318")

# Define the visualization parameters.
vizParams = {"bands": ["B5", "B4", "B3"], "min": 0, "max": 0.5, "gamma": [0.95, 1.1, 1]}

# Center the map and display the image.
Map.setCenter(-122.1899, 37.5010, 10)
# San Francisco Bay
Map.addLayer(image, vizParams, "False color composite")

Exercise 1 - Creating cloud-free imagery

Create a cloud-free imagery of Texas for the year of 2022. You can use either Landsat 9 or Sentinel-2 imagery. Relevant Earth Engine assets:

Break 1 (10 mins)

Visualizing Earth Engine data (30 mins)

Geemap Inspector tool, plotting tool, interactive GUI for data visualization

Using the inspector tool

Map = geemap.Map(center=(40, -100), zoom=4)

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

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

Map.addLayer(dem, vis_params, "SRTM DEM")
Map.addLayer(
    landsat7,
    {"bands": ["B4", "B3", "B2"], "min": 20, "max": 200, "gamma": 2.0},
    "Landsat 7",
)
Map.addLayer(states, {}, "US States")
Map.add_inspector()
Map

Using the plotting tool

Map = geemap.Map(center=[40, -100], zoom=4)

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

landsat_vis = {"bands": ["B4", "B3", "B2"], "gamma": 1.4}
Map.addLayer(landsat7, landsat_vis, "Landsat")

hyperion = ee.ImageCollection("EO1/HYPERION").filter(
    ee.Filter.date("2016-01-01", "2017-03-01")
)

hyperion_vis = {
    "min": 1000.0,
    "max": 14000.0,
    "gamma": 2.5,
}
Map.addLayer(hyperion, hyperion_vis, "Hyperion")
Map.add_plot_gui()
Map

Map.set_plot_options(add_marker_cluster=True, overlay=True)

Legends, color bars, and labels

Built-in legends

from geemap.legends import builtin_legends

for legend in builtin_legends:
    print(legend)

Add NLCD WMS layer and legend to the map.

Map = geemap.Map(center=[40, -100], zoom=4)
Map.add_basemap("Google Hybrid")
Map.add_basemap("NLCD 2019 CONUS Land Cover")
Map.add_legend(builtin_legend="NLCD", max_width="100px")
Map

Add NLCD Earth Engine layer and legend to the map.

Map = geemap.Map(center=[40, -100], zoom=4)
Map.add_basemap("HYBRID")

nlcd = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019")
landcover = nlcd.select("landcover")

Map.addLayer(landcover, {}, "NLCD Land Cover 2019")
Map.add_legend(
    title="NLCD Land Cover Classification", builtin_legend="NLCD", height="460px"
)
Map

Custom legends

Add a custom legend by specifying the colors and labels.

Map = geemap.Map(add_google_map=False)

labels = ["One", "Two", "Three", "Four", "etc"]

# colors can be defined using either hex code or RGB (0-255, 0-255, 0-255)
colors = ["#8DD3C7", "#FFFFB3", "#BEBADA", "#FB8072", "#80B1D3"]
# legend_colors = [(255, 0, 0), (127, 255, 0), (127, 18, 25), (36, 70, 180), (96, 68 123)]

Map.add_legend(labels=labels, colors=colors, position="bottomright")
Map

Add a custom legend by specifying a dictionary of colors and labels.

Map = geemap.Map(center=[40, -100], zoom=4)
Map.add_basemap("Google Hybrid")

legend_dict = {
    "11 Open Water": "466b9f",
    "12 Perennial Ice/Snow": "d1def8",
    "21 Developed, Open Space": "dec5c5",
    "22 Developed, Low Intensity": "d99282",
    "23 Developed, Medium Intensity": "eb0000",
    "24 Developed High Intensity": "ab0000",
    "31 Barren Land (Rock/Sand/Clay)": "b3ac9f",
    "41 Deciduous Forest": "68ab5f",
    "42 Evergreen Forest": "1c5f2c",
    "43 Mixed Forest": "b5c58f",
    "51 Dwarf Scrub": "af963c",
    "52 Shrub/Scrub": "ccb879",
    "71 Grassland/Herbaceous": "dfdfc2",
    "72 Sedge/Herbaceous": "d1d182",
    "73 Lichens": "a3cc51",
    "74 Moss": "82ba9e",
    "81 Pasture/Hay": "dcd939",
    "82 Cultivated Crops": "ab6c28",
    "90 Woody Wetlands": "b8d9eb",
    "95 Emergent Herbaceous Wetlands": "6c9fb8",
}

nlcd = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019")
landcover = nlcd.select("landcover")

Map.addLayer(landcover, {}, "NLCD Land Cover 2019")
Map.add_legend(title="NLCD Land Cover Classification", legend_dict=legend_dict)
Map

Creating color bars

Add a horizontal color bar.

Map = geemap.Map()

dem = ee.Image("USGS/SRTMGL1_003")
vis_params = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}

Map.addLayer(dem, vis_params, "SRTM DEM")
Map.add_colorbar(vis_params, label="Elevation (m)", layer_name="SRTM DEM")
Map

Add a vertical color bar.

Map.add_colorbar(
    vis_params,
    label="Elevation (m)",
    layer_name="SRTM DEM",
    orientation="vertical",
    max_width="100px",
)

Split-panel map and linked maps

Split-panel maps

Create a split map with basemaps.

Map = geemap.Map()
Map.split_map(left_layer="Google Terrain", right_layer="OpenTopoMap")
Map

Create a split map with Earth Engine layers.

Map = geemap.Map(center=(40, -100), zoom=4, height=600)

nlcd_2001 = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2001").select("landcover")
nlcd_2019 = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019").select("landcover")

left_layer = geemap.ee_tile_layer(nlcd_2001, {}, "NLCD 2001")
right_layer = geemap.ee_tile_layer(nlcd_2019, {}, "NLCD 2019")

Map.split_map(left_layer, right_layer)
Map

Linked maps

Create a 2x2 linked map for visualizing Sentinel-2 imagery with different band combinations.

image = (
    ee.ImageCollection("COPERNICUS/S2")
    .filterDate("2018-09-01", "2018-09-30")
    .map(lambda img: img.divide(10000))
    .median()
)

vis_params = [
    {"bands": ["B4", "B3", "B2"], "min": 0, "max": 0.3, "gamma": 1.3},
    {"bands": ["B8", "B11", "B4"], "min": 0, "max": 0.3, "gamma": 1.3},
    {"bands": ["B8", "B4", "B3"], "min": 0, "max": 0.3, "gamma": 1.3},
    {"bands": ["B12", "B12", "B4"], "min": 0, "max": 0.3, "gamma": 1.3},
]

labels = [
    "Natural Color (B4/B3/B2)",
    "Land/Water (B8/B11/B4)",
    "Color Infrared (B8/B4/B3)",
    "Vegetation (B12/B11/B4)",
]

geemap.linked_maps(
    rows=2,
    cols=2,
    height="300px",
    center=[38.4151, 21.2712],
    zoom=12,
    ee_objects=[image],
    vis_params=vis_params,
    labels=labels,
    label_position="topright",
)

Timeseries inspector and time slider

Timeseries inspector

Check the available years of NLCD.

Map = geemap.Map(center=[40, -100], zoom=4)
collection = ee.ImageCollection("USGS/NLCD_RELEASES/2019_REL/NLCD").select("landcover")
vis_params = {"bands": ["landcover"]}
years = collection.aggregate_array("system:index").getInfo()
years

Create a timeseries inspector for NLCD.

Map.ts_inspector(
    left_ts=collection,
    right_ts=collection,
    left_names=years,
    right_names=years,
    left_vis=vis_params,
    right_vis=vis_params,
    width="80px",
)
Map

Time slider

Create a map for visualizing MODIS vegetation data.

Map = geemap.Map()

collection = (
    ee.ImageCollection("MODIS/MCD43A4_006_NDVI")
    .filter(ee.Filter.date("2018-06-01", "2018-07-01"))
    .select("NDVI")
)
vis_params = {
    "min": 0.0,
    "max": 1.0,
    "palette": "ndvi",
}

Map.add_time_slider(collection, vis_params, time_interval=2)
Map

Create a map for visualizing weather data.

Map = geemap.Map()

collection = (
    ee.ImageCollection("NOAA/GFS0P25")
    .filterDate("2018-12-22", "2018-12-23")
    .limit(24)
    .select("temperature_2m_above_ground")
)

vis_params = {
    "min": -40.0,
    "max": 35.0,
    "palette": ["blue", "purple", "cyan", "green", "yellow", "red"],
}

labels = [str(n).zfill(2) + ":00" for n in range(0, 24)]
Map.add_time_slider(collection, vis_params, labels=labels, time_interval=1, opacity=0.8)
Map

Visualizing Sentinel-2 imagery

Map = geemap.Map(center=[37.75, -122.45], zoom=12)

collection = (
    ee.ImageCollection("COPERNICUS/S2_SR")
    .filterBounds(ee.Geometry.Point([-122.45, 37.75]))
    .filterMetadata("CLOUDY_PIXEL_PERCENTAGE", "less_than", 10)
)

vis_params = {"min": 0, "max": 4000, "bands": ["B8", "B4", "B3"]}

Map.add_time_slider(collection, vis_params)
Map

Exercise 2 - Creating land cover maps with a legend

Create a split map for visualizing NLCD land cover change in Texas between 2001 and 2019. Add the NLCD legend to the map. Relevant Earth Engine assets:

Analyzing Earth Engine data (30 mins)

Image descriptive statistics

Use a sample Landsat image.

Map = geemap.Map()

centroid = ee.Geometry.Point([-122.4439, 37.7538])
image = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").filterBounds(centroid).first()
vis = {"min": 0, "max": 3000, "bands": ["B5", "B4", "B3"]}

Map.centerObject(centroid, 8)
Map.addLayer(image, vis, "Landsat-8")
Map.addLayer(centroid, {}, "Centroid")
Map

Check image properties.

image.propertyNames()

Check image property values.

image.toDictionary()

Get specific image properties.

image.get("CLOUD_COVER")  # 0.05

Get image properties with easy-to-read time format.

props = geemap.image_props(image)
props

Compute image descriptive statistics.

stats = geemap.image_stats(image, scale=30)
stats

Zonal statistics

Zonal statistics

Add Earth Engine data to the map.

Map = geemap.Map(center=[40, -100], zoom=4)

# Add NASA SRTM
dem = ee.Image("USGS/SRTMGL1_003")
dem_vis = {
    "min": 0,
    "max": 4000,
    "palette": ["006633", "E5FFCC", "662A00", "D8D8D8", "F5F5F5"],
}
Map.addLayer(dem, dem_vis, "SRTM DEM")

# Add 5-year Landsat TOA composite
landsat = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
landsat_vis = {"bands": ["B4", "B3", "B2"], "gamma": 1.4}
Map.addLayer(landsat, landsat_vis, "Landsat", False)

# Add US Census States
states = ee.FeatureCollection("TIGER/2018/States")
style = {"fillColor": "00000000"}
Map.addLayer(states.style(**style), {}, "US States")
Map

Compute zonal statistics. In this case, we compute the mean elevation of each state and export the results to a CSV file.

out_dem_stats = "dem_stats.csv"
geemap.zonal_stats(
    dem, states, out_dem_stats, stat_type="MEAN", scale=1000, return_fc=False
)

Display the csv file as a table.

geemap.csv_to_df(out_dem_stats).sort_values(by=["mean"], ascending=True)

Compute zonal statistics of mean spectral values of each state.

out_landsat_stats = "landsat_stats.csv"
geemap.zonal_stats(
    landsat,
    states,
    out_landsat_stats,
    stat_type="MEAN",
    scale=1000,
    return_fc=False,
)

geemap.csv_to_df(out_landsat_stats)

Zonal statistics by group

Compute zonal statistics by group. In this case, we compute the area of each land cover type in each state and export the results to a CSV file.

Map = geemap.Map(center=[40, -100], zoom=4)

# Add NLCD data
dataset = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019")
landcover = dataset.select("landcover")
Map.addLayer(landcover, {}, "NLCD 2019")

# Add US census states
states = ee.FeatureCollection("TIGER/2018/States")
style = {"fillColor": "00000000"}
Map.addLayer(states.style(**style), {}, "US States")

# Add NLCD legend
Map.add_legend(title="NLCD Land Cover", builtin_legend="NLCD")
Map

Compute zonal statistics by group and convert the area unit from square meters to square kilometers.

nlcd_stats = "nlcd_stats.csv"

geemap.zonal_stats_by_group(
    landcover,
    states,
    nlcd_stats,
    stat_type="SUM",
    denominator=1e6,
    decimal_places=2,
)

geemap.csv_to_df(nlcd_stats)

Calculate the percentage of each land cover type in each state.

nlcd_stats = "nlcd_stats_pct.csv"

geemap.zonal_stats_by_group(
    landcover,
    states,
    nlcd_stats,
    stat_type="PERCENTAGE",
    denominator=1e6,
    decimal_places=2,
)

Zonal statistics with two images

geemap.csv_to_df(nlcd_stats)

Zonal statistics by zone

The zonal statistics by zone algorithm is similar to the zonal statistics by group algorithm, but it takes an image as the zone input instead of a feature collection.

Map = geemap.Map(center=[40, -100], zoom=4)
dem = ee.Image("USGS/3DEP/10m")
vis = {"min": 0, "max": 4000, "palette": "terrain"}
Map.addLayer(dem, vis, "DEM")
Map

landcover = ee.Image("USGS/NLCD_RELEASES/2019_REL/NLCD/2019").select("landcover")
Map.addLayer(landcover, {}, "NLCD 2019")
Map.add_legend(title="NLCD Land Cover Classification", builtin_legend="NLCD")

Computer the mean elevation of each land cover type.

stats = geemap.image_stats_by_zone(dem, landcover, reducer="MEAN")
stats

stats.to_csv("mean.csv", index=False)

Compute the standard deviation of each land cover type.

geemap.image_stats_by_zone(dem, landcover, out_csv="std.csv", reducer="STD")

geemap.csv_to_df("std.csv")

Exercise 3 - Zonal statistics

Find out which state has the highest mean temperature on in the United States on June 28, 2023. Relevant Earth Engine assets:

Coordinate grids and fishnets

Creating coordinate grids

Create a latitudinal grid with a 5-degree interval.

lat_grid = geemap.latitude_grid(step=5.0, west=-180, east=180, south=-85, north=85)

Map = geemap.Map()
style = {"fillColor": "00000000"}
Map.addLayer(lat_grid.style(**style), {}, "Latitude Grid")
Map

df = geemap.ee_to_df(lat_grid)
df.head()

Create a longitudinal grid with a 5-degree interval.

lon_grid = geemap.longitude_grid(step=5.0, west=-180, east=180, south=-85, north=85)

Map = geemap.Map()
style = {"fillColor": "00000000"}
Map.addLayer(lon_grid.style(**style), {}, "Longitude Grid")
Map

Create a rectangular grid with a 10-degree interval.

grid = geemap.latlon_grid(
    lat_step=10, lon_step=10, west=-180, east=180, south=-85, north=85
)

Map = geemap.Map()
style = {"fillColor": "00000000"}
Map.addLayer(grid.style(**style), {}, "Coordinate Grid")
Map

Creating fishnets

Create a fishnet based on an Earth Engine geometry.

Map = geemap.Map()
Map

Use the drawing tools to draw a polygon on the map above. If no polygon is drawn, the default polygon will be used.

roi = Map.user_roi

if roi is None:
    roi = ee.Geometry.BBox(-112.8089, 33.7306, -88.5951, 46.6244)
    Map.addLayer(roi, {}, "ROI")
    Map.user_roi = None

Map.centerObject(roi)

Create a fishnet based on a user-drawn polygon with a 2-degree interval.

fishnet = geemap.fishnet(roi, h_interval=2.0, v_interval=2.0, delta=1)
style = {"color": "blue", "fillColor": "00000000"}
Map.addLayer(fishnet.style(**style), {}, "Fishnet")

Create a new map.

Map = geemap.Map()
Map

Draw a polygon on the map.

roi = Map.user_roi

if roi is None:
    roi = ee.Geometry.Polygon(
        [
            [
                [-64.602356, -1.127399],
                [-68.821106, -12.625598],
                [-60.647278, -22.498601],
                [-47.815247, -21.111406],
                [-43.860168, -8.913564],
                [-54.582825, -0.775886],
                [-60.823059, 0.454555],
                [-64.602356, -1.127399],
            ]
        ]
    )
    Map.addLayer(roi, {}, "ROI")

Map.centerObject(roi)
Map

Create a fishnet based on a user-drawn polygon with specified number of rows and columns.

fishnet = geemap.fishnet(roi, rows=6, cols=8, delta=1)
style = {"color": "blue", "fillColor": "00000000"}
Map.addLayer(fishnet.style(**style), {}, "Fishnet")

Land use and land cover change analysis

Forest cover mapping

We will use the Hansen Global Forest Change v1.10 (2000-2022) dataset.

dataset = ee.Image("UMD/hansen/global_forest_change_2022_v1_10")
dataset.bandNames()

Select the imagery for 2000.

Map = geemap.Map()
first_bands = ["first_b50", "first_b40", "first_b30"]
first_image = dataset.select(first_bands)
Map.addLayer(first_image, {"bands": first_bands, "gamma": 1.5}, "Year 2000 Bands 5/4/3")

Select the imagery for 2022.

last_bands = ["last_b50", "last_b40", "last_b30"]
last_image = dataset.select(last_bands)
Map.addLayer(last_image, {"bands": last_bands, "gamma": 1.5}, "Year 2022 Bands 5/4/3")

Select the tree cover imagery for 2000.

treecover = dataset.select(["treecover2000"])
treeCoverVisParam = {"min": 0, "max": 100, "palette": ["black", "green"]}
name = "Tree cover (%)"
Map.addLayer(treecover, treeCoverVisParam, name)
Map.add_colorbar(treeCoverVisParam, label=name, layer_name=name)
Map.add_layer_manager()
Map

Extract tree cover 2000 by using the threshold of 10%.

threshold = 10
treecover_bin = treecover.gte(threshold).selfMask()
treeVisParam = {"palette": ["green"]}
Map.addLayer(treecover_bin, treeVisParam, "Tree cover bin")

Forest loss and gain mapping

Visualize forest loss.

Map = geemap.Map()
Map.add_basemap("Google Hybrid")
treeloss_year = dataset.select(["lossyear"])
treeLossVisParam = {"min": 0, "max": 22, "palette": ["yellow", "red"]}
layer_name = "Tree loss year since 2000"
Map.addLayer(treeloss_year, treeLossVisParam, layer_name)
Map.add_colorbar(treeLossVisParam, label=layer_name, layer_name=layer_name)
Map.add_layer_manager()
Map

Compare forest loss and gain.

Map = geemap.Map()
Map.add_basemap("Google Hybrid")
treeloss = dataset.select(["loss"]).selfMask()
treegain = dataset.select(["gain"]).selfMask()
Map.addLayer(treeloss, {"palette": "red"}, "Tree loss")
Map.addLayer(treegain, {"palette": "yellow"}, "Tree gain")
Map.add_layer_manager()
Map

Zonal statistics by country

Compute zonal statistics to find out which country has the largest forest area in 2000.

Add a country boundary layer to the map.

Map = geemap.Map()
countries = ee.FeatureCollection(geemap.examples.get_ee_path("countries"))
style = {"color": "#000000ff", "fillColor": "#00000000"}
Map.addLayer(countries.style(**style), {}, "Countries")
Map

Compute zonal statistics by country.

geemap.zonal_stats(
    treecover_bin,
    countries,
    "forest_cover.csv",
    stat_type="SUM",
    denominator=1e6,
    scale=1000,
)

Create a pie chart to visualize the forest area by country.

geemap.pie_chart(
    "forest_cover.csv", names="NAME", values="sum", max_rows=20, height=400
)

Create a bar chart to visualize the forest area by country.

geemap.bar_chart(
    "forest_cover.csv",
    x="NAME",
    y="sum",
    max_rows=20,
    x_label="Country",
    y_label="Forest area (km2)",
)

Calculate the forest loss area by country.

geemap.zonal_stats(
    treeloss.gt(0),
    countries,
    "treeloss.csv",
    stat_type="SUM",
    denominator=1e6,
    scale=1000,
)

Create a bar chart to visualize the forest loss area by country.

geemap.pie_chart("treeloss.csv", names="NAME", values="sum", max_rows=20, height=600)

Create a bar chart to visualize the forest loss area by country.

geemap.bar_chart(
    "treeloss.csv",
    x="NAME",
    y="sum",
    max_rows=20,
    x_label="Country",
    y_label="Forest loss area (km2)",
)

Exercise 4 - Analyzing forest cover gain and loss

Find out which US state has the largest forest gain and loss between 2000 and 2022. Create pie charts and bar charts to show the results. Relevant Earth Engine assets:

Break 2 (10 mins)

Exporting Earth Engine data (30 mins)

Exporting images

Add a Landsat image to the map.

Map = geemap.Map()

image = ee.Image("LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318").select(
    ["B5", "B4", "B3"]
)

vis_params = {"min": 0, "max": 0.5, "gamma": [0.95, 1.1, 1]}

Map.centerObject(image)
Map.addLayer(image, vis_params, "Landsat")
Map

Add a rectangle to the map.

region = ee.Geometry.BBox(-122.5955, 37.5339, -122.0982, 37.8252)
fc = ee.FeatureCollection(region)
style = {"color": "ffff00ff", "fillColor": "00000000"}
Map.addLayer(fc.style(**style), {}, "ROI")

To local drive

geemap.ee_export_image(image, filename="landsat.tif", scale=30, region=region)

Check image projection.

projection = image.select(0).projection().getInfo()
projection

crs = projection["crs"]
crs_transform = projection["transform"]

Specify region, crs, and crs_transform.

geemap.ee_export_image(
    image,
    filename="landsat_crs.tif",
    crs=crs,
    crs_transform=crs_transform,
    region=region,
)

To Google Drive

geemap.ee_export_image_to_drive(
    image, description="landsat", folder="export", region=region, scale=30
)

Exporting image collections

point = ee.Geometry.Point(-99.2222, 46.7816)
collection = (
    ee.ImageCollection("USDA/NAIP/DOQQ")
    .filterBounds(point)
    .filterDate("2008-01-01", "2018-01-01")
    .filter(ee.Filter.listContains("system:band_names", "N"))
)

collection.aggregate_array("system:index")

To local drive

geemap.ee_export_image_collection(collection, out_dir=".", scale=10)

To Google Drive

geemap.ee_export_image_collection_to_drive(collection, folder="export", scale=10)

Exporting feature collections

Map = geemap.Map()
fc = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017").filter(
    ee.Filter.eq("country_na", "Germany")
)
Map.addLayer(fc, {}, "Germany")
Map.centerObject(fc)
Map

To local drive

geemap.ee_to_shp(fc, filename="Germany.shp", selectors=None)

geemap.ee_export_vector(fc, filename="Germany2.shp")

geemap.ee_to_geojson(fc, filename="Germany.geojson")

geemap.ee_to_csv(fc, filename="Germany.csv")

gdf = geemap.ee_to_gdf(fc)
gdf

df = geemap.ee_to_df(fc)
df

To Google Drive

geemap.ee_export_vector_to_drive(
    fc, description="germany", fileFormat="SHP", folder="export"
)

Exercise 5 - Exporting images by a fishnet

Create a fishnet with a 4-degree interval based on the extent of [-112.5439, 34.0891, -85.0342, 49.6858]. Use the fishnet to download the Landsat 7 image tiles by the fishnet using the geemap.download_ee_image_tiles() and geemap.download_ee_image_tiles_parallel() functions. Relevant Earth Engine assets:

  • ee.Image('LANDSAT/LE7_TOA_5YEAR/1999_2003')

Map = geemap.Map(center=[40, -100], zoom=4)
image = ee.Image("LANDSAT/LE7_TOA_5YEAR/1999_2003").select(["B4", "B3", "B2"])
Map.addLayer(image, {"min": 20, "max": 200, "gamma": 2.0}, "Landsat")
region = ee.Geometry.BBox(-112.5439, 34.0891, -85.0342, 49.6858)
Map.addLayer(region, {}, "ROI")
Map.centerObject(region)
Map

fishnet = geemap.fishnet(region, h_interval=4.0, v_interval=4.0, delta=0.5)
style = {"color": "ffff00ff", "fillColor": "00000000"}
Map.addLayer(fishnet.style(**style), {}, "Fishnet")

geemap.download_ee_image_tiles(
    image, fishnet, out_dir="tiles", scale=1000, crs="EPSG:3857"
)

geemap.download_ee_image_tiles_parallel(
    image, fishnet, out_dir="tiles", scale=1000, crs="EPSG:3857"
)

Creating Satellite timelapse animations (30 mins)

Creating satellite timeseries

Use the Harmonized Sentinel-2 MSI dataset to create a timeseries for a given location.

collection = ee.ImageCollection("COPERNICUS/S2_HARMONIZED").filterMetadata(
    "CLOUDY_PIXEL_PERCENTAGE", "less_than", 10
)

Specify the location of interest and date range.

start_date = "2016-01-01"
end_date = "2022-12-31"
region = ee.Geometry.BBox(-122.5549, 37.6968, -122.3446, 37.8111)

Create an annual composite.

images = geemap.create_timeseries(
    collection, start_date, end_date, region, frequency="year", reducer="median"
)
images

Display the timeseries.

Map = geemap.Map()

vis_params = {"min": 0, "max": 4000, "bands": ["B8", "B4", "B3"]}
labels = [str(y) for y in range(2016, 2023)]

Map.addLayer(images, vis_params, "Sentinel-2", False)
Map.add_time_slider(images, vis_params, time_interval=2, labels=labels)
Map.centerObject(region)
Map

Creating satellite timelapse animations

NAIP timelapse

Map = geemap.Map(center=[40, -100], zoom=4)
Map

roi = Map.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-99.1019, 47.1274, -99.0334, 47.1562)
    Map.addLayer(roi, {}, "ROI")
    Map.centerObject(roi)

collection = geemap.naip_timeseries(roi, start_year=2009, end_year=2022, RGBN=True)

years = geemap.image_dates(collection, date_format="YYYY").getInfo()
print(years)

size = len(years)
images = collection.toList(size)
for i in range(size):
    image = ee.Image(images.get(i))
    Map.addLayer(image, {"bands": ["N", "R", "G"]}, years[i])
Map

timelapse = geemap.naip_timelapse(
    roi,
    out_gif="naip.gif",
    bands=["N", "R", "G"],
    frames_per_second=3,
    title="NAIP Timelapse",
)
geemap.show_image(timelapse)

Landsat timelapse

Map = geemap.Map()
Map

Pan and zoom the map to an area of interest. Use the drawing tools to draw a rectangle on the map. If no rectangle is drawn, the default rectangle shown below will be used.

roi = Map.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-74.7222, -8.5867, -74.1596, -8.2824)
    Map.addLayer(roi)
    Map.centerObject(roi)

timelapse = geemap.landsat_timelapse(
    roi,
    out_gif="landsat.gif",
    start_year=1984,
    end_year=2022,
    start_date="01-01",
    end_date="12-31",
    bands=["SWIR1", "NIR", "Red"],
    frames_per_second=5,
    title="Landsat Timelapse",
    progress_bar_color="blue",
    mp4=True,
)
geemap.show_image(timelapse)

Map = geemap.Map()
roi = ee.Geometry.BBox(-115.5541, 35.8044, -113.9035, 36.5581)
Map.addLayer(roi)
Map.centerObject(roi)
Map

timelapse = geemap.landsat_timelapse(
    roi,
    out_gif="las_vegas.gif",
    start_year=1984,
    end_year=2022,
    bands=["NIR", "Red", "Green"],
    frames_per_second=5,
    title="Las Vegas, NV",
    font_color="blue",
)
geemap.show_image(timelapse)

Map = geemap.Map()
roi = ee.Geometry.BBox(113.8252, 22.1988, 114.0851, 22.3497)
Map.addLayer(roi)
Map.centerObject(roi)
Map

timelapse = geemap.landsat_timelapse(
    roi,
    out_gif="hong_kong.gif",
    start_year=1990,
    end_year=2022,
    start_date="01-01",
    end_date="12-31",
    bands=["SWIR1", "NIR", "Red"],
    frames_per_second=3,
    title="Hong Kong",
)
geemap.show_image(timelapse)

Sentinel-2 timelapse

Map = geemap.Map()
Map

Pan and zoom the map to an area of interest. Use the drawing tools to draw a rectangle on the map. If no rectangle is drawn, the default rectangle shown below will be used.

roi = Map.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-74.7222, -8.5867, -74.1596, -8.2824)
    Map.addLayer(roi)
    Map.centerObject(roi)

timelapse = geemap.sentinel2_timelapse(
    roi,
    out_gif="sentinel2.gif",
    start_year=2016,
    end_year=2021,
    start_date="01-01",
    end_date="12-31",
    frequency="year",
    bands=["SWIR1", "NIR", "Red"],
    frames_per_second=3,
    title="Sentinel-2 Timelapse",
)
geemap.show_image(timelapse)

MODIS timelapse

MODIS vegetation indices

Map = geemap.Map()
Map

roi = Map.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-18.6983, -36.1630, 52.2293, 38.1446)
    Map.addLayer(roi)
    Map.centerObject(roi)

timelapse = geemap.modis_ndvi_timelapse(
    roi,
    out_gif="ndvi.gif",
    data="Terra",
    band="NDVI",
    start_date="2000-01-01",
    end_date="2022-12-31",
    frames_per_second=3,
    title="MODIS NDVI Timelapse",
    overlay_data="countries",
)
geemap.show_image(timelapse)

MODIS temperature

Map = geemap.Map()
Map

roi = Map.user_roi
if roi is None:
    roi = ee.Geometry.BBox(-171.21, -57.13, 177.53, 79.99)
    Map.addLayer(roi)
    Map.centerObject(roi)

timelapse = geemap.modis_ocean_color_timelapse(
    satellite="Aqua",
    start_date="2018-01-01",
    end_date="2020-12-31",
    roi=roi,
    frequency="month",
    out_gif="temperature.gif",
    overlay_data="continents",
    overlay_color="yellow",
    overlay_opacity=0.5,
)
geemap.show_image(timelapse)

GOES timelapse

roi = ee.Geometry.BBox(167.1898, -28.5757, 202.6258, -12.4411)
start_date = "2022-01-15T03:00:00"
end_date = "2022-01-15T07:00:00"
data = "GOES-17"
scan = "full_disk"

timelapse = geemap.goes_timelapse(
    roi, "goes.gif", start_date, end_date, data, scan, framesPerSecond=5
)
geemap.show_image(timelapse)

roi = ee.Geometry.BBox(-159.5954, 24.5178, -114.2438, 60.4088)
start_date = "2021-10-24T14:00:00"
end_date = "2021-10-25T01:00:00"
data = "GOES-17"
scan = "full_disk"

timelapse = geemap.goes_timelapse(
    roi, "hurricane.gif", start_date, end_date, data, scan, framesPerSecond=5
)
geemap.show_image(timelapse)

roi = ee.Geometry.BBox(-121.0034, 36.8488, -117.9052, 39.0490)
start_date = "2020-09-05T15:00:00"
end_date = "2020-09-06T02:00:00"
data = "GOES-17"
scan = "full_disk"

timelapse = geemap.goes_fire_timelapse(
    roi, "fire.gif", start_date, end_date, data, scan, framesPerSecond=5
)
geemap.show_image(timelapse)

Exercise 6 - Creating timelapse animations

Use the geemap timelapse GUI to create a timelapse animation for any location of your choice. Share the timelapse on social media and use the hashtagd such as #EarthEngine and #geemap. See this example.

Map = geemap.Map()
Map.add_gui("timelapse")
Map

Break 3 (10 mins)

Building interactive web apps (30 mins)

This section is optional. We might not have enough time to cover this section.

Follow the instructions here to build an interactive Earth Engine web app with Solara and geemap. You need to sign up for a free Hugging Face account to create the web app. It is free and easy to sign up.

Exercise 7 - Building an interactive web app for visualizing land cover change

After following the instructions above, you should have a web app that looks like this:

The web app URL should look like this: https://giswqs-solara-geemap.hf.space/split-map.